Gravitational Wave Signals from Multiple Hidden Sectors

We explore the possibility of detecting gravitational waves generated by first order phase transitions in multiple dark sectors. Nnaturalness is taken as a sample model that features multiple additional sectors, many of which undergo phase transitions that produce gravitational waves. We examine the cosmological history of this framework and determine the gravitational wave profiles generated. These profiles are checked against projections of next generation gravitational wave experiments, demonstrating that multiple hidden sectors can indeed produce unique gravitational wave signatures that will be probed by these future experiments.

In this work, we explore the possibility of having multiple decoupled hidden sectors. Large numbers of hidden sectors can solve the hierarchy problem as in the Dvali Redi model [37], in the more recently explored N naturalness [38] framework, or in orbifold Higgs models [39,40]. They can also be motivated by dark matter considerations [22,41,42]. Motivated by solutions to the hierachy problem, we consider hidden sectors with the same particle content as the Standard Model that have all dimensionless couplings (defined at some high scale) equal to those of the Standard Model. The only parameter that varies across sectors is the dimension two Higgs mass squared parameter, m 2 H . This simple ansatz can lead to very rich phenomenology and interesting gravitational wave spectra, but we stress that it is only a starting point for exploring the space of theories with multiple hidden sectors.
In this setup, there are two qualitatively different kinds of sectors: • Standard Sectors: Those with m 2 H < 0 where electroweak symmetry is broken by the vacuum expectation value (vev) of a fundamental scalar. As * Paul.Smith3@carleton.ca Dylan.linthorne@carleton.ca stolar@physics.carleton.ca in [38], we assume that the standard sector with the smallest absolute value of m 2 H is the Standard Model.
• Exotic Sectors: Those with m 2 H > 0. In this case, electroweak symmetry is preserved below the mass of the Higgs, and broken by the confinement of QCD [43].
Cosmological observations, particularly limits on extra relativistic degrees of freedom at the time of Big Bang Nucleosynthesis and the time of the formation of the cosmic microwave background (CMB) [44], require that most of the energy in the universe is in the Standard Model sector as we will quantify. Therefore, the hidden sectors cannot be in thermal equilibrium at any time, and the physics of reheating must dump energy preferentially in the Standard Model sector. This can be accomplished with primodial axion-like particle (ALP) models [45,46] and with the reheaton method [38]. We will also explore alternative parameterizations of reheating that satisfy this condition.
In all the above models, there is some energy in the hidden sectors, and these sectors undergo thermal evolution independent of the SM sector. If their initial reheating temperature is above their weak scale, the standard sectors will undergo phase transitions associated with the breaking of electroweak symmetry and with confinement of QCD. The exotic sectors will also undergo a phase transition when QCD confines and electroweak symmetry is broken simultaneously. The condition for these transitions to leave imprints on the stochastic gravitational wave spectrum is that they strongly first order phase transitions (SFOPT) [12][13][14][15]. This does not occur at either the electroweak or QCD phase transition in the SM, but as we will show, it does happen for the QCD phase transition in some standard sectors and in all exotic sectors that reheat above the QCD phase transition.
This work is organized as follows: section II introduces the particle content of the model, section III discusses the phase transition behaviour of the both the standard and exotic sectors present, section IV lays out hidden sector reheating, section V applies constraints from cosmological observables allowing for the calculation of gravitational wave signatures in section VI, and, finally, section VII ties everything up.

II. PARTICLE SETUP
We consider the following Lagrangian as in [38]: with L 0 = L SM being the Standard Model Lagrangian, and L i being a copy of the SM Lagrangian with different fields, but with all dimensionless parameters the same. Each of the Lagrangians does contain a dimensionful operator: where H i is a Higgs field in each sector, and the mass term is parametrically given by where Λ is some high scale cutoff, N is the number of sectors, and r is the mass parameter in the SM in units of Λ 2 H /N . We view the parameterization of Eq. (3) as a random distribution in theory space up to the cutoff Λ, therefore this setup solves the hierarchy problem if r ∼ O(1) [38] 1 and our sector is the one that that has the smallest absolute value of the Higgs mass parameter. We have taken for simplicity that there equal numbers sectors with positive and negative m 2 H , but this assumption does not affect our analysis. This N naturalness framework can be generalized: the various sectors can possess a wide range of particle content that can be freely selected by the model builder. The one exception to this is that "our" sector must consist of the Standard Model.
From the above Lagrangians, the Higgs in sectors with i ≥ 0 will get a vev given by: λ i is the quartic coefficient of the scalar potential and is the same across all sectors, λ i = λ. This is another way to see how this framework can solve the hierarchy problem: the Higgs vev is parametrically smaller than the cutoff for N 1. The "standard sectors" with i > 0 feature electroweak symmetry breaking just like in the SM, however the vevs scale with the changing mass parameter: v i ∼ v SM √ i. This means that the masses of the fermions and the W and Z will also increase proportional to √ i. The consequences of this scaling on the 1 Constraints require r to be somewhat smaller than 1.
confinement scale of QCD in the i ≥ 1 sectors is further discussed in Sec. III. The "exotic sectors" with i < 0 provide a radical departure from our own. m 2 H > 0 leads to no vev for the Higgs, and electroweak symmetry is only broken at very low scales due to the phase transition from free quarks to confinement at the QCD scale Λ QCD [43], and the masses of the W and Z are comparable to those of QCD resonances. The mass of fundamental fermions are produced via four-fermion interactions generated after integrating out the SU(2) Higgs multiplet. This leads to very light fermions: with y f representing the Yukawa coupling to fermion f . As we will see, the extremely light quarks that appear in these sectors dramatically change the nature of the QCD phase transition -unlike the SM, the transition is strongly first order. Again, this is further developed in Sec. III. Crucially, this results in the production of gravitational waves. This is the physical signature we explore in this paper; the calculation and results are presented in Sec. VI.

III. QCD PHASE TRANSITION
We now study the nature of the QCD phase transition across the different sectors. Due to the confining nature of QCD, the exact nature of the phase transition is often difficult to ascertain analytically and requires the study of lattice simulations. In the SM, it is known that the phase transition is a crossover and does lead to gravitational wave signals [47,48]. In the general case with 3 or more colours, the phase transition can be strongly first order in two regimes [49][50][51]: • Three or more light flavours.
Light indicates a mass small compared to the confinement scale Λ QCD , but what that means quantitatively is not precisely determined. In the SM, the up and down quarks are light, but the strange is not sufficiently light for an SFOPT. For the standard sectors in our setup, the quark masses increase with increasing vev, so for sufficiently large i, all the quarks will be heavier than Λ QCD , 2 and those large i sectors will undergo an SFOPT if they are reheated above the the confinement scale. Conversely, exotic sectors with zero vev feature six very light quarks, so all the exotic sectors undergo SFOPT at the temperature of QCD confinement.
We now calculate the QCD confinement scale for each sector following the same procedure as [52]. First, due to the parameters of each sector being taken to be identical save for the Higgs mass squared (thus v = v i , where v is the SM vev), we assume that the strong coupling of every sector is identical at some high scale. Using the one-loop running, the β function can be solved: where n i f is the number of quark flavours with mass less than µ/2 and Λ i is the scale where it would confine if all quarks remain massless. In the SM defined at scales well above all the quark masses, we have Λ QCD = 89 ± 5 in M S [53]. Because we have set the strong couplings equal at high scales, Λ = Λ i for all i at high scales for all sectors. However, since the masses of the quarks in each sector are different, we end up with a unique running of the coupling for each sector. At every quark mass threshold for a given sector, we match the coupling strengths above and below the threshold and determine the new Λ i for the lower scale. For example, at the mass of the top quark, we match a five flavour coupling with the six flavour one: and thus Suppressing the i's for notational cleanliness, we can arrive at similar relations at the bottom and charm thresholds Λ (4) = (m b ) 2/25 (Λ (5) ) 23/25 , These can be combined to show that This type of matching procedure can be done as many times as necessary for a given sector. The process terminates when Λ i for a given scale is larger than the next quark mass threshold (i.e running the scale down arrives at the Λ QCD phase transition before reaching the next quark mass scale). In cosmological terms, we can envision a sector's thermal history unfolding, whereas the plasma cools below each quark mass threshold and said quarks are frozen out. At a certain point, the sector arrives at the QCD phase transition and confinement occurs -if this occurs when ≥ 3 quarks are at a much lower scale or all quarks have already frozen out, we get the desired phase transition.

A. Standard Sectors
As shown in Eq. (4), for standard sectors with increasing index i the vevs of said sectors increase v i ∝ √ i. This leads to increasingly heavy particle spectra for higher sectors -eventually leading to sectors that are essentially pure Yang-Mills that feature strong first order phase transitions. This, of course, prompts the question: at what index i do said phase transitions begin? Using the methods outlined in the prior section we determine Λ QCD to have a relevant value of at the energy scale we're interested in. Λ i (6) is identical for all sectors and is taken to have a standard model value of Λ (6) M S = (89 ± 6) MeV [53]. Rewriting Eq. (11) in terms of standard model variables, where m q without a superscript is the mass of q in the SM. We take the sector with SFOPT to be the ones when the mass of the up quark, down quark, and QCD phase transition scale are all comparable: This can be solved for i: As we will see in Sec. IV, in the original N naturalness setup [38], the energy dumped into the ith sector scales as i −1 , so there will will not be enough energy in the sectors with i > i c to see a signature of these phase transitions. However, if we move away from the original N naturalness reheating mechanism and begin exploring mirror sectors with large vevs and with relative energy densities ρ i /ρ SM ∼ 10%, a possibility allowed by current constraints, we can have sectors with relatively high dark QCD scales that produce detectable gravitational waves. From Eq. (11) we can determine the confinement scale of an arbitrary mirror sector. If we take Higgs vevs as high as the GUT scale ∼ 10 16 GeV, then we can use Eq. (12) to get confinement scales as high as ∼ 38 GeV . The signals of this sector and other test cases like it are explored in Sec. VI.

B. Exotic Sectors
In every exotic sector the fermion masses are exceptionally light: their masses are generated by dimension six operators with the Higgs integrated out as shown in Eq. (5), and are therefore all below the confinement scale. The exotic sectors all have identical oneloop running of the QCD gauge coupling, and thus all have approximately the same confinement scale given by Λ ex ∼ 90 MeV. These sectors all have six light fermions, so a strong first order phase transition occurs for all exotic sectors at this temperature. The confinement of these sectors directly leads to the production of both baryons and mesons as we have the spontaneous breaking of SU(6) × SU(6) → SU(6) and thus 35 pseudo-Goldstone bosons (pions). The masses obtained through the phase transition can be approximated through the use of a generalization of the Gell-Mann-Oakes-Renner relation [54,55], where V ∼ Λ QCD , F π is the pion decay constant. One expects that within a given sector F π ∼ V ∼ Λ QCD [55] and as exotic sectors have Λ ex ∼ 90 MeV while the SM features Λ QCD = (332 ± 17) MeV [53] we expect at most coefficient relative to the SM value. So, for pions in exotic sector i: Here, a and b denote the component quark flavours.

IV. REHEATING N SECTORS
A key issue within N naturalness is how to predominantly gift energy density to our own sector so as to not be immediately excluded by cosmological constraints, particularly those from number of effective neutrinos (N ef f ). Here we review the results of [38]. Reheating occurs through the introduction a "reheaton" field. After inflation, the reheaton field possess the majority of the energy density of the Universe. Although this field can generically be either bosonic or fermionic, we reduce our scope to a scalar reheaton φ. Our focus is primarily the production of gravitational waves from multiple sectors and a fermion reheaton does not change the scaling of the energy density of the exotic sectors and thus does not affect expected gravitational wave profiles.
In order to maintain the naturalness of our SM sector, the reheaton coupling is taken to be universal to every sector's Higgs. However, a large amount of the Universe's energy density must ultimately be deposited in our own sector for N naturalness to avoid instant exclusion. In order to accomplish this, the decay width of the reheaton into each sector must drop as |m H | grows. If we insist that the reheaton is a gauge singlet that is both the dominant coupling to every sector's Higgs and lighter than the naturalness cutoff Λ H / √ N , then we construct a model that behaves as desired. The appropriate Lagrangian for a scalar reheaton φ is: Note that cross-quartic couplings of the form κ|H i | 2 |H j | 2 that could potentially ruin the spectrum of N naturalness are absent, taken to be suppressed by a very small coupling. Effective Lagrangians for the two different types of sectors present in this theory can be obtained by integrating out of the Higgs bosons in every sector: with C i representing numerical coefficients, g the weak coupling constant, and W µν the SU(2) field strength tensor. Immediately from Eq. (18), we can see that the matrix element for decays into standard sectors is inversely proportional to that sectors Higgs mass, The loop decay of φ → γγ is always sub-leading and can be neglected. It should be noted that as one goes to sectors with larger and larger vevs, the increasing mass of the fermions (m f ∼ v i ∼ v SM √ i) eventually leads to situations where the decay to two on-shell bottom or charm quarks is kinematically forbidden, m φ < 2m q . For sectors where this kinematic threshold is passed for charm quarks, the amount of energy in these sectors becomes so small that contributions to cosmological observables can be safely ignored. All in all, we end up with a decay width that scales as Since we can expect energy density to be proportional to the decay width, ρi ρ SM ≈ Γi Γ SM , this indicates that energy density of standard sectors falls: with r s being the ratio of the energy density of the first additional standard sector over the energy density of our sector. For the exotic sectors, Eq. (18) indicates a matrix element scaling M m 2 H >0 ∼ 1/m 2 Hi and is also loop suppressed. This leads to a significantly lower energy density than the standard sectors. Both the decay width and energy density for these sectors scale as: As a final note, in this setup the reheating temperature of the SM, T RH , has an upper bound on the order of the weak scale. If this bound is not observed, the SM Higgs mass would have large thermal corrections -leading to the branching ratios into other sectors being problematically large [38]. Thus we only consider relatively low reheating temperatures 100 GeV.
Ultimately, after examining the gravitational wave case produced by standard N naturalness, we also consider a more generic parameterization where the reheating temperature of each sector is a free parameter and is in general uncorrelated with the Higgs mass parameter. This allows us to explore a broader model space with multiple dark sectors at a huge range of scales. For these models, the reheating mechanism remains unspecified.

V. CONSTRAINTS
In general, the multi-hidden sector models explored feature a huge number of (nearly) massless degrees of freedom. Dark photons and dark neutrinos abound in these sectors and, assuming a relatively high reheat temperature, the leptons, quarks, and heavy bosons of these sectors can also be relativistic. In N naturalness this feature is realized quite dramatically: each of the N sectors possess relativistic degrees of freedom. The presence of these particles can have two main effects: extra relativistic particles can alter the expansion history of the universe through changes to the energy density or hidden sectors can feature annihilations that reheat the photons or neutrinos of our sector near Big Bang Nucleosynthesis (BBN) and affect the light element abundances. The effective number of neutrino species, N ef f , is impacted by these contributions and, as such, is the strictest constraint that must be dealt with when studying these type of multi-phase transition models. The SM predicts that N SM ef f = 3.046 [56]. This is in good agreement with the 2σ bounds from studies of the Cosmic Microwave Background (CMB) by Planck combined with baryon acoustic oscillations (BAO) measurements [44]: Various different assumptions about the history of the universe can be made and different data sets can be chosen to obtain slightly different results [32] -for the purposes of this exploratory work, wading through this landscape is unnecessary. Additionally, for any decoupled hidden sector [38]. Because the constraints on N ef f are stronger at photon decoupling than at BBN, we can focus purely on the constraints provided by the former. Future CMB experiments [57] will improve the bound from Eq. (21) by about an order of magnitude. This could significantly reduce the allowed temperature ratio of any hidden sector, or alternatively could provide evidence for such sectors in a way that is complementary to the gravitational wave signatures described below. For fully decoupled sectors that never enter (or reenter) thermal equilibrium with our sector, we obtain additional contributions to N SM ef f [32] ∆N ef f = 4 7 Here, g h represents the effective number of relativistic degrees of freedom for the hidden sector 3 , and we parameterize the hidden sector temperature by [32] ξ 3 g h = N boson + 7N fermion /8. and these should be evaluated at the time of photon decoupling. We take this approach and generalize it to include many additional sectors: For a dark sector with one relativistic degree of freedom, its temperature must be T DS ∼ 0.6 T SM to not be excluded. Applying the energy density formula [58], to both said dark sector and the SM then taking the ratio indicates that the dark sector would have an energy density ρ ∼ 0.038 ρ SM .

A. Exotic Sector Contributions
We begin by computing the constraints on exotic sectors; these are significantly weaker than those for standard sectors [38]. At the time of photon decoupling, T γ ∼ 0.39 eV while the temperature of the exotic sectors is lower. This means that for sectors with small and moderate i, we can use Eqs. (5) and (16) to see that the pions will be non-relativistic leaving at most 11.25 effective degrees of freedom per sector from photons and neutrinos. For very large i, the pions can be much lighter, but those sectors also have very little energy in them in the standard reheating scenario. Coupling the number of effective degrees of freedom per sector with the energy density scaling of ∼ 1/m 4 H as in Eq. (20) means that the zero vev sectors have small temperature ratios. Assuming a reheating temperature of 100 GeV and a completely uniform distribution of sectors, the temperature of the first exotic sector is slightly more than 6% of our sector at reheating. Applying Eq. (25) to this particular situation gives us: with T RH E1 /T RH being the ratio of the reheat temperatures of the 1st exotic sector and our own sector (0.06 in standard N naturalness with r = 1). This sum is dominated by i = 1, the sector with the lowest Higgs mass (and thus the most energy density) gives us a contribution of O(10 −4 ) to ∆N ef f . Evolving the sectors thermal histories forward in time to the recombination era gives us a slightly larger value, but still of order O(10 −4 ), well below current CMB bounds. It should be noted that modifying the exotic sectors' structure (e.g. adjusting the exotic sectors to have a lower Higgs mass squared or clustering multiple hidden sectors close to the first exotic one) leads to a ∆N ef f contribution that is larger than the base N naturalness case. This increase is typically not excluded by current bounds, indicating a large degree of liberty in the structure and number of exotic hidden sectors.

B. Standard Sector Contributions
Within the context of vanilla N naturalness, the majority of contributions arise from standard sectors. This is explored in detail in [38]; here we briefly summarize these arguments. All additional standard sectors are very similar to our own: they have the same particle content and couplings and differ only by the Higgs mass. As our sector is taken to be the lightest so as to be preferentially reheated, every other standard sector features an earlier freeze-out of their respective particles. This ultimately leads to each sector having at most the same number of relativistic degrees of freedom as the SM.
In [38], the standard sector contributions are expressed as: In the case that the reheaton is lighter than the lightest Higgs (ours), this can be expressed as with y c,b representing the charm and bottom yukawa couplings and with m φ being the mass of the reheaton. Application of these results indicates that for a majority of the parameter space, vanilla N naturalness requires mild fine-tuning (r in Eq. (4) set to a value 1). Numerical results for the fine-tuning required for various reheaton masses were presented in [38].

C. Generalized Reheating Scenarios
The generalization of possible reheating mechanisms mentioned in section IV -where the reheating mechanism no longer depends on the Higgs' mass parameter of a given sector -opens up a wide range of hidden sectors for study. Specifically, this allows mirror sectors with large Higgs vevs to be reheated to significant energy densities and thus produce gravitational waves with enough power to be detected. Crucially, despite this analysis being limited to mirror sectors with large Higgs masses, this analysis pertains to any strong, confining phase transition at high scales. Since N ef f constraints remain our strongest cosmological bounds for massive standard sectors, our starting point for exploring the limits of high transition temperatures is Eq. (25). Assuming heavy, standard sectors (with the only relativistic particles being photons and neutrinos) we can saturate the bounds of Eq. (21) and solve for the maximum temperature allowed for any number of sectors: T i ∼ 0.38 T SM 1 hidden sector, T i ∼ 0.25 T SM 5 hidden sectors, T i ∼ 0.21 T SM 10 hidden sectors, T i ∼ 0.12 T SM 100 hidden sectors, (31) where all the hidden sectors have the same temperature as one another.
Using these restrictions, we can examine the behaviour of standard sectors with a much larger vev than our own. In terms of the N naturalness framework, this means we can get an SFOPT for QCD if we look at sectors with i greater than the critical index of Eq. (14) where all the quark masses are above the QCD confinement scale, as long as their temperatures are below the bounds presented here.

VI. GRAVITATIONAL WAVE SIGNALS
We now turn to the gravitational wave signatures of our setup. At high temperatures, each of the hidden sectors has QCD in the quark/gluon phase, but at temperatures around Λ QCD,i , the i th sector undergoes a phase transition into the hadronic phase that we computed for the different sectors in Sec. III. As discussed in that section, this phase transition will be strongly first order (SFOPT) for certain numbers of light quarks, which will generate gravitational waves. This differs from QCD in the SM sector, as the PT is a crossover and not first order [59]. A SFOPT proceeds through bubble nucleation, where bubbles of the hadronic phase form in the vacuum of the quark phase. These bubbles will expand, eventually colliding and merging until the entire sector is within the new phase. These bubbles are described by the following Euclidean action [60]: where the time component has been integrated out due to nucleation occuring not in vaccum but in a finite temperature plasma. φ is the symmetry breaking scalar field with a non-zero vev. In the case of the chiral phase transition, the scalar field breaking the SU(N f ) R × SU(N f ) L chiral symmetry is the effective quark condensate φ i ∼ qq i of the respective sector. We leave the thermalized potential V (φ, T ) general. As previously stated, an exact QCD potential at the time of the chiral phase transition is not well understood outside of lattice results. In [31] chiral effective Lagrangian was used to calculate a low energy thermalized potential for confining SU(N ). The amount of energy density dumped into the individual sectors dictates the energy budget for the PT and hence for the gravitational waves. Assuming that the SM sector is radiation dominated, a quantity that characterizes the strength of the PT is the ratio of the latent heat of the phase transition, , to the energy density of radiation, at the time of nucleation [61], with being calculable from the scalar potential. Assuming that there is a negligible amount of energy being dumped back into the SM, which would cause significant reheating of ρ γ , the latent heat should correspond to the energy density of the hidden sector going through the PT. The parameter g * in the denominator of Eq. (33) is the number of relativistic degrees of freedom at the time of the phase transition, with contributions from species in both the visible and dark sectors. It has weak temperature dependence in a single sector, but when dealing with multiple hidden sectors, g * gains contributions from all N sector's relativistic degrees of freedom, weighted by their respective energy densities with ξ being the temperature ratio defined in Eq. (24). The bounds from effective number of neutrinos [44] mean that ξ i 1 for all i, so g * ≈ g * ,γ . In the case of dark QCD-like chiral phase transitions, the temperature of the phase transition is on the the order of the symmetry breaking scale of the respective sector, T i h ∼ O(Λ QCD,i ). The work of [31] calculated α with an effective chiral Lagrangian have found various upper bounds. We take the optimistic scenario where the numerator is bounded above by the symmetry breaking scale where T nuc γ is the temperature of the SM photon bath at the time of the phase transition. Another important parameter to characterize the phase transition is its inverse timescale, β [16]. The inverse timescale can be calculated using the action in Eq. (32), The ratio of β and the Hubble constant, at the time of nucleation, H controls the strength of the GW signal, Due to the lack of a general analytic QCD potential, it is not possible to use Eq. (37) to calculate β/H.
There are dimensional arguments [13,14] that predict β/H ∼ 4Log(M p /Λ QCD,i ), although these arguments make specific assumptions about the potential. In more recent work, some authors [31,34] have attempted to estimate it using first-order chiral effective theories and the Polyakov-Nambu-Jona-Lasinio (PNJL) models. These studies claim a large range of values with no consensus reached on the precise order of the scaled inverse timescale. Under these circumstances, we take the optimistic case in which β/H is O(10).

A. Production of Gravitational Waves
Gravitational waves are produced with contributions from different components of the SFOPT's evolution. It is commonplace to parameterize the spectral energy density in gravitational waves by [62] Ω where ρ c = 3H 2 /(8πG) is the critical energy density. The total gravitational wave signal is a linear combination of three leading contributions: Each components is scaled by its own unique efficiency factor, κ. The three leading order contributions to the GW power spectrum are as follows: • Scalar field contributions Ω φ : Caused by collisions of the bubble walls, the solutions being completely dependent on the scalar field configuration, with efficiency factor κ φ = 1 − α ∞ /α [63,64].
• Magnetohydrodynamical contributions Ω B : Turbulence within the plasma, left over from the sound wave propagation, will produce gravitational waves with efficiency factor κ turb ≈ 0.1κ v [66].
The parameter α ∞ denotes the dividing line between the runaway regime (α > α ∞ ) and the non-runaway regime (α < α ∞ ). Explicitly [16,32,61], for particles with n i degrees of freedom that obtain mass through the phase transition.
The exotic sectors have essentially massless degrees of freedom pre-phase transition and pions with negligible masses post-phase transition. Other composite particles, such as baryons, do gain a mass of the order of Λ ex ; this is, however, still much smaller than the order of ρ R leading to small α ∞ according to Eq. (40). 4 Heavy standard sectors that undergo SFOPT for QCD feature no baryons due to all quarks being above their respective QCD scales. They do, however, feature glueballs that obtain a mass of the order of the SFOPT and, just as in the case for exotic sectors above, feature small α ∞ . Ultimately, calculating the critical phase transition strength for each sector using Eq. (40) indicates that every sector has a phase transition in the runaway regime. Each component of the spectral energy density in Eq. (38) is proportional to a power of their respective efficiency factors κ. In this case, the efficiency factors for the sound wave and MHD contributions are small and the GWs produced are dominantly from bubble collisions, The form of the GW energy density at the time of nucleation is given by [32], where we use v ∼ 1 for runaway bubbles. Quantities such as Ω * GW that are calculated at the time of nucleation are denoted with an asterisk, and they must then be evolved to relate to their values at the time of observation. S(f ) is the spectral shape function for the signal and a parametric from has been found through numerical simulations [64] of bubble wall collisions: The peak frequency, f p is a function of the temperature of the SM at the time of nucleation. The various hidden sectors can phase transition at different scales and therefore temperatures, causing a shift in the GW spectrum's peak frequency given by [64], where g * is calculated using Eq. (34), although, due to the lack of substantial reheating into the hidden sectors, the SM contribution is dominant. Now that the framework has been laid out for the creation of GW from a single SFOPT, we generalize to multiple sectors going under independent, coherent SFOPT. In the models presented in this paper, we consider a subset of N hidden sectors that undergo a phase transition at a SM temperature of T i γ . As the GWs propagate in free space, the energy density and frequency spectrum, at the time of production Ω * GW (f ), will redshift to today's value Ω 0 GW (f ) = A Ω * GW ((a 0 /a)f ). The redshifting factor, A, accounts for the redshifting of both ρ GW and ρ c [32,67], where a (a 0 ) and H (H 0 ) are the scale factor and Hubble constant at the time of nucleation (observation), respectively. Assuming that the sectors are completely decoupled before and after their respective SFOPT, the total GW signal that would be measured today is given by the coherent sum, We assume that the parameters of the SFOPT do not differ between sectors: the relativistic degrees of freedom, phase transition rate, and the dark QCD scale, are all similar. This makes the redshifting factor A i independent of sector number. Applying this to the standard reheating scenario of N natrualness, introduced in IV, we get GW signals as seen in Fig. 1. Plotted are the individual contributions to the signal from each phase transitioned sector, as well as the coherent sum of all sectors. Future GW interferometers and pulsar timing array sensitivity curves are shown in comparison to the signal. The sensitivity curves are interpreted as the region of possible detection if intersected with the GW signal, and the construction of these curves is detailed in Section VI B. Notice that the total signal is dominated by the first sector's contribution. This is caused by the quartic temperature ratio suppression in Eq. (33) and the large temperature gaps between adjacent sectors. Such a suppression leads to standard N natrualness evading future detector thresholds by a few orders of magnitude in units of energy density. This is not the case if we consider more generalized reheating scenarios. Once the restriction that sectors with small Higgs masses are preferentially reheated has been lifted, we can explore a much more vast landscape of hidden sectors than are allowed in the reheaton case. Here, we construct several different scenarios that are both detectable and demonstrate a variety of gravitational wave profiles. Specifically, we explore benchmarks that lead to a deviation in the peak behaviour of the total GW signal (the superposition of stochastic GW from individual SFOPT) from a standard power law signal.
It should be noted that the key phenomenological constraint on all of these models is ∆N ef f , giving us a maximum allowed temperature ratio (when compared to the SM) for each reheated hidden sector: Eq. (31) shows the maximum temperature ratios for specific numbers of additional hidden sectors. Due to the rather harsh scaling of the GW strength, α, with temperature ratio shown in Eq. (35), we take the optimistic approach of keeping the temperature ratio as high as allowed by CMB data for all of the hidden sectors.
In the following, we focus on heavy standard sectors -pure Yang-Mills sectors with much heavier particles (specifically quarks) and, as shown in Sec. III, the SFOPT these entail. The reason for this arises from Eq. (43): every exotic sector features a phase transition that occurs at Λ ex ∼ 90 MeV. If we maximize the allowed temperature ratio, this gives us a (SM) photon temperature, T γ , that places our signal directly in the frequency void between the detection region of pulsar timing arrays and space-based interferometers (see Sec. VI B). The location of the peak can be changed by dropping the temperature ratio, but the adjustment required to end up with a signal with an appropriate peak frequency makes the overall signal too weak to detect. As shown in Sec. III, standard sectors can have much higher temperature phase transitions. As such, maintaining the maximum allowed temperature ratio between the hidden sector(s) and the SM gives a much larger photon temperature and a proportionally larger peak frequency; ultimately allowing for detection by space-based interferometers.
There are four scenarios that we examine, with key parameters presented in Tab. I.
• Maximized signal: A single additional heavy hidden sector reheated to a temperature that saturates current experimental bounds. The SM photon bath temperature at the time of the hidden sector PT is 87 GeV. In the N naturalness framework this is equivalent to reheating a standard sector with i ∼ 10 16 up to the maximum allowed temperature ratio.
• Large split scenario: A scenario where two additional hidden sectors have been reheated -these sectors have Higgs vevs that are split by a factor of This results in a difference in the scale of the SFOPTs leading to the SM photon bath temperature changing a large amount during the time between the PTs. This, in turn, leads to a large separation in the peak frequency of their gravitational wave signals. In the N naturalness framework this is equivalent to reheating two standard sectors, one with i ∼ 10 12 and another with i ∼ 10 15 up to the maximum allowed temperature ratio.
• Medium split scenario: Similar to the previous case: these sectors have Higgs vevs that are split by a factor of resulting in a much smaller difference in the peak frequency of their gravitational wave signals. In the N naturalness framework this is equivalent to reheating two standard sectors, one with i ∼ 10 12 and another with i ∼ 10 13 up to the maximum allowed temperature ratio.
• Five sector scenario: Five sectors are reheated to the maximum allowed temperature ratio, each with vevs that are larger than the previous sector.
In all cases where multiple sectors are reheated, we assume for simplicity that all the hidden sectors are reheated to the same temperature. The GW results of these cases are presented in Fig. 2. In all cases, the summed GW signal is detectable by one or more proposed interferometers. When changing the assumptions on β/H, the scenarios in Fig. 2 are still detectable for values ranging between O(1) and O(100). As β/H increases (decreases) the peak frequency moves to higher (lower) frequencies, dictated by Eq. (43), where as the amplitude decreases (increases) shown in Eq. (41).
The frequency dependence in Eq. (42) takes the form of f /f p , this causes a cancellation between the redshifting factors. As multiple sectors phase transition at different times, and therefore different SM photon temperatures, the peaks will shift relative to each other, purely from the linear temperature dependence of the peak frequency f p ∼ T γ given in Eq. (43). This is seen in Fig. 2, where the spectrum peaks are shifted causing a peak broadening of the summed spectrum. The broadening can be substantial if the hidden sectors transition between a large   (4)). It should be noted that although the various sectors undergo phase transitions at different temperatures, they are all assumed to be reheated to the same initial temperature.
gap of time (temperature). Eventually, a temperature limit will be reached where two (or multiple) distinct peaks will be visible, provided that the amplitudes are comparable.

B. Detection of Stochastic Graviational Waves
A stochastic gravitational wave background could be detectable if the the signal-to-noise ratio (SNR) is above some threshold value, ρ > ρ th , dictated by the capabilities of future interferometers and pulsar timing arrarys (PTA). These interferometers/PTAs quote their experimental sensistivies in terms of spectral noise curves, S eff (f ), which can be translated into units of energy density through h 2 Ω eff (f ) = 2π 2 3H 2 f 3 S eff (f ). If the experiment uses a single (multiple) detector, the autocorrelated (cross-correlated) SNR is used in comparing to the threshold value ρ th . The auto-correlated and crosscorrelated SNR are explictely given as [69], (auto-correlated) where T is the exposure time of the experiment. The integration covers the entire broadband range of frequencies (f min , f max ). LISA [8] and B-DEICIGO [11] are proposed to be single detector interferometers, where as BBO [5] and DEICIGO [10] would be built from an array of multiple interferometers. GW signals produced from an early cosmological phase transition would be seen as a stocastic background. Assuming that the GW follows a power law background in frequency, it is commonplace to quote the power law integrated (PLI) sensistivity curves [68]. The PLI curves are constructed using information from the power law form of the signal, where γ is the spectral index of the power law, and f ref is an arbitrary reference frequency which has no effect on the PLI sensetivities. h 2 Ω γ is the energy density calculated using Eq. (49) with spectral index γ and reference frequency f ref . The method of calculating the PLI curves involves plotting h 2 Ω GW (f ), using Eq. (50), for various spectral indices γ and for some fixed threshold value of ρ th . Each curve will lay tangent to the PLI curve, more formally, FIG. 2. Gravitational wave spectral energy density for the various scenarios found in Tab. (I). All contributions are assumed to be purely from runaway bubble collisions Ω φ , with β/H = 10. The inset is a closer look at the region around the peaks. The shaded curves are the same as Fig. 1.
The spectral noise curves used to create the PLI curves shown in Figs. 1 and 2 were taken from [11,32,[70][71][72] for the interferometers, and [7,32] for the SKA pulsar timing array. We have assumed an observation time of T = 4 years for the interferometers and T = 5, 10, 20 years for the various stages of SKA. In the case of the PTA experiments, the sensitivity curves are dependent on how frequently the pulsar's timing residuals, δt, are measured. When using Eq. (49) to construct the PLI curves for SKA, the upper integration bound is inversely proportional to pulsar's timing residual, f min = 1/δt. In this work, it is assumed that δt = 14 days, but this may underestimate the capabilities of SKA as well as the cadences of the pulsar populations. If the timing residuals are lowered the maximum frequency reach of SKA increases, and the corresponding PLI curves in Figs. 1 & 2 are shifted to the right, possibly giving the PTAs sensitivity to some of the scenarios considered here.

VII. CONCLUSION
As detection capabilities increase, gravitational wave signals continue to grow in importance as phenomenological signatures that can offer us a unique glimpse into the universe as it was in the early epochs. The space-based interferometers planned for the next generation of GW experiments will be sensitive enough to begin searching for signals of the cataclysmic disruption of space-time due to SFOPT. As we inch closer to these measurements becoming available, it becomes important to develop ways to analyze and understand this data.
Here, we examined scenarios, including N naturalness, that involve multiple hidden sectors and calculated the GW profiles present. Our GW projections demonstrate that although N naturalness with the reheaton scenario presented in [38] is not projected to be detectable in the near future, more generalized scenarios with multiple hidden sector SFOPTs are in an observable region and will begin to be probed by next generation space experiments. Both cases feature important parts of their GW signals in the void between frequencies detectable by pulsar timing arrays and space-based interferometers -providing theoretical impetus for new experiments capable of probing this region of frequency space.
Further, our results provide a framework for understanding and using GW signals in two different ways: first as a unique signal for specific theories featuring multiple SFOPTs and also as a challenge to broaden the understanding of GW detector sensitivity.
In the former case, this demonstrates the power of GW signals to probe deep into the unknown arena of com-plex hidden sectors. Individual SFOPT are understood to create GW that are assumed to follow an approximate power law. If a model predicts of the presence of 2, 5, or more additional sectors, or features a single extra sector with multiple PTs, deviations from a standard power law can occur. The multiple transitions that occur in the models outlined here create signals that follow this trend: although the individual GW do obey approximate power laws, their sum does not -leading to a unique signal indicating so-called dark complexity. Explicitly, a broadening or distortion of the signal around the peak frequency, precisely where the signal has the most energy, could point to a multi-SFOPT scenario and gently guide us in the direction of multiple hidden sectors.
Shifting to the other part of our framework, our results leads to the question "how well can experiments probe non power law signals?" For frequency ranges away from the peak of the total, GW signals the quoted detection thresholds should hold: the signals fall off as a power law to a very good approximation. However, for areas around the peak frequency the answer is less clear; the PLI curves are built under the assumption of a power law. This points to the need for future work to better understand detection prospects for the next generation of GW detectors.