Dark Holograms and Gravitational Waves

Spectra of stochastic gravitational waves (GW) generated in cosmological first-order phase transitions are computed within strongly correlated theories with a dual holographic description. The theories are mostly used as models of dark sectors. In particular, we consider the so-called Witten-Sakai-Sugimoto model, a $SU(N)$ gauge theory coupled to different matter fields in both the fundamental and the adjoint representations. The model has a well-known top-down holographic dual description which allows us to perform reliable calculations in the strongly coupled regime. We consider the GW spectra from bubble collisions and sound waves arising from two different kinds of first-order phase transitions: a confinement/deconfinement one and a chiral symmetry breaking/restoration one. Depending on the model parameters, we find that the GW spectra may fall within the sensibility region of ground-based and space-based interferometers, as well as of Pulsar Timing Arrays. In the latter case, the signal could be compatible with the recent potential observation by NANOGrav. When the two phase transitions happen at different critical temperatures, characteristic spectra with double frequency peaks show up. Moreover, in this case we explicitly show how to correct the redshift factors appearing in the formulae for the GW power spectra to account for the fact that adiabatic expansion from the first transition to the present times cannot be assumed anymore.


Introduction
The measurement of the first direct gravitational wave (GW) signal by LIGO in 2015 [1] has started a new era in observational astrophysics. Not only the observation of black hole and neutron star mergers are tremendously important discoveries, but current and future experiments are now expected to be able to measure GW signals from several different sources. This promises to give experimental access to physics which would be challenging to investigate with other types of observations. Not surprisingly, there are currently several experiments in the developing phase, which will considerably extend the accessible GW frequency and sensitivity ranges in the near future. In this situation, it is of clear interest to study possible sources of GWs which could be detected in these facilities.
In this paper, we consider stochastic GW spectra produced in first-order cosmological phase transitions. The generation of GWs, in this case, is determined by the dynamics of bubbles of true vacuum nucleated in the metastable phase once the temperature of the Universe descends below the phase transition temperature [2][3][4][5][6]. The bubbles can generate GWs either by their collisions or by their interaction with the plasma medium, through sound waves or turbulence. We refer to [7][8][9] for reviews. It is a challenging task to connect the qualitative picture of the bubble dynamics to solid predictions for the power spectra of GWs that can be observed in experimental devices. Luckily, there are general formulae in the literature that estimate the GW spectra once some parameters characterizing the phase transition are known. These parameters depend on the details of the microscopic model describing the transition. The evaluation of the parameters and the formulae for the spectra typically rely on a series of controlled and less controlled approximations. It is a crucial goal to reduce to zero the number of uncontrolled approximations such that the theoretical predictions can be reliably tested in experiments.
This paper makes a step in this direction for cosmological transitions in sectors described by strongly-coupled Yang-Mills or QCD-like theories. The latter appear in many dark matter models (see, e.g., [10][11][12]). We consider scenarios where the dark matter is constituted e.g. by dark glueballs, pions or baryons.
Whenever the theory is confining, one expects a confinement/deconfinement transition as the Universe cools below the theory's dynamical scale. If the transition is first order, it may generate GWs, which are the study objects of this paper.
When the gauge theory includes (approximately) massless quarks, the strongly-coupled dynamics is such that the (approximate) chiral symmetry is broken at a scale that might or might not coincide with the gauge theory's dynamical scale. We consider both the case in which the confinement phase transition implies the chiral symmetry phase transition and the case in which it does not. The first case also includes the Peccei-Quinn transition in the simplest composite axion model with hidden gauge group [13][14][15]. The second case includes the Peccei-Quinn first-order phase transition of the recently proposed holographic axion model [16,17], where the axion appears as a pseudo-Nambu-Goldstone boson associated with the chiral symmetry breaking of an extra pair of quark/antiquark fields.
To be more specific, we consider theories where the rank of the gauge group is sufficiently large such that the planar approximation is reliable. In this case, a class of very interesting models are the ones admitting a holographic description. 1 As a prototype, we consider the top-down theory that, in the deep IR, better resembles planar Yang-Mills, or planar QCD if we consider the flavored version, known as the Witten-Sakai-Sugimoto (WSS) model [18,19]. In the case of a YM or QCD-like dark sector outside the regime where the holographic description is completely reliable, the latter can be employed as an effective approach to the strong dynamics of the theory.
The WSS model has been widely used to study various aspects of QCD at low energy, with notable success. In the present context, the model is interesting because it features two first-order phase transitions, the first one associated with confinement/deconfinement and the second one associated with chiral symmetry breaking/restoration.
In most of the cases that we investigate, we employ the WSS model not as a proxy for QCD but as a model for a dark sector. Being a so-called top-down model, the WSS has the advantage that computations performed in the planar limit at strong coupling are reliable, in the sense that there is a precise control on the validity regime of the various approximations, something which usually does not occur in effective phenomenological models or bottom-up holographic theories. In fact, this property eliminates one of the sources of uncertainty in the calculation of the parameters for the GWs spectra when dealing with strongly-coupled theories, and it constitutes the main motivation for this paper.
In a previous work [20], we have addressed the problem of the nucleation of bubbles of true vacuum associated with both the confinement/deconfinement phase transition and the chiral symmetry breaking/restoration phase transition in the WSS model. 2 In the present work, we use those results to compute the stochastic GW spectra, due to bubble collisions and sound waves, in several beyond Standard Model scenarios featuring the WSS model. As we will see, the main conclusion of our analysis is that there is a large window of the WSS parameter space where the GW signals may be accessible in near-future experiments. Moreover, the model allows for the generation of GWs compatible with the possible observation recently reported by NANOGrav [21].
The paper is organized as follows. In section 2, we introduce the WSS model and summarize the steps of the analysis needed to find the GW spectra. 1 It is a widespread belief that every gauge theory in the planar limit admits a perturbative string theory description. The latter can or cannot have a low-energy limit where classical gravity is reliable, depending on the theory's details. For example, pure Yang-Mills theory does not have a classical gravitational description. Nevertheless, there are infinite classes of theories, which we call holographic, admitting such a description. 2 In [20], the Randall-Sundrum scenario has been briefly discussed as well. In that context, we were able to compute the derivative term in the effective bounce action in the high temperature regime. That term was missing in previous literature on the subject.
In section 3, we consider three different dark matter scenarios. These are cases where the chiral symmetry transition, if present, is implied by the confinement one. In subsection 3.4, we discuss the results for the GW spectra. Figure 1 encodes in a global view some benchmark results of the investigation.
In section 4, we consider two scenarios where GWs come from the chiral symmetry breaking/restoration phase transition. In one of them, the chiral transition is followed by a separated confinement/deconfinement one. We thus investigate the fascinating possibility of detecting a GW spectrum with two peaks. In this case, moreover, we outline the fact that the usual assumption of adiabatic expansion of the Universe from the first phase transition to present times cannot be used anymore: the presence of a second phase transition requires a refinement of the usual redshift factors in the formulae for the GW spectra. The results for the GW spectra are reported in subsection 4.3 and in figure 3.
We will conclude with a summary and some observations in section 5. We collect all the holography-related details of the WSS model in appendix A. Appendix B provides an overview of all of the relevant formulae used to obtain the GW spectra. In particular, in B.1.4, we discuss how the occurrence of two separated phase transitions affects the quantities that determine the GW spectra, providing explicit formulae for the modified redshift factors advocated in section 4. Finally, in appendix C, we review the results of [20] that are useful for the present paper and provide approximate analytical expressions of the relevant GW parameters for the confinement/deconfinement transition in the small temperature regime.

The WSS model and its embedding in cosmology
In this section, we describe the features of the Witten-Sakai-Sugimoto model that are needed in order to understand the calculation of the GW spectra. More details on the model and on the bubble configurations nucleated in the phase transitions are reported in appendices A and C.
The WSS model is a (3+1)-dimensional non-supersymmetric Yang-Mills theory with gauge group SU (N ) coupled to a tower of adjoint Kaluza-Klein (KK) fields and to N f fundamental flavors (quarks) [18,19] (see also [22] for a review). The model possesses five independent parameters. Two of them are actually dimensional quantities: M KK , which represents the dynamically generated scale providing the mass of the first glueball and that of the first KK field, and L which gives the scale of chiral symmetry breaking f χ,L , as we will discuss in a moment. The other three, dimensionless parameters are given by N , N f , and the 't Hooft coupling λ at the scale M KK . We will consider the regime The properties of the model at low energies are very similar to the real-world QCD ones since they include confinement, mass gap, and chiral symmetry breaking. We can actually write more precisely the last condition in (2.1) as (see e.g. [23]) which holds in the confined regime. 3 One of the main motivations for studying the model in this paper is that it exhibits two first-order phase transitions. The first one separates the low temperature confined phase of the theory from the high temperature deconfined one. The critical temperature for the transition is [24] 3) The second first-order phase transition separates the chirally symmetric phase from the phase where chiral symmetry is broken [24]. In the general case, L is a free parameter of the model that can be used to separate the confinement scale from the chiral symmetry breaking one. When L > 0.97M −1 KK the confinement/deconfinement transition implies the chiral symmetry breaking/restoration one. In contrast, when L < 0.97M −1 KK , the two transitions are independent, with the chiral symmetry breaking/restoration one occurring at the temperature The parameter L has the maximal value L = πM −1 KK , when the scale of chiral symmetry breaking reads (2.5) In the opposite limit L πM KK , we have [16,25] (2.6) So far, we have been assuming that all the N f quarks condense at the same scale, dictated by the same value of L. But of course we can actually have several classes of quarks with different values of L.
To summarize, the phase diagram of the model is the following: , the theory is confining and chiral symmetry is broken; 2π , the theory is deconfined and: if T < 0.1538 L , chiral symmetry is broken; if T > 0.1538 L , chiral symmetry is preserved.

Cosmological WSS phase transitions
In this subsection, we describe the general framework needed to calculate the GW spectra, also fixing our notation. We leave most of the technical details, which are quite standard, to appendix B, for the benefit of the reader who is not familiar with this type of computations. We will consider a cosmological setting where the Universe starts at some high temperature, in which the WSS is in the deconfined phase, and then cools down. Depending on the scenario that we consider, the WSS sector will undergo one or two first-order phase transitions. They are triggered by the nucleation of bubbles of true vacuum (confined phase or chirally broken phase, depending on the transition) in the plasma, which is in the metastable false vacuum (deconfined or chirally symmetric). These bubbles will expand and eventually fill all the Universe, leaving it in the true vacuum state. The percolation temperature T p is defined as the temperature of the Universe when this process completes. We will compute it case by case, using the formulae discussed in appendix B.1.2.
The cosmological evolution of the Universe is described, as usual, by the Friedmann-Lemaitre-Robertson-Walker (FLRW) metric 5 where R(t) is the cosmic scale which defines the Hubble scale H(t) =Ṙ(t)/R(t). The latter is determined by the total energy density through the Friedmann equation with M P l ≈ 2.4 · 10 18 GeV. The energy density ρ takes contributions from the standard model and from the dark sector. In the sector described by the WSS model, the energy density in the deconfined and in the confined phase at order O(N 2 ) reads, respectively, In the limit (2.2), the contribution of N f quarks to the energy density in the high-temperature regime and in the low-temperature one at order O(N f N ), in the case L = πM −1 KK , read (see e.g. [23,26]) As discussed above, when L πM −1 KK , there is an intermediate phase where the gauge theory is deconfined and the quarks are condensed. In this case, the energy density is not known analytically. However, it can be computed numerically starting from the energy density of the chirally-unbroken configurations. In particular, it reads, where ∆S is defined in (C.23) and In fact, T P ∆S gives exactly the difference of free energies of the flavors in the broken and unbroken phases. Using the fact that the energy is the derivative of the free energy w.r.t. the temperature, the second term on the r.h.s. of the first relation (2.11) is the difference of the energies in the two phases, so that adding the known contribution of the unbroken phase, one is left with that of the broken phase. As we will comment on in section 4.2, the energy density of condensed quarks with L = πM −1 KK in the confined phase will always be subleading and can be neglected.
From (2.9) and (2.10), we see that the confined phase of the WSS model carries a temperature-independent contribution to the energy, which would act as a cosmological constant after the phase transition. Since the measured cosmological constant almost vanishes, the zero-point energy has to be shifted accordingly. As a result, the energy density in the deconfined and chirally symmetric phase reads 6 ρ deconf unbroken = ρ rad,glue + ρ rad,SM + ρ rad,χ + ρ 0,glue + ρ 0,χ , (2.12) where ρ rad,SM = π 2 30 g SM * (T ) is the Standard Model contribution, given by the temperature-dependent number of relativistic degrees of freedom g SM * . The factor is defined as the ratio between the temperature T of the dark sector and that of the Standard Model T V . As we will see, ξ can (and in some cases must) be different from 1. The energy density in the deconfined and chirally broken phase reads ρ deconf broken = ρ rad,glue + ρ rad,SM + ρ b,χ + ρ 0,glue , (2.15) 6 In the most general case, we have quarks of both L = πM −1 KK and L πM −1 KK kind. Hence, the contribution ρ 0,χ is not simply given by (2.10b), because the latter holds only for the L = πM −1 KK . The L πM −1 KK contribution is suppressed by a M KK /f χ,L factor and can be usually neglected.
whereas in the confined and chirally broken phase it is where g * (T ) receives contributions from the Standard Model degrees of freedom and eventually from relativistic particles from the dark sector. We will investigate several scenarios where (2.12), (2.15) and (2.16) will be used. The cases will differ for the values of the parameters N f , N , λ, and the number of degrees of freedom involved.
Away from the phase transitions, the universe evolves adiabatically, i.e. according to the conservation of the entropy where, in general, g S * (T ) = g * (T ), see appendix B.1.4. During the phase transition, an amount of energy is released and the plasma gets heated up. The temperature T R of the plasma at the end of the transition is called reheating temperature and is found via the conservation of energy. This point will play an important role in section 4, where we will consider the case in which the universe undergoes two first-order phase transitions. As we will see, the presence of the second phase transition modifies the redshift of the GW signal compared to the adiabatic evolution one, usually assumed to be valid after the single phase transition.
As we detail in appendix B, the efficiency of the phase transition depends on the ratio Γ/H 4 , where Γ is the bubble nucleation rate. In the case in which a single field describes the transition, the bubble nucleation rate Γ can be computed in the semiclassical approximation using the formalism developed in [2][3][4][5][6]. The confining phase transition of the WSS model involves several fields. In [20], as reviewed in appendix A, we took an effective approach inspired by [27] where only a single field is involved. The formula for the bubble nucleation rate is reported in (B.1), which involves a comparison between the efficiency of quantum and thermal fluctuations. The former are given by the O(4)-symmetric solution, and the latter by the O(3)-symmetric one. In the analysis, we always have to verify which kind of bubble dominates.
Depending on the phase transition's efficiency, the universe may remain trapped in the false vacuum for a long time after it reaches the critical temperature, featuring supercooling. In this case, the energy density may include a temperature-independent contribution, which may start to dominate, acting as an effective cosmological constant that makes the universe inflate. 7 As a result, it is not guaranteed that the phase transition completes, hence in the analysis, we will always have to check that it actually does. Technically, this is done through formula (B.21) discussed in appendix B. Depending on whether percolation enters the vacuum-dominated phase or not, the percolation temperature is computed, respectively, by (B.13) or (B.16). In performing these and the following calculations, we use the Chapman-Jouguet formula (B.25) for the velocity of the bubble. 8 Gravitational waves are produced during the propagation of nucleated bubbles in the plasma in three ways: collisions among bubbles, collisions of plasma sound waves, and turbulence in the plasma. Unfortunately, the turbulence contribution to the gravitational waves spectra is currently not well-understood. Typically it is deemed as subdominant. We will only consider the contributions coming from bubble collision and from the sound waves for these reasons. The formulae for the spectra in these two cases are given, respectively, by (B.38) and (B.39).
As we discuss in appendix B, it is not easy to estimate how the energy is distributed among the various contributions. Comprehension of the bubble dynamics and, most importantly, interaction with the plasma is one of the major open problems in the field so that the results are affected by huge incertitudes. For this reason, in this paper, the results for the spectra are presented separately for the bubble collision and sound waves contributions, pretending that all of the energy is concentrated in one of them in turn. The true spectra will obviously be in between these two "extremal" cases.
The GW spectrum depends crucially on a parameter, usually called α, which accounts for the amount of energy released in the transition. We are going to use its expression in terms of the trace of the energy-momentum tensor (formula (B.24)), adjusting in any place the number of relativistic d.o.f. at the relevant temperature scale. 9 As we will see, the spectrum with a larger magnitude is that associated with sound waves.
In the next two sections, we present the analysis in the various scenarios. From the WSS model perspective, the main difference among them is given by the choice of the parameters M KK , f χ,L , N f , N , and λ. Actually, for what concerns the next section, the latter two enter through the combination g ≡ λN 2 . (2.18) As a general framework, although both N and λ are required to be large parameters, it is natural not to introduce a huge hierarchy of scales. Thus, we tend to prefer (but not limit ourselves) to consider not-too-large values of the parameter g, starting from g 100.

GWs from deconfinement/confinement phase transition
In this section, we present the GW spectra produced in three possible dark scenarios, which we name Dark HQCD 1, Dark glueballs and Dark axion. The "H" in HQCD stands for "Holographic", to underline the fact that there are extra modes w.r.t. standard QCD-like theories. In these scenarios, gravitational waves are always associated with the confinement/deconfinement phase transition. It is important to outline that the WSS model realizes explicitly, in a specific regime of parameters, scenarios which have been previously proposed in the literature (see e.g. [10][11][12] for reviews). While it would be very interesting to further study the phenomenological implications of this regime of parameters, in this paper we just concentrate on the gravitational waves spectra. Thus, in the following subsections we are going to sketch the different scenarios, discussing the main information needed for the computation of the GW spectra. The latter are determined with the formulae collected in appendix B and the results are presented in subsection 3.4.

Dark HQCD 1
QCD-like theories with N f flavors can provide different dark matter candidates. Depending on the details of the models, the main fraction of dark matter can come from dark baryons, nuclei, mesons, and so on. Analogously, the dynamically generated scale, which in the WSS model is denoted as M KK , varies considerably among the various theories, typically from about 100 MeV to about 100 TeV. In this subsection we consider the WSS model with N f flavors, in the regime (2.1), as providing a strongly-correlated large N dark QCD-like sector. Previous studies of gravitational wave spectra in similar scenarios include [31][32][33][34][35][36].
We have analyzed the spectra of GW produced in the phase transition for the dynamical scale values M KK = 10 n GeV , n = −1, 0, ..., 6 , (3.1) and for g = 10 m , m = 2, 3, 6, 10 . (3. 2) The case g = 10 2 is the only one where the Universe at the time of bubble percolation is in a radiation domination phase, hence we employ formula (B.16) to determine the percolation temperature; in all the other cases, the Universe is in a vacuum domination era and we have to employ formula (B.13). For g = 10 2,3,6 , the relevant bounce solution is the O(3)-symmetric one, while for g = 10 10 the O(4)-symmetric configuration dominates.
In determining the reheating temperature according to formula (B.23), care must be taken to count the correct number of degrees of freedom both in the Standard Model and in the dark sector. In fact, in the confined phase of the dark sector there can be glueballs, KKmodes and mesons which become relativistic at the reheating temperature. This happens for g = 10 6 and g = 10 10 . In the first case, only the lightest glueball and KK mode must be included, together with the lightest mesons. In contrast, in the second case, the reheating temperature is about seven times M KK . At this scale, many glueballs from Table 2 in [37] as well as many mesons must be included, giving hundreds of d.o.f. Unfortunately, the spectrum of KK modes is not known in detail. The first KK modes have mass of one M KK , but we have no definite information on the number of degrees of freedom at 7M KK . We give a very rough estimate of this number assuming that the density of KK modes has the same dependence on the energy as the spectrum of glueballs. We then double the number of degrees of freedom to account for the fermionic glueballs and KK modes. The same is done for the mesons. However, we underline that the incertitude associated to the number of degrees of freedom introduces an error that does not spoil the order of magnitude of our results.

Dark Glueballs
Another well-motivated class of dark matter candidates is represented by stable glueballs, the bound states of SU (N ) Yang-Mills theory. The WSS model with N f = 0 is therefore suitable for describing such a scenario and for performing in this context reliable calculations. Being derived in the quenched approximation, the results of section 3.1 can be seen as also concerning a scenario where the non-interacting dark sector is constituted by a SU (N ) Yang-Mills theory without flavors.
The latter can also model the case where the dark matter is actually self-interacting, a possibility which helps softening the problems of the ΛCDM model with small-scale structures [38]. In this scenario, phenomenology can be satisfied for glueball masses ranging from keV to fraction of GeV. When the order of the latter is around one MeV or smaller, one has to take care of phenomenological constraints related to the effective number of neutrino species and coming from Cosmic Microwave Background (CMB) measurements, and from measurements of the relative abundance of elements in Big Bang Nucleosynthesis. They imply that the dark sector cannot be in thermal equilibrium with the visible sector. In particular, the dark sector temperature T has to be smaller than T V , the visible sector one [39,40]. As a result, non-gravitational couplings among the two sectors have to be absent or extremely small. Whenever this is the case, gravitational waves produced in first-order transitions can be one of the few means at our disposal in order to observe direct signals coming from the dark sector. Previous studies of the GW spectra in similar cases within the context of simple effective models can be found in [39,40].
In this section we will investigate cases with dynamical scale values M KK = 10 n keV , n = 0, 1, 2, 3, 4, 5 . The other main difference with respect to the analysis performed in section 3.1 is given by the fact that the ratio ξ = T /T V can be smaller than one. We assume that ξ stays constant during the bubble nucleation and GW observation process. We explicitly explore benchmark cases where ξ = 10 −1 , g = 10 3,10 .
We have also checked that the smaller the value of ξ is, the more the signal is suppressed. For ξ = 10 −5 , for example, the signal will be completely invisible in near future facilities. As for the Dark HQCD 1 scenario, the cases with not-so-large values of g are likely to be more phenomenologically fit (e.g. they do not suppress too much the dark matter self-interactions). Let us briefly describe the main features of the calculation. The decoupling of the dark and visible sectors implies that whenever we consider plasma effects, the plasma in question is just the one of the dark sector. As a consequence, there are two relevant α parameters (formula (B.26)), denoted as α and α D , measuring respectively the energy released in the transition w.r.t. the visible sector energy density only and w.r.t. the dark sector energy density only. The velocity of the bubble wall is determined by formula (B.25) with α replaced by α D . The same is true for the efficiency parameter κ v (formula (B.40)) for the sound wave spectra.
For g = 10 3 , in all the cases the Universe is found to be in a radiation domination era at the time of percolation. In fact, values of ξ < 1 enhance the contribution of the SM energy density of radiation against the dark sector vacuum energy density. The bubbles in these cases have O(3) symmetry. Only for the case of g = 10 10 , ξ = 10 −1 the Universe is in a vacuum domination era, the percolation temperature is very small due to supercooling and O(4) bubbles dominate. Also, in the cases of g = 10 3 , ξ = 10 −5 and g = 10 10 , ξ = 10 −1 the reheating temperature is considerably different from T p , so that we have to consider many glueball and KK modes from the dark sector. However, due to the damping factor ξ, the contribution of the dark degrees of freedom is quite suppressed w.r.t. the contribution of the SM particles.

Dark Axion
In this section, we analyze a third range of dark sector dynamical scales, relevant for composite QCD axion models. The benchmark model is the one discussed in [14] building on the model in [13] (see also [15], and [41,42] for recent reviews).
In its simplest realization, the model comprises a dark SU (N ) Yang-Mills sector and four massless flavors in its fundamental representation. Three of them form a triplet of the QCD SU (3) c gauge group, whereas the fourth constitutes a singlet. The global symmetry includes an axial U (1) A , which plays the role of the Peccei-Quinn symmetry. In fact, the latter is anomalous and spontaneously broken by the flavor condensation due to the strong dynamics of the dark SU (N ). The associated pseudo-Nambu-Goldstone boson is then a composite axion. In this scenario, the confinement/deconfinement transition of the dark SU (N ) theory implies the Peccei-Quinn phase transition, which is of the first order. Previous studies of GW spectra from Peccei-Quinn transitions in effective theories (possibly of bottom-up Randall-Sundrum type) can be found in [29,43,44].
In the model of [13], the axion decay constant f a is related to f χ by Thus, from (2.5), we read Consistency with phenomenology requires f a 10 8 GeV. Moreover, formula (2. 2) with and therefore we are led to consider dynamical scales M KK 10 9 GeV. We will consider two benchmark values of g, g = 10 3 , g = 10 8 . (3.8) The details of the calculations are very similar to the ones in section 3.1. In all the cases, the Universe is in an energy domination era at the time of percolation. For g = 10 3 (g = 10 8 ) the O(3) (O(4)) bounce dominates. In order to determine the reheating temperature for g = 10 8 , we have to take into account glueball, KK and mesonic degrees of freedom.

Results for the spectra
In this section, we describe the results for the GW spectra generated by the first-order confinement/deconfinement transition of the holographic model. As we have already mentioned, we do not consider the contribution from turbulence in the plasma and we separately consider the contributions from bubble collisions and sound waves.
For what concerns the sound waves contribution, there is a further incertitude due to the unknown source duration. Until very recently, the source was expected to last for a long time in Hubble units. Under this assumption, most of the literature has employed the formulae reviewed in [7,45]. However, it has been recently pointed out that the source can be quite short, see e.g. [9,[46][47][48]. Accordingly, the power spectrum is quenched by the short time factor (B.41). In this paper, an agnostic attitude is taken and both spectra, with and without quenching factor, are presented. This allows us to have an idea of the possible range of the signal and to compare the results with previous literature.
In summary, three types of spectra are calculated: the one from bubble collisions Ω c , the one from sound waves without quenching factor Ω sw and the one from sound waves with quenching factor Ω sw,q . As a general trend, Ω c is found to give the smallest peak signal.  Figure 1: Examples of GW power spectra h 2 Ω GW due to bubble collisions (Ω c , dashed lines) and sound waves in the case of short source duration (Ω sw,q , continuous lines) and long source duration (Ω sw , dotted lines). Expected sensitivities (PLISCs) for a number of experimental facilities are reported for comparison [49]. From left to right, the spectra correspond to the following parameters: (M KK /GeV, g) = (10 −6 , 10 10 ) (blue lines), (10 2 , 10 6 ), (10 6 , 10 6 ) (green lines), (10 9 , 10 3 ) (red lines).
Moreover, the peak frequency increases with M KK and the amplitude of the signal increases with g.
In figure 1 we report examples of power spectra. In the plot, a few benchmark values of the parameters M KK , g are chosen to show the detectability potential of the GW emissions. A number of experimental sensitivities are shown for comparison.
The first clear result is that in various cases the GW signals are going to be detectable in near future experiments, with the possible exception of the composite axion model.
Notice that Ω c and Ω sw,q approximately span an order of magnitude in power of the signal around the peak, represented in the figure by the regions in between the dashed and continuous curves. Notice, moreover, that the upper value of the signal for the Ω sw spectra (dotted lines) is greatly amplified w.r.t. the quenched case Ω sw,q (continuous lines); the true signal from sound waves is expected to be in between the two types of lines. The total signal is expected to be a combination of the one from sound waves and the one from collisions.
The blue lines at the left of the plot show a representative case for a small dynamical scale value, M KK = 1 keV, relevant for the Dark Glueballs scenario, for g = 10 10 and for the value ξ = 0.1 of the ratio between the dark and the visible sector temperatures. It is clear that the signal is potentially detectable by pulsar timing array experiments such as IPTA and SKA. Actually, the most "optimistic" scenario where almost all the energy of the process goes into GWs from sound waves is of great experimental interest. In fact, in this case the signal could be visible in current single experiments such as NANOGRAV, EPTA and PTTA. Actually, very recently, the results of 12.5 years observations by NANOGRAV have been reported in [21], showing strong evidence for a stochastic spectrum compatible with GW signals with frequency peak around 10 −9 − 10 −8 Hz and average energy density h 2 Ω GW ∼ 10 −10 . If, among the possible sources of this signal, there is space for a cosmological strongly first-order phase transition in a dark sector -as it has been recently suggested in [50][51][52] -our Dark Glueball model could be viewed as a possible candidate.
Although it is not shown in the figure, the same possibility of detection happens if the scale is raised up to around M KK ∼ 100 keV. Moreover the signal falls in the same experiment sensitivity curves also for the minimal value of g = 10 3 (again for ξ = 0.1), although with a smaller magnitude and for a shifted value of scale M KK ∼ 10 2 − 10 3 keV (otherwise the peak frequency is too small).
The two sets of green lines at the center of the plot correspond to the parameter value g = 10 6 and energies respectively of M KK = 10 2 and M KK = 10 6 GeV, relevant for the Dark HQCD 1 scenario. The first case is going to be detectable already by LISA and clearly by the more sensitive experiments such as BBO and DECIGO. The same remains true down to M KK ∼ 10 GeV and g ∼ 10 2 (not shown in the plot). The second case of M KK = 10 6 GeV is detectable by ET or CE facilities. Of course, all the intermediate energies can be detected, and this remains true even for smaller values of g down to 10 2 and larger values of M KK 10 7 GeV. For g = 10 10 the signal is visible at LISA starting from M KK ∼ 1 GeV. Thus, a few near future experiments (LISA and ET for example) are going to be able to fully probe strongly coupled dark QCD-like sectors (with large ranks) in the energy range M KK ∼ 1 − 10 7 GeV.
Finally, the three red lines at the right of the plot correspond to g = 10 3 and M KK = 10 9 GeV, and are relevant for the Dark Axion scenario with f a ∼ 10 8 GeV. Only in the optimistic case in which the duration of the sound waves' source is long, the spectrum falls within the sensitivity curve of CE. Since we expect the real signal to be in the region between the three curves, this case is unlikely to be detectable in near-future experiments. Moreover, if M KK increases, such that f a > 10 8 GeV, the curves are shifted to larger values of the peak frequencies. As a result, the Dark Axion scenario is not favorable for producing detectable gravitational waves.
The picture that emerges is that the upcoming generation of experimental facilities is going to probe a large parameter space of dark sectors admitting a holographic description. Moreover, in various (optimistic) scenarios, stochastic GW background generated in this type of models can be detectable by the advanced version of currently running experiments.

GWs from chiral phase transition
In this section, we consider scenarios that display a chiral symmetry breaking/restoration phase transition separated from the eventual confinement/deconfinement one. This implies the fascinating consequence of having two distinct peaks in the spectrum of stochastic GWs. Firstly, we discuss the possible scenarios and then we present the results for the spectra.

Dark HQCD 2
The scenario that we consider in this subsection is a close cousin of the Dark HQCD 1 scenario of subsection 3.1: the WSS model describes a dark sector, very weakly interacting with the Standard Model (in the most extreme case, even interacting with the Standard Model only gravitationally). The difference with respect to what has been discussed in section 3.1 concerns the choice of the WSS parameter L. In section 3.1, the latter was taken to be L = πM −1 KK , corresponding to the chiral symmetry breaking scale f χ given in (2.5). In contrast, here we will consider cases with L πM −1 KK , for which the chiral symmetry breaking scale f χ,L is given in (2.6). As said, this implies that the chiral phase transition is separated from the confinement one.
An important difference with respect to the scenario of section 3.1 is that the evolution of the Universe cannot be considered to be adiabatic from the time of the chiral symmetry breaking transition to the present time, since there is a second first-order phase transition. This calls for a correction of the standard formulae for the redshift of the signal, which are derived under the assumption of adiabatic evolution. In fact, the adiabatic assumption holds from the time of the chiral symmetry breaking transition to the percolation time of the confinement transition. Then, assuming fast reheating in the confinement transition, the temperature has a sudden jump from the percolation temperature to the reheating temperature. Finally, from this time to the present day, the Universe continues to evolve adiabatically. In appendix B.1.4 this behavior is reflected in formulae (B.34), (B.35) for the frequency and power spectrum redshifts.
A consequence of these formulae is that the magnitude of the chiral symmetry breaking transition signal decreases if the value of g = λN 2 increases. This is due to powers of the ratio of the percolation and reheating temperatures of the confining transition, T p,conf , T R,conf , appearing in formulae (B.34), (B.35) (in a coefficient which we have called δ in (B.36)). As we semi-analytically estimate in appendix C.3, an increase of g implies more supercooling, hence T p,conf decreases and at the same time T R,conf increases, resulting in a suppression of the GW signal. For this reason, in the present scenario we are going to describe the case where λ and N are such that g has a "small " value. In particular, we will investigate the representative case λ = N = 10 , g = 10 3 . It is convenient to introduce dimensionless quantities, such that the critical temperature for the chiral symmetry breaking transition corresponds toT = 1. The condition that the chiral symmetry breaking transition happens above the deconfinement transition gives the constraint that with the choice (4.1) corresponds tof χ > 0. 13. In fact, the signal is enhanced if the chiral symmetry breaking scalef χ is large. The validity of the quenched approximation we are assuming for the flavors constrains the magnitude of this parameter. In particular, the requirement that the approximation works at the percolation temperature and at the reheating temperature sets the limitf χ ≤ 60 for the choice of parameters (4.1). This comes from the requirement that the energy density of the flavors is subleading with respect to the one of the gluonic degrees of freedom, see section 2.1. Thus, we will consider the benchmark values A noticeable difference with respect to the cases analyzed in section 3 is that the energy released in the transition is much smaller than the energy of radiation, since the former comes from the flavors, which are quenched, while the latter mostly comes from the gluons. As a result, the parameter α is much smaller than one in this case and the bubble velocity sometimes is not very close to unity. Since the energy released in the transition is small as compared to the total energy, we expect the reheating temperature to be close to the percolation temperature.
Regarding the counting of degrees of freedom, in the case at hand, by normalizing the entropy density as at the time of emission we have the three contributions from the Standard Model, gluons and flavors g S * = g * ,SM + g S * ,glue + g S * ,χ , (4.6) with [24] g S * ,glue = 5 · 2 6 π 2 From the energy density of section 2.1 we read g * = g * ,SM + g * ,glue + g * ,χ , (4.8) with (4.9b)

Holographic Axion
Another scenario where a chiral symmetry breaking/restoration takes place is the holographic QCD axion model of [16], which we call HoloAxion in the following. The WSS theory is considered as a model for the strong interactions of the Standard Model, including the QCD axion physics. The axion arises as a composite particle, analogous to the η , coming from an extra flavor with L πM −1 KK so that it condenses at a large scale f a = f χ,L Λ QCD . In contrast, the SM quarks are embedded in such a way that the related chiral symmetry breaking scale is given by f χ in (2.5). The condensation of the axion is a Peccei-Quinn first-order transition which can therefore generate gravitational waves.
The energy density of the false vacuum configuration in this case reads formally as (2.12). Let us briefly comment on each contribution. Since the QCD sector of the theory, gluons and quarks, is described by the WSS model, the related relativistic degrees of freedom are not counted in g SM * (which then has 27.75 as its maximal value) in ρ rad,SM . Concerning ρ rad,χ , the number of flavors in (2.10a) is N f = 7, because we have six QCD quarks plus an extra flavor that provides the axion. The contribution ρ 0,χ is given by (2.10b) with N f = 6, because the latter holds only for the case L = πM −1 KK . The remaining flavor gives a contribution analogous to (2.10b) but suppressed by a factor of M KK /f χ,L , hence it can be neglected.
Since in this scenario the WSS model describes the strong sector, the usual, uncontrolled extrapolation of the regime of validity of these formulae to the real world parameter values is performed. This amounts to quitting the planar regime by setting N = 3. Then, the parameters λ and M KK are determined by fitting the ρ-meson mass and the value of f π = f χ , giving [19] λ = 33.26 , M KK = 0.949 GeV .  coming from axion phenomenological constraints.

Results for the spectra
Let us comment on the behavior of the spectra that we find in the scenarios where a chiral symmetry breaking/restoration transition occurs. In the Dark HQCD 2 scenario, two separated phase transitions occur and the signal is given by the sum of the signals of the two phase transitions. Since we work in the quenched approximation (2.2), the chiral symmetry phase transition is characterized by smaller released energies and therefore smaller signal magnitudes with respect to the confinement/deconfinement one. The peak of the signal of the chiral symmetry transition is at higher frequencies than that due to the confinement/deconfinement transition. Being smaller, the former might be negligible with respect to the tail of the confinement signal and therefore the chiral symmetry phase transition would be effectively unobservable. 10 Since the signals associated to bubble collisions are suppressed with respect to the ones due to sound modes, we discuss only the latter.
Examples of the signals for different values of the parameters, with and without the correcting factor (B.41) for the duration of the transition, are reported in figure 2. Clearly, larger values off χ are more effective in separating the peak due to the chiral symmetry transition from that due to confinement. Figure 3 offers an example of the scenarios that we have been discussing in this section, presenting the comparison of the computed spectra with the sensitivity curves of experiments. The green curves correspond to a representative two-peak case in the Dark HQCD 2 scenario, namely that where M KK = 100 GeV andf χ = 60. It displays a large peak due to the confinement/deconfinement transition at frequency f ∼ 10 −3 Hz which fits into the sensitivity curve of LISA, and a smaller peak due to the chiral symmetry transition at frequency f ∼ 10 −1 Hz which does not fit into the LISA sensitivity curve but is expected to be visible by the next generation facilities such as BBO and DECIGO. The conclusion is that the two-peak signal is certainly within reach of the next generation facilities at least for a certain region of parameter space.
Concerning the HoloAxion case, the result for the extremal case where the axion decay constant takes the lower allowed value f a ∼ 10 8 GeV is displayed in red in figure 3. The frequencies of the peak of these curves are too large and their magnitudes are too small to be captured by near-future facilities like ET or CE, even in the optimistic case in which we do not include the quenching factor (B.41) due the duration of the sound waves. Moreover, as f a increases, the peak frequencies increase as well, hence going further away from the sensitivity curves of the experiments. Thus, we conclude that the Peccei-Quinn transition in the HoloAxion scenario cannot be seen in near future experiments.

Conclusions
Cosmological first-order phase transitions generate stochastic gravitational wave backgrounds potentially visible in present and next generation experimental facilities. Dark sectors, as hidden sectors interacting with the standard model very weakly, are good candidates where to explore such transitions. Many dark sector models in the literature are Yang-Mills or QCD-like theories. If the rank of the gauge groups of these theories is sufficiently large, the planar limit constitutes a good approximation to their dynamics. In this paper, we have considered the scenario where a dark sector admits a top-down holographic dual description in the gravity regime. This means that the theory is in the planar limit and there is a gap in the spectrum of hadron masses. When the theory admits such a dual description, we have full control on its strongly-coupled dynamics, without the need to employ effective models and uncontrolled approximations. But even if the theory is not exactly in this regime, one can view the holographic description as an effective tool to model the strong coupling dynamics -this latter approach has been used extensively for QCD.
Describing dark sectors by means of dual gravitational theories opens up the possibility of studying their dynamics at strong coupling. In this paper, we concentrated on the production of gravitational waves in first-order transitions. Using the well-known Witten-Sakai-Sugimoto holographic model, we have investigated two types of transitions. The first type is the confinement transition, possibly implying a chiral symmetry breaking transition. The second type is a chiral symmetry breaking transition separated from the confinement one -the latter happening at a later time in the cosmological evolution.
Making use of the bubble configurations studied in the companion paper [20], we have been able to calculate all the relevant parameters necessary for the determination of the GW spectra. The latter are usually affected by a number of assumptions and sometimes uncontrolled approximations. The holographic approach allowed us to erase from this number the use of uncontrolled approximations to the strong dynamics of the dark theory.
The results of our investigation are partially in line with other studies in the literature. In table 1 we report the benchmark cases displayed in figures 1 and 3. In the case of the single confinement transition, there is a large part of the parameter space of the theory where the GW signal is going to be detectable in the next generation facilities (see figure  1 for examples). These include pulsar timing arrays as well as space-and ground-based interferometers, depending on the dynamical scale of the theory. Interestingly, a window of parameter space can produce a signal within the current NANOGrav sensitivity, explaining the recent potential observation in this experiment.

Summary of the benchmark cases Scenario
Dynamical scale Chiral scale Experiment Dark HQCD 1 10 2 ,  When the chiral symmetry breaking transition is separated from the confinement one, the model predicts two distinct peaks in the GW spectra. Detection of both peaks would represent an exciting smoking gun for the models with two transitions. The gravity regime allows to explore faithfully a branch of parameter space where the chiral symmetry signal is smaller than the confinement one. Nevertheless, we have shown that there are certain values of parameters allowing for observation of the two peaks, for example by space-based interferometers ( figure 3). It would be interesting to study the correlations of the two peaks, which could distinguish the holographic model from other models with two phase transitions.
Finally, we have considered Peccei-Quinn transitions in two distinct axion models: a standard composite axion from a hidden sector [13][14][15] and the recently introduced holographic axion model [16]. Unfortunately, in both cases, the lower bound on the axion decay constant around 10 8 GeV corresponds to a peak frequency which is too large for detection in the near future. In this respect, the model is distinct from the holographic bottom-up (phenomenological) ones recently investigated in [29,43], where the possibility of tuning a very small parameter, measuring the departure from conformality, allows to produce signals within the sensitivity of ET or CE.
In this paper we have started to use top-down holographic models to study dark (hidden) sectors. It will be clearly interesting to employ this approach to first characterize the model parameter space compatible with current observational constraints, and then produce predictions for observables in the strong coupling regime of the theory.

A More on the Witten-Sakai-Sugimoto model
In this appendix, we review the holographic construction of the Witten-Sakai-Sugimoto model employed in the present paper.
Let us briefly describe the string theory embedding of the WSS model. A stack of N D4-branes wrapped on a circle S 1 x 4 with coordinate x 4 ∼ x 4 + 2π/M KK give rise to the fields that transform in the adjoint representation of the gauge group. Fields transforming in the fundamental representation of the gauge group are introduced through pairs of D8/anti-D8branes. These are transverse to the circle S 1 x 4 in such a way that the D8-branes and anti-D8-branes are separated by a distance L ≤ πM −1 KK along S 1 x 4 . When the N f fundamental fields are massless, the theory exhibits a U (N f ) L × U (N f ) R classical global symmetry, the chiral symmetry, which is realized by the gauge symmetry of the D8/anti-D8-branes. When L = πM −1 KK , the scale of chiral symmetry breaking coincides with the confinement scale. This is the choice of parameters that is useful to model QCD, and that was, indeed, considered in the original version of the model [19]. In the general case, L can be considered as a free parameter of the model. This latter case has been considered in the recently proposed Holographic QCD axion scenario [16,17].
In the regime (2.1), the gauge plus matter adjoint sector can be studied by considering the near-horizon limit of the backreaction of the D4-branes, that is a solution of Type IIA supergravity dual to a theory known as Witten-Yang-Mills (WYM). It is given by a curved metric, a non-trivial dilaton, and a four-form Ramond-Ramond (RR) field strength. If we consider the theory at finite temperature, the time direction is compactified on the circle t ∼ t + 1/T , and therefore we have two circles: S 1 t and S 1 x 4 . As a result, depending on the temperature, there are two competing solutions. Let us briefly present them, working with the Euclidean signature.
The background that dominates in the high-temperatures regime is the black hole one, where ω 4 is the four-sphere volume form and g s , l s are the string coupling and string length. The parameter u T is related to the Hawking temperature T by The background that dominates in the low-temperature regime is called "solitonic" and reads with the dilaton and F 4 fields keeping precisely the same form as in the previous case. The holographic dictionary that relates the string theory quantities and the field theory ones reads It can be shown that, since in the first case g 00 (u T ) = 0 and in the second one g 00 (u 0 ) = 0, the high-temperature solution is dual to the deconfined phase and the low-temperature one is dual to the confined phase. By computing the free energy of the two backgrounds, one finds that the system exhibits a first-order phase transition at temperature T c = M KK /2π.
Let us consider now the fundamental matter sector. In the regime (2.1), the backreaction of the D8/anti-D8-branes can be neglected. As a result, they can be treated in the probe approximation, namely by means of the Dirac-Born-Infeld action for the D8 branes on the original backgrounds. The embedding of the branes on the geometry will then be a solution x 4 = x 4 (u) of the equation of motion coming from this action and found by asking the D8/anti-D8-brane pair to be separated by a distance L on the circle S 1 x 4 , as mentioned above. Let us describe the solutions in the two phases and how they depend on the distance L.
In the confined phase, each D8/anti-D8 branes pair is actually bound to join into a single U-shaped configuration. From the field theory perspective, this fact is interpreted as a realization of chiral symmetry breaking. When L = πM −1 KK , the branes are antipodal and join at a value u J of the holographic coordinate that coincides with the smallest value of the coordinate range, that is u J = u 0 . This means that chiral symmetry breaking occurs at the confinement scale. In contrast, when L < πM −1 KK , the branes join at u J > u 0 , meaning that chiral symmetry breaking and confinement can occur at different scales. In the QCD-like setup, N f coincident pairs of D8/anti-D8-branes are placed in the antipodal configuration and the model realizes the breaking of U (N f )×U (N f ) to the diagonal U (N f ). At low energies, the effective action on the D8-branes reproduces the chiral Lagrangian (with pion decay constant f π ∼ √ N M KK ) including the Skyrme term and the axial anomaly term that gives mass to the η particle. The model has been generalized in [16] so that the effective Lagrangian also includes the Peccei-Quinn axion. 11 This can be easily obtained by considering an extra D8/anti-D8 pair placed in a non-antipodal configuration in order to achieve a separation between the axion scale f a and the QCD confinement scale. The U (1) gauge symmetry of the extra pair of branes is holographically interpreted as the Peccei-Quinn U (1) P Q global symmetry whose breaking gives the axion as a pseudo-Nambu-Goldstone boson.
In the deconfined phase, branes and anti-branes are not bound to join, because they can terminate on the horizon. As a result, depending on L, there are two possible D8-brane embeddings. If L > 0.97M −1 KK , the branes remain disconnected: the embedding x 4 = x 4 (u) reduces to a constant. From the field theory side, this corresponds to chiral symmetry restoration. In contrast, if L < 0.97M −1 KK , both the connected and the disconnected embeddings are allowed and, depending on temperature, only one is energetically favored. As a result, the model features a further first-order transition, occurring at a critical temperature T χ c different from the confinement/deconfinement critical temperature T c and given by T χ c 0.1538/L. For T > T χ c , the disconnected solution is preferred and thus chiral symmetry is restored, while for T < T χ c , the connected one is favored and thus chiral symmetry is broken.

B Calculation of the gravitational wave spectra
In this appendix, we review all the formulae needed to calculate the gravitational wave spectra produced by cosmological first-order phase transitions. The formulae for the spectra are reported in section B.2. They require the knowledge of some crucial parameters which we discuss in section B.1. These are essentially given by the temperature (and hence the related value of the Hubble parameter) at which the phase transition completes, the phase transition duration β −1 , computed starting from the bubble nucleation rate Γ, the strength α, i.e., the energy budget of the transition and the bubble wall speed v.

B.1.1 Bubble nucleation rate
First-order phase transitions are triggered by the nucleation of true vacuum bubbles on the false vacuum state. Such nucleation can occur through thermal or quantum fluctuations. As we will discuss in the following subsections, whether the transition actually takes place depends on the ratio Γ/H 4 , where Γ is the bubble nucleation rate per unit of volume and H is the Hubble scale. The latter is determined by the energy density ρ through the Friedmann equation H 2 = ρ/3M 2 P l , where M P l ≈ 2.4 · 10 18 GeV. The bubble nucleation rate can be computed using the well-known formalism developed in [2,5,6] for models where the transition is described by a single field Φ. One has to find a particular solution Φ B of the Euclidean equation of motion usually called bounce. The latter satisfies the following boundary conditions: it approaches the false vacuum Φ f at Euclidean infinity and a constant Φ 0 at the center of the bubble. 12 When the transition from the false to the true vacuum is due to quantum tunneling, the bounce is O(4) symmetric: in this case Φ B only depends on the radial coordinate ρ = √ t 2 + x i x i , where t is the Euclidean time and x i are the space coordinates. When the transition is (mostly) driven by thermal fluctuations, the bounce is O(3) symmetric: in this case Φ B = Φ B (ρ), with ρ = √ x i x i . The configuration which dominates the process is the one for which the rate Γ has the larger value. As a result, the formula for the bubble nucleation rate reads where ρ w is the size of the O(4) bubble. The bounce action S 3 appearing in (B.1) is defined by Euclidean action for the scalar field. The action S 4 is defined analogously.

B.1.2 The relevant temperatures
In order to calculate the spectrum of GWs, the first datum to determine is the temperature at which the waves are produced. Since from the time of nucleation, which happens at plasma temperature T n , to the time where most of the collisions take place and most of the sound waves collide there could be a sizable difference, the percolation temperature T p is considered to be the relevant one for the production of gravitational waves [9]. In the following, we are going to discuss both T n and T p .

Nucleation temperature
The nucleation time t n is defined as the time at which the total number of nucleated bubbles per Hubble patch from t = t c (the time when the Universe is at the critical temperature T c ) to t = t n is order one, where f (T ) is a polynomial function, usually assumed to be T 4 from dimensional analysis. Let us write the Taylor expansion of the exponent as Thus, the condition (B.4) can be approximately computed as where we extended the integration domain to infinity, and therefore it reads Percolation temperature The percolation temperature T p is defined as the Universe temperature when the fraction of space sitting in the true vacuum takes a benchmark conventional value. We choose the latter to be one. 14 In order to compute the percolation temperature, we have to estimate the size of a bubble as a function of time, which involves the knowledge of the bubble wall speed v. We follow [46]. The fraction of space in the true vacuum reads where r b (t, t ) is the size of the bubble in comoving coordinates as a function of time, which can be obtained by Here (B.12) We therefore define the percolation temperature T p by In the scenarios considered in this paper, the energy density includes a radiation term and a vacuum term. Adopting the notation of [46], we therefore write H = H R + H V . Let us consider approximate solutions to (B.12). Firstly, let us consider the case in which the vacuum contribution H V can be neglected. This is expected to give a good approximation when supercooling is not significant. Assuming H R = c R T s , and constant velocity v, we obtain (B.14) The formulae of reference [46] are retrieved putting s = 2, r = 1, and c R = ( √ 3M P l ξ g ) −1 . In the WSS model (see (2.9)), space occupied by the false vacuum, V f alse ∝ R(t) 3 exp (−I RV ), is decreasing at the supposed percolation temperature. This translates into the condition Once the percolation temperature has been determined, one can derive a crucial parameter for the spectrum

Reheating temperature
During a first-order phase transition, entropy is released and therefore the Universe gets heated. Assuming that the entropy release is approximately instantaneous, we define the reheating temperature T R as the temperature of the Universe after the release. By exploiting the conservation of the energy density during the transition, we find the reheating temperature through the formula where ρ f and ρ t are, respectively, the total energy density in the false and true vacua. Especially in the case of strong first order transitions, like the ones examined in this paper, the reheating temperature may be greater than the critical temperature T c (see e.g. [53]).
In this case, one should check whether the inverse phase transition could take place or not. In the cases examined in this paper this does not happen essentially because the distance in field space between the two minima of the effective potential at T = T R T c is "large" enough to drastically suppress the rate of the inverse transition w.r.t. the Hubble scale.

B.1.3 Released energy and wall velocity
Another crucial parameter for the gravitational wave spectra is the ratio of the energy released in the transition to the energy of the radiation bath [9]. In particular, the formulae for the spectra include the parameter α defined as where θ = (ρ − 3p)/4 is the trace of the energy-momentum tensor, and the ∆ indicates the difference between the false and true vacua. The knowledge of the parameter α allows us to estimate the velocity of the bubble walls according to the Chapman-Jouguet formula Formula (B.25) has a limited range of validity; in particular, it has to be corrected when the friction in the bubble interactions with the plasma is significant. However, to provide better estimates of the wall speed is still one of the big open problems in determining the bubble dynamics.
In the case in which we consider a dark sector that is not in thermal equilibrium with the visible one, we have to define two separated α parameters for the two sectors, which take into account the fact that the relevant radiation could be only the one of the visible (dark) sector, ρ rad,SM (ρ rad,glue ). Note that if the Standard Model plasma is not interacting with the dark sector one, in formula (B.25) one has to replace α with α D .

B.1.4 Redshift
Once we know the parameters that we have discussed so far, we can compute the spectrum of gravitational waves as it appears at the time of production. From this time to the time of detection, the signal gets redshifted due to the cosmological expansion. We are going to discuss how to take into account the redshift of the signal in two different circumstances, both occurring in the scenarios that we study in the present paper. Let us start with the case in which the Universe evolves adiabatically from gravitational waves emission to the detection time [54]. This is the case in which only one first-order phase transition occurs. Hence, it includes all the scenarios that we consider in this paper but the Dark HQCD 2 one. Let us call T e and T d the temperature of the Universe, respectively, at the emission and at the detection times. The detection temperature is T d ∼ 2.35 · 10 −13 GeV. The adiabatic evolution is characterized by the conservation of the entropy from which we find the ratio of the scale factors between the two temperatures In this expression, g S * ,e and g S * ,d are the number of relativistic degrees of freedom at the time of emission and detection, respectively; they are computed in the free case using the general formula where T i represents the temperature of the i-th species.
The frequency f and the energy density 15 Ω of the GWs get redshifted as R −1 and R −4 respectively, hence . (B.30b) The Hubble scale H is given by the energy density via the Friedmann equation Here, g * is defined in the free case as 16 As a result, using (B.28) and (B.31), we find Let us now consider the case in which two first-order phase transitions occur. Among the scenarios studied in this paper, this happens in the Dark HQCD 2 scenario of section 4.1, where a chiral symmetry breaking/restoration transition is followed by a confinement/deconfinement one. We will refer to this case, even though the discussion will be valid for two generic separated first-order phase transitions.
When we compute the redshift of the gravitational waves spectrum associated with the chiral symmetry transition, we have to take into account that adiabaticity is violated during the confinement/deconfinement one. As a result, conservation of entropy can be used from the time of the chiral symmetry breaking transition, where the temperature T e is taken to be the reheating temperature, to the percolation time of the confinement transition T p,conf . Then, assuming fast reheating in the confinement transition, the temperature has a sudden 15 As customary in cosmology, Ω is defined as the energy density divided by the critical density ρ crit,0 = 3M 2 P l H 2 0 , where H 0 is the Hubble scale computed in the present epoch. 16 Notice that if some species are decoupled from the bath, g * = g S * . This, notoriously, occurs in the cosmological evolution because of neutrino decoupling when electrons and positrons become non-relativistic. Neutrino and photon relic temperatures do not coincide. They are related by T ν = (4/11) 1/3 T γ . As a result, today g * ≈ 3.36 is different from g S * = 3.91.
jump from the percolation temperature T p,conf to the reheating temperature T R,conf . Finally, from this time to the present, the Universe evolves adiabatically, and we can again use the conservation of entropy. All in all, the redshifted frequency and energy density read and With respect to the single-transition case, the difference is encoded in the parameter In models with multiple, separated phase transitions, a δ factor for each transition after the first one must be included in the formulae for the GW spectra.

B.2 Formulae for the spectra
Let us finally discuss the formulae that allow us to find the gravitational wave spectrum. In linear approximation, the spectrum is given by the sum of three contributions, coming from the collisions of the bubbles, from collisions of plasma sound waves and plasma turbulence, Here, h is defined from today's value of the Hubble scale through H 0 = 100 hKm/s/Mpc. Following [9], we are going to neglect the turbulence contribution because it is still not wellunderstood and because it is expected that only a small fraction of the transition energy is converted to turbulence. Let us first consider the collision contribution. Using the so-called envelope approximation, a formula for the signal of gravitational waves coming from bubble collisions was numerically found in [54]. An improved version of such a formula (see, e.g., the review [45]) reads where f denotes the frequency of the waves, and v the average velocity of the bubbles. The factor κ quantifies the fraction of available energy converted into gravitational waves coming from bubble collision. Finally, the spectral form S env and the peak frequency f env are given where we are going to use the approximate value z p ∼ 10, and the efficiency factor in the case v ∼ 1 is If the Standard Model plasma is not in thermal equilibrium with the dark sector, in this formula one has to use α D instead of α [40]. In fact, formula (B.39a) is valid under the assumption that the source of GWs lasts for a period longer than a Hubble time. If the source's duration is short, turbulence effects can be sizable and one can estimate that the net effect is to multiply formula (B.39a) by a factor [9] (8π This term tends to reduce the amplitude of the signal. On the other hand, one should add the contribution due to turbulence, which is very uncertain and as stated it is ignored in this paper. Thus, Ω sw including the term (B.41) really corresponds to a lower bound on the contribution of the plasma to the GW spectrum.

C Holographic bubbles
In this appendix, we review the results of [20] that are used in order to compute the spectrum of gravitational waves. In [20], the confinement/deconfinement phase transition occurring in the WSS model was studied using an effective approach inspired by [27], and deployed in order to reduce the problem to a single-scalar one (a recent alternative approximation can be found in [58]). Moreover, the chiral symmetry breaking/restoration transition was described by the deformation of the D8 embedding, again encoded in a single scalar mode. We start with the case of the confinement/deconfinement phase transition and then we discuss the chiral symmetry breaking/restoration one.

C.1 Confinement/deconfinement phase transition
The main idea for studying the confinement/deconfinement phase transition is to promote the parameters u 0 and u T in (A.1) and (A.3) to fields depending on a radial coordinate ρ in such a way that the background displays a conical singularity. For the case in which the transition is driven by thermal fluctuations and the bounce is O(3) symmetric, this can be achieved by taking the ansatz Once the solution is found, one can plug it back into the action. As we have already pointed out, what really enters the formula for the nucleation rate is the difference between the on-shell action on the bounce solution and the action evaluated on the false vacuum, (C.10) From the numerical results and the functional form of the thin and thick wall approximations studied in Appendix A of [20], a continuous analytic approximation to the action for the O(3) bubble can be provided as follows, For small temperatures, one could also have O(4) symmetric bounces. The action for this case is simply (C.6) generalized such that it enjoys O(4)-symmetry, From this action, we find the bounce solution imposing the same boundary conditions as above. For the O(4) bubble, since it is only defined for small temperatures, it is sufficient to consider the functional form of the thick wall approximation, giving Let us recall that the O(4) configuration can be admitted only if the bubble radius is much smaller than 1/T , [5,6]. In [20] the following convention has been adopted: the maximal allowed radius for the O(4) configuration to be considered is set by ρ w = 1/2πT (the radius of the thermal circle). In the present setup, this corresponds toT ≈ 0.06, which explains the condition in parenthesis in (C.13).

C.2 Chiral symmetry breaking/restoration phase transition
As we have reviewed in section 2 and appendix A, in the regime M KK 2π < 0.1538 L (C.14) the WSS model displays a first-order phase transition associated with chiral symmetry breaking in the deconfined phase. In the probe regime N f N , the phase transition can be studied just by considering the Dirac-Born-Infeld action on the fixed black hole background (A.1) which, using spherical coordinates for the 3d Euclidean physical space reads Considering an ansatz in which the embedding of the D8-branes in the above background is described by a function x 4 (u, ρ), the DBI action reads S DBI = T 8 g s d 9 xρ 2 u R −3/2 where T 8 is the D8-brane tension. Let us now rescale the coordinates as follows. First, define , u = y u T , u J = y J u T , (C. 17) so that the periodicity of the cigar coordinate now reads Then, let us define . (C.20) Using the definitions above, the DBI action can be rewritten as y 5/2 1 + (y 3 − 1)(∂ y x) 2 + (∂ σ x) 2 − 1 dy − 2 7 (y J (σ) 7/2 − 1) , (C.23) where we have subtracted the contribution of the straight brane/antibrane pair configuration. The Euler-Lagrange equation for x(y, σ) reads ∂ y σ 2 y 5/2 (y 3 − 1)(∂ y x) 1 + (y 3 − 1)(∂ y x) 2 + (∂ σ x) 2 + ∂ σ σ 2 y 5/2 (∂ σ x) 1 + (y 3 − 1)(∂ y x) 2 + (∂ σ x) 2 = 0 . (C. 24) This is a non-linear partial differential equation, which is extremely hard to solve even numerically. An escape strategy has then been put forward in [20]: it amounts to look for approximate solutions by using a reasonable variational ansatz,  ) where ρ is the 4d Euclidean radial coordinate is not consistent with the equations of motion. In [20], just as a reference, we have considered a "naive O(4) configuration" obtained by simply considering the measure d 4 x to be given by dΩ 3 dρρ 3 , where dΩ 3 is the measure of the three-sphere. We have not found convincing indications that a real O(4) configuration can actually be achieved for the chiral symmetry breaking transition. Thus, in this case, we have decided to focus only on the O(3) symmetric one.

C.3 GW parameters
It is instructive to give an estimate of how the relevant parameters entering the computation of the stochastic GW spectra depend on the WSS parameters. In this subsection we focus on the confinement/deconfinement phase transition and neglect the flavor contributions.
In the small temperature regime, when the O(4) symmetric bounce dominates the transition, the bubble nucleation rate (B.1) can be easily computed using the relations (C.13) giving Γ(T ) = M 4 Again,T n decreases as g and M KK increase.
In the same limits ad before, the percolation temperature approximately coincides with the nucleation temperature and where g * SM = g * SM (T R ) = O(100). The reheating temperature is thus independent from M KK and increases with g. When g 1, it is parametrically larger than the critical temperature. Finally, using the definition (B.24), it is possible to estimate the parameter α, measuring the relative energy released during the transition. In the small temperature regime it reads α(T p ) ≈ 1 5T 6 p 1 . (C.43)