Gravitational Waves from Nnaturalness

We study the prospects for probing the Nnaturalness solution to the electroweak hierarchy problem with future gravitational wave observatories. Nnaturalness, in its simplest incarnation, predicts $N$ copies of the Standard Model with varying Higgs mass parameters. We show that in certain parameter regions the scalar reheaton transfers a substantial energy density to the sector with the smallest positive Higgs squared mass while remaining consistent with bounds on additional effective relativistic species. In this sector, all six quarks are much lighter than the corresponding QCD confinement scale, allowing for the possibility of a first-order chiral symmetry-breaking phase transition and an associated stochastic gravitational wave signal. We consider several scenarios characterizing the strongly-coupled phase transition dynamics and estimate the gravitational wave spectrum for each. Pulsar timing arrays (SKA), spaced-based interferometers (BBO, Ultimate-DECIGO, $\mu$Ares, asteroid ranging), and astrometric measurements (THEIA) all have the potential to explore new regions of Nnaturalness parameter space, complementing probes from next generation cosmic microwave background radiation experiments.


I. INTRODUCTION
The naturalness puzzle associated with the Higgs mass has for several decades inspired a vision of rich dynamics underlying the electroweak scale, possibly involving supersymmetry, new strong dynamics, or extra spatial dimensions, with a host of new states within reach of high energy colliders.However, the key lessons of the Large Hadron Collider (LHC), including the existence of a Higgs boson with properties in agreement with the Standard Model (SM) predictions and the absence thus far of new degrees of freedom at the TeV scale, have led physicists to question this traditional vision and pursue new lines of attack on the hierarchy problem, see for example Refs., which also include earlier relevant literature, as well as the recent reviews [23,24].Among these new ideas, one particularly interesting approach, known as Nnaturalness [25], posits N mutually non-interacting copies of the SM with Higgs mass parameters distributed over the range of the cutoff of the theory.In this way, some sectors will accidentally have Higgs mass parameters that are parametrically smaller than the cutoff, and our SM is identified with the sector having the smallest (negative) squared Higgs mass.Finally, a light 'reheaton' with universal portal couplings to each sector will naturally transfer most of its energy to our sector and only slight fractional energy densities to the other sectors, allowing for a viable cosmology with small, but potentially testable, departures from ΛCDM.
By construction, experimental and observational tests of Nnaturalness are scarce, with the most robust probes coming from cosmology [25].The extra energy deposited in the other sectors leads to dark radiation, which can be probed in current and future generation CMB experiments [26][27][28].Also, the slightly heavier neutrinos from other sectors may free stream around matter-radiation equality, which can suppress the matter power spectrum to a level that is potentially measurable [29,30].On the other hand, the possibility of probing the additional sectors or the reheaton at accelerator experiments is remote.It is therefore of great interest to find additional probes of the scenario.
In this work we investigate the prospects for probing Nnaturalness through gravitational wave (GW) signatures.The basic idea we will explore concerns the dynamics of QCD in the exotic sectors having positive squared Higgs masses.In such exotic sectors, all six quarks are light in comparison to the corresponding QCD confinement scale.Therefore, these exotic sectors may undergo a first-order phase transition (FOPT) associated with the breaking of the corresponding SU (6) L × SU (6) R × U (1) A chiral symmetry [31], which in turn generates a stochastic GW signal [32][33][34][35][36][37][38][39][40][41][42][43][44][45][46][47].This possibility was discussed previously in Ref. [48] which, however, concluded that for a particular point in parameter space the energy density in the exotic sectors was too small to lead to a detectable GW signal.We revisit this possibility and identify regions of the Nnaturalness parameter space where the exotic sector with the lightest Higgs contains enough energy density to yield a detectable GW signal, yet is still consistent with current bounds on additional relativistic degrees of freedom.Depending on the details of the QCD phase transition in this sector, which are unfortunately obscure due to strong dynamics, the corresponding GW wave signal is predicted to lie in nHz − Hz frequency range, with an amplitude that is potentially detectable by several current and planned GW observatories.
The discovery of GWs by the LIGO-Virgo-KAGRA collaborations [49][50][51][52] has opened a new observational window to the universe in the Hz − kHz frequency range.Existing and new observatories planned in the next decade and beyond will be capable of measuring GWs over a much wider frequency range and with significantly smaller amplitudes.While the GWs observed by LIGO-Virgo-KAGRA are sourced by mergers of compact objects such as O(1 − 100 M ⊙ ) black holes, new observatories offer the promise of probing a variety of GW signals and sources.This includes stochastic GWs, which may be generated by FOPTs, as well as other exotic sources [53][54][55][56][57][58][59][60][61][62][63][64][65][66][67][68][69][70][71].Recently, the NANOGrav pulsar timing array (PTA) reported evidence for a signal in the O(1 − 10) nHz frequency range with a Hellings-Downs correlation [72] that is characteristic of stochastic GWs in their 15-yr survey [73].This result is supported by the marginal evidence for a stochastic GW signal from the European PTA [74] and is also consistent with results from the Parkes PTA [75] and China PTA [76].This milestone, if eventually confirmed with future data, marks a new era of exploration of novel astrophysical and cosmological GW sources.
We explore the capability of several existing and future GW experiments, including PTAs, spaced-based interferometers, and astrometric measurements, to probe the Nnaturalness exotic sector QCD phase transition.To this end, we first identify the most promising regions of Nnaturalness parameter space for a GW signal where the first exotic sector contains a substantial energy density while remaining consistent with constraints from the CMB on new effective relativistic degrees of freedom.We consider several scenarios characterizing the dynamics of the exotic sector QCD phase transition, which allow us to explore a range of possible GW signals in a model-agnostic fashion.Depending on these assumptions, we find that experiments such as the Square Kilometer Array (SKA) PTA [77], the spacedbased interferometers LISA [78], BBO [79], Ultimate-DECIGO [80], and µAres [81], and asteroid laser ranging [82], and future astrometric measurements [83] by the proposed THEIA experiment [84] have the potential to probe Nnaturalness.On the other hand, we find that Nnaturalness is unlikely to account for all of the stochastic GW signal recently reported by NANOGrav due to stringent CMB constraints on new relativistic degrees of freedom.
The rest of this paper is organized as follows: In Section II we review the minimal Nnaturalness model with a scalar reheaton, focusing on the salient features, particularly of the first exotic sector, that will be used in our subsequent analysis of the cosmology and GW signal.Next, in Section III, we discuss the cosmology of the scenario and estimate the contribution from the other sectors to the extra effective relativistic degrees of freedom at late time.In Section IV we discuss our estimate of the GW signal under different assumptions regarding the nature of the exotic sector QCD phase transition.Our main results are presented in Section V, which include a delineation of the Nnaturalness parameter space that may potentially be probed by future GW observatories.Our conclusions and outlook are presented in Section VI.Appendices A and B contain technical details on the reheaton decays and example estimates of the effective relativistic degrees of freedom in the different sectors at several stages of the cosmological history, respectively.

II. NNATURALNESS
The minimal1 Nnaturalness model contains N copies of the SM that are mutually decoupled.The Higgs squared mass parameters are assumed to vary uniformly from one sector to another according to the relation Here i labels the sector, with our SM identified with the i = 0 sector such that m 2 The parameter r controls the relative distance of m2 H from zero, with 0 ≤ r < 2, and Λ H is the cutoff of the theory.The SM-like sectors with i > 0 have negative squared mass parameters, while the exotic sectors with i < 0 have positive squared mass parameters.
Besides the N sectors, the other crucial ingredient in Nnaturalness is the reheaton which is assumed to dominate the energy density of the universe at some time following inflation.
Ref. [25] considered models with a scalar reheaton and models with a fermionic reheaton.
For concreteness, in this work we focus on the real scalar ϕ reheaton, with Lagrangian where m ϕ is the reheaton mass and a is a universal dimensionful coupling of the reheaton to the Higgs fields.The couplings in Eq. ( 2) cause the reheaton to decay into all sectors.
As we will discuss in detail below, a light reheaton, with a mass near the electroweak scale, can dominantly decay to and populate the i = 0 SM sector, thereby allowing for a viable cosmology and a solution to the hierarchy problem.Besides the variation in their Higgs mass parameters, the sectors are assumed to be identical in all respects, which implies the theory has a softly broken sector permutation symmetry.The SM-like sectors, due to their large, negative Higgs squared masses, undergo electroweak symmetry breaking in the familiar way, ⟨H i ⟩ ̸ = 0, with the Higgs fields obtaining large vacuum expectation values (VEVs) , with λ the universal Higgs quartic coupling and v = 246 GeV the SM Higgs VEV.Instead, in the exotic sectors the Higgs squared masses are positive, and electroweak symmetry breaking is triggered by QCD strong dynamics through the formation of a quark condensate ⟨qq⟩ i ̸ = 0.The exotic sector quarks receive masses of order m q i ∼ y q y t Λ 3 QCD i /m 2 H i and therefore are all much lighter than corresponding confinement scale Λ QCD i ∼ O(100 MeV).Hence, these exotic sectors, with six light quark flavors, may undergo FOPTs associated with the breaking of the corresponding SU (6) L × SU (6) R × U (1) A chiral symmetry [31], which then generates a stochastic GW signal.We note that there is still some debate in the literature on the order of this phase transition, and we comment further on this in Sec.III.Assuming the phase transition is first order, the detectability of the GW signal depends on how much energy density is contained in the exotic sectors, and only the first exotic sector (i = −1) may have a substantial energy density in the cosmologically allowed regions of parameter space.To understand this, we must carefully examine the cosmological evolution of the model, which, in any case, is of central importance in the Nnaturaless solution to the hierarchy problem.
This will be discussed in detail in the Section III.
The minimal Nnaturalness model considered here is thus characterized by four parameters: N , r, m ϕ , and a.The GW signal will be relatively insensitive to the value of N because the signal originates solely from the i = −1 exotic sector.For concreteness, for the rest of this work we fix N = 10 4 which allows for a solution to the little hierarchy problem, with Λ H ∼ 10 TeV, and evades any potential issues with overclosure from massive stable states from the other sectors.Likewise, the Nnaturalness mechanism and the GW signal are not sensitive to the precise value of the universal coupling a since it cancels out in the reheaton decay branching ratios.The only requirement is that a is small enough so that the reheating temperature is below the electroweak scale, which can always be satisfied.The GW signals will therefore be controlled by m ϕ and r.

A. Reheaton decays
The fraction of the reheaton energy density transferred to each sector is proportional to the reheaton partial decay width into each sector, which we denote by Γ i .If the reheaton is light, with mass of order the electroweak scale, it can dominantly decay to the SM sector.
Depending on the values of m ϕ and r, there may also be a significant energy density stored in the other sectors, leading to constraints and probes from additional relativistic degrees of freedom and GWs, to be discussed in the next sections.Some details related to the reheaton decays to the SM-like and exotic sectors are provided in Appendix A; we now summarize their basic properties.
In the SM and SM-like sectors electroweak symmetry breaking causes ϕ to mix with the corresponding physical Higgs boson h i , with mixing angle θ i ≃ a v i /m 2 h i ≈ a/m h i , where m h i is the physical Higgs mass for the sector.Thus, the ϕ partial decay widths to SM-like sectors scale as Γ i ∝ 1/m 2 h i , with decays to the SM being the largest.Eq. (1) implies that m h i decreases as r is increased.Thus, the fractional energy densities deposited in the i ≥ 1 SM-like sectors tend to increase as r increases.
In the exotic sectors, the small effects of electroweak symmetry breaking from QCD can be neglected as far as the decays of the reheaton are concerned.Except for perhaps the lightest exotic sectors, we expect m ϕ < m H i , in which case the reheaton decays to the exotic sector i mainly proceed through a loop via ϕ → W i W i , B i B i , with a decay width given by Γ ∝ 1/m 4 H i , or through a four-body decay ϕ → H * i H * i .Therefore, the energy density stored in the heavier exotic sectors is generally insignificant in the viable regions of parameter space.Only the first exotic sector (i = −1) may potentially receive a substantial portion of the reheaton's energy density.As can be seen from Eq. ( 1), as r is increased the first exotic sector Higgs mass m H −1 decreases, and for m H −1 ∼ m ϕ /2 or below the reheaton may have a sizable or even dominant branching ratio into the lightest exotic sector via the two-or three-body decay Besides the general trends outlined above, when the reheaton mass is close to the SM Higgs mass, m ϕ ∼ m h , there is a resonant enhancement in ϕ−h mixing, θ SM ≃ av/(m 2 h −m 2 ϕ ), which enhances the decays of the reheaton to the SM sector.Another such enhancement occurs for m ϕ ≳ 2m W , 2m Z , when reheaton decays to on-shell SM weak bosons open up.

B. Properties of the First Exotic Sector
As explained above, among the exotic sectors only the first one (i = −1) may have a significant fraction of the reheaton energy density in the cosmologically viable regions of parameter space.Thus, it is only this sector which may furnish a potentially detectable stochastic GW signal from its corresponding QCD FOPT.
In this sector, the Higgs squared mass is positive.Hence, in the absence of QCD strong dynamics, electroweak symmetry would not be spontaneously broken, and all fermions and gauge bosons would be massless.However, as in the SM sector, QCD in this sector becomes strongly interacting at scales of order 1 GeV, and a quark condensate forms, ⟨qq⟩ −1 ∼ 4πf 3 This in turn generates masses for the leptons and quarks (though the latter are confined into hadrons at low energy).The chiral symmetry is explicitly broken by the Yukawa couplings and the electroweak gauge interactions, and this explicit breaking will cause the remaining 32 pions to obtain masses, i.e., they are pseudo-Nambu Goldstone bosons (pNGBs).
We now provide some results for the mass spectrum of the exotic sector states lighter than confinement scale, which will be relevant in our discussion of the cosmology and GW signal.For our quoted numerical estimates in the following, we choose a benchmark m H −1 = 70 GeV and take f π −1 = 30 MeV. 3 The electroweak gauge boson masses are given by m ing numerical estimates of order 20 MeV.The 32 pions receive masses of parametric size from the explicit chiral symmetry breaking by the quark Yukawa couplings.We find the pions range in mass between about 1 keV and 100 keV.The leptons ) yielding estimates of 0.1 meV, 20 meV, 0.3 eV for the electron, muon, and tau, respectively.Neutrinos are expected to be extremely light, m ν −1 < 10 −11 eV, while the photon is massless.

III. COSMOLOGICAL EVOLUTION
The cosmological history of the Nnaturalness model starts when the reheaton dominates the energy density of the universe.The reheaton then decays into all available channels, reheating the universe such that each sector is populated with an energy density that scales with the reheaton's partial decay width in that sector ρ i /ρ SM ≃ Γ i /Γ SM .Each sector thermalizes within its own sector with corresponding energy and entropy densities given by where T denotes the SM temperature, ξ i ≡ T i /T is the ratio of ith sector temperature to that of the SM, and g * ρ,i (g * s,i ) denotes the the effective number of relativistic (entropy) degrees of freedom in sector i.We refer the reader to Appendix B for example estimates of the effective relativistic degrees of freedom at various cosmological epochs.A lower energy density equates to a colder temperature (relative to the temperature of the SM bath) because Thus sectors with larger |i| will be increasingly cold. 4  3 We assume that f πi scales linearly with Λ QCD i .The value of Λ QCD i changes weakly with i due to the change in the quark mass thresholds.In particular, taking α S (m Z ) SM ≃ 0.118 as input and using one-loop running, we find that Λ QCD −1 /Λ QCD SM ≈ 0.3. 4We assume that the baryon asymmetry in all SM-like and all exotic sectors is negligible otherwise the additional matter would overclose the universe [25].
The reheat temperature T RH of the SM bath can be taken as a free parameter since it is governed by the free coupling a, and we consider Λ QCD −1 ≲ T RH ≲ v.The upper bound must be imposed to avoid finite temperature corrections to the Higgs potential which would spoil the Nnaturalness mechanism, while we impose the lower bound so that the first exotic sector is reheated above its corresponding confinement scale such that the sector experiences a cosmological FOPT.We will fix T RH = 100 GeV for concreteness.The temperature ratio for sector i at reheating is given by The SM sector then proceeds following the usual cosmological evolution, while each of the SM-like sectors evolves in a similar way to the SM sector.In particular, the ordering of neutrino decoupling, electron-positron annihilation, and photon recombination in these sectors is the same as in the SM [85].The radiation in each of these sectors, in the form of free-streaming photons and neutrinos, will contribute to ∆N eff which measures additional relativistic degrees of freedom relative to the SM.As discussed in the previous section, of the exotic sectors only the first one may have a substantial energy density, and the cosmological evolution of this sector features several qualitative differences from that of the SM-like sectors.We will outline these differences in detail below, but we point out here that because the spectra of particles in the first exotic sector are much lighter than in the SM sector, the photons can be interacting until much later times.These sectors also contribute to ∆N eff but behave as interacting radiation rather than free streaming.The dominant cosmological signal is thus an unavoidable contribution to ∆N eff whose size is determined by the relevant partial width into the SM sector, compared to the sum of all other sectors.5

A. Cosmology of the First Exotic Sector
Once populated by the decay of the reheaton, the first exotic sector thermalizes and cools as the universe expands.During this initial period of evolution, all degrees of freedom are essentially massless except for the Higgs doublet H −1 .
As the exotic sector cools to temperatures of order Λ QCD −1 ∼ O(100) MeV, a chiral symmetry breaking phase transition is precipitated by the formation of the quark condensate.
The conventional wisdom, due to an argument of Pisarski and Wilczek [31], is that this phase transition is first order.Employing a linear sigma model description of the quark bilinear order parameter, they performed a renormalization group analysis using a perturbative ϵ expansion and noted the absence of infrared stable fixed points for N f ≥ 3 light quark flavors, which is indicative of a first-order phase transition.Subsequent studies using phenomenological models have confirmed this result, see, e.g., Ref. [87].This question has also been studied at various points on the lattice over the past decades, with some confirming the claim of a first-order transition for N f ≥ 3 [88,89] and others challenging it [90].
Thus, the question of the order of the phase transition is still an open one and further study is needed to settle the issue, see Ref. [91] for some perspectives in this direction.We will follow the conventional wisdom and assume that the exotic sector phase transition with six light flavors is first order.
The phase transition commences at the critical temperature T crit −1 , the point at which the potential energy of the true and false vacua are the same.For the exotic sector QCD with six massless quarks, we take T crit −1 = 85 MeV, which is estimated with order 20% uncertainty [92].Starting from the exotic sector in the symmetric phase, bubbles of true vacuum nucleate, expand, and merge, such that eventually the sector ends up in the broken phase.The nucleation temperature, T nuc −1 ≲ T crit −1 , marks the point at which the first bubbles nucleate.After nucleation, the bubbles take time to percolate until 34% end up in the true vacuum corresponding to temperature T perc −1 [93].Thus it is reasonable to assume T perc −1 ≲ 85 MeV, and we consider the range 50 MeV ≤ T perc −1 ≤ 85 MeV.The temperature of the SM at percolation is T perc = T perc −1 /ξ perc −1 where ξ perc −1 is the temperature ratio between the exotic sector and SM right before percolation, i.e., in the unbroken phase.
The strength of the exotic sector phase transition is characterized by the the parameters α −1 and α tot , which are defined as where ∆θ −1 is the difference in the trace of the energy momentum tensor in the unbroken and the broken phases.Due to the strongly-coupled dynamics during the exotic sector QCD phase transition, we will not be able to provide a first principles calculation of α −1 .Instead, we consider several distinct scenarios for the phase transition dynamics with varying choices for α −1 .We defer a detailed discussion of the considerations underlying these assumptions to Sec.IV.
Once the phase transition concludes, the exotic sector is reheated to a temperature T rh −1 .6Assuming an instantaneous transition from percolation to reheating and using energy conservation, we may write where ∆V −1 is the difference in free energy.Using Eqs.(5,7) and assuming ∆V −1 ≃ ∆θ −1 (see e.g., Ref. [94]), we may write The QCD FOPT results in entropy production in the exotic sector, which may be encoded in the ratio of entropy densities before and after the QCD phase transition, where we have used Eq. ( 8).The temperature of the SM at the end of the phase transition is T rh = T rh −1 /ξ rh −1 where ξ rh −1 is the temperature ratio between the exotic sector and SM right at the end of the phase transition, i.e., in the broken phase.Assuming instantaneous reheating we have T rh = T perc .
Following the QCD phase transition, entropy is conserved in the exotic sector throughout its subsequent evolution.As discussed in the previous section, the light degrees of freedom with masses below Λ QCD −1 consist of the electroweak gauge bosons, pions, charged leptons, neutrinos, and photons.As the temperature drops below their masses, the electroweak gauge bosons and pions leave the exotic sector bath.Interestingly, neutrinos in this sector typically decouple while both the muon and tau are relativistic.To see this, we estimate the neutrino scattering rate as Γ ν,−1 ∼ G 2 F −1 (ξ −1 T ) 5 and compare it to the Hubble rate.Noting that G F −1 ∼ f −2 π −1 in the exotic sector, this gives the decoupling temperature as Given that the charged lepton masses discussed in the previous section are typically below the eV scale, we see from Eq. ( 10) that neutrino decoupling in the exotic sector typically happens before electrons, muons, and taus annihilate.Following neutrino decoupling, and somewhat before recombination, the taus annihilate and heat the photon bath relative to the neutrinos by a factor which can be derived in the standard way using entropy conservation arguments.Thus, near recombination the exotic sector relativistic species comprise photons, neutrinos, electrons, and muons.The muons eventually annihilate at late times while the electrons and photons remain in equilibrium until today.

B. ∆N eff in Nnaturalness
The most important constraint on Nnaturalness comes from bounds on ∆N eff during the epoch of recombination.Bounds from Planck, namely Planck + Lensing + BAO [95], constrain free streaming ∆N CMB eff ≤ 0.3 (all bounds quoted here are at the 95% confidence level).The exotic sector more closely corresponds to an interacting fluid which results in a slightly weaker bound of ∆N CMB eff ≤ 0.45 [96].There is a well-known tension between data from Planck and data from SH0ES, but when the data from SH0ES [97] is incorporated the bound on interacting radiation is further relaxed to ∆N CMB eff ≤ 0.7 [98].We use ∆N CMB We can relate ξ CMB i to ξ RH i in Eq. ( 4), which depends on the reheaton partial decay width ratio Γ i /Γ SM and thus on the Nnaturalness model parameters m ϕ and r.
We first consider the contribution to Eq. ( 12) from the SM-like sectors.Using the fact that the total entropy in both the SM sector and the SM-like sectors is conserved between the epochs of reheating and recombination, along with Eq. ( 4), we obtain For the exotic sectors, only the first such sector may potentially give a significant contribution to ∆N eff , so we focus our discussion on that sector.In comparison to the SM-like sectors, one important difference in this sector is the entropy production due to the QCD FOPT.To account for this, we first relate ξ CMB −1 to ξ rh −1 using the fact that entropy is conserved between the end of the exotic sector QCD phase transition and CMB epoch.Next, we account for the change in the exotic sector temperature during the phase transition, from the time of percolation to reheating, given by Eq. ( 8).This equation gives a relation between ξ rh −1 and ξ perc −1 .Finally, we may again use entropy conservation to relate the ξ perc −1 to ξ RH −1 .The final result for the first exotic sector contribution ∆N CMB eff,−1 is In comparison to the contributions from SM-like sectors in Eq. ( 13), one key difference in Eq. ( 14) is the presence of the factor D s,−1 given in Eq. ( 9), which encodes the entropy production in the exotic sector due to the FOPT.Using the benchmark values for the relativistic degrees of freedom given in Appendix B, Eq. ( 14) gives the relation Γ −1 /Γ SM ≈ 0.1(1 + α −1 ) −1 ∆N CMB eff,−1 /0.7 .Using Eqs.(13,14), in Fig. 1  The dashed contours show ∆N CMB eff,−1 and demonstrate that there is even viable space where the energy density in the i = −1 exotic sector is larger than the sum over all of the SM-like sectors.The primary reason this is possible is illustrated by the light blue shaded region which shows where the two-body decay ϕ → H −1 H −1 goes on-shell increasing its branching ratio substantially.It is this region where ∆N eff is as large as possible, while not violating existing constraints, and is dominated by the i = −1 exotic sector that a GW signal may be observable.

IV. GRAVITATIONAL WAVE SIGNAL
The first exotic sector with positive Higgs squared mass is predicted to have a substantial energy density, consistent with bounds on ∆N eff , in certain regions of the Nnaturalness parameter space.This sector may experience a cosmological QCD FOPT, and we investigate its associated stochastic GW signal.
A cosmological FOPT can produce a stochastic GW signal due to several effects.During the phase transition, GWs are produced through collisions of bubble walls [35,36,40].Additionally, production of GWs following the phase transition occurs due to sounds waves [43,44,46] and magnetohydrodynamic turbulence [38,39,42] in the plasma.The relevant physical quantity characterizing the GW signal is the differential GW density parameter, Ω GW (f ) = (1/ρ c ) dρ GW /d log f , where f is the frequency of the GW and ρ c is the critical density.Sophisticated numerical simulations have been performed to properly model these dynamical processes and predict the resulting GW spectrum.The GW spectrum at emission can then be conveniently described by a semi-analytical parameterization that is fit to the results of numerical simulations.Following Ref. [99], which makes use of the results from Refs.[41,42,44], we employ the following parameterization: where I = BW, SW denotes bubble walls and sound waves, respectively, and f em is the GW frequency at emission.We do not consider the contribution from turbulence in this work since this source suffers from significant uncertainties [100][101][102][103][104].
We see from Eq. ( 15) that the GW spectrum depends on the phase transition strength parameters α −1 and α tot , defined in Eq. ( 5), the phase transition duration parameter β/H, the bubble wall velocity v w , and the efficiency factors κ BW (fraction of the vacuum energy carried by the bubble walls during collision) and κ SW (energy fraction transferred to plasma bulk motion).We will return shortly to discuss our assumptions for these quantities, which depend on the detailed nature of the exotic sector QCD phase transition.For the other quantities appearing in Eq. ( 15) we use the results from Refs.[41,105].For the normalization factors, we have (N BW , N SW ) = (1, 0.159).The velocity factor takes into account a potential suppression due to the wall velocity, with (∆ BW , ∆ SW ) = ((0.11v 3 w )/(0.42 + v 3 w ), 1).The exponents are given by (p BW , p SW ) = (2, 2) and (q BW , q SW ) = (2, 1).The spectral shape functions and corresponding peak frequencies are taken to be For the sound wave contribution, we also include a additional suppression factor [106,107] for large β/H given by where c s = 1/ √ 3 is the speed of sound in the relativistic plasma.
In obtaining the observed spectrum today, one must account for the expansion of the universe from the time of GW emission until today, which redshifts both the energy density and the GW frequency: Here Ω 0 GW (Ω em GW ) denotes the spectrum today (at emission), f is the frequency today, a 0 (a perc ) is the scale factor today (at percolation), and R is a redshift factor.We define the time of emission to coincide with the time of percolation, when a substantial fraction of the universe is filled with bubbles of the true vacuum.Neglecting the small effect of entropy production during the exotic sector QCD phase transition, the relevant factors in Eq. ( 18) are given by a 0 a perc = g perc * s,tot g 0 * s,tot where T 0 = 2.725 K ≈ 0.235 meV is the present temperature of the CMB, H 0 (H perc ) is the Hubble parameter today (at percolation) with H 0 = 100 h km Mpc −1 s −1 and h 2 Ω 0 γ ≈ 2.47×10 −5 is the present photon density parameter.The factors counting relativistic degrees of freedom are given by We now return to discuss our assumptions regarding nature of the exotic sector QCD phase transition as well as the key parameters governing the GW spectrum, namely, α −1 , β/H, v w , and the efficiency factors.In principle, if the temperature-dependent effective potential describing the phase transition is known, one can compute these quantities.However, in our scenario, due to the associated strong dynamics, we are not able to provide a first principles analysis of the phase transition properties.Ideally, the phase transition could be studied using lattice methods (see Ref. [91] for some perspectives), though there are no existing studies which map on to our scenario.Attempts have been made in literature to model the effective potential and resulting GW signal in certain strongly-coupled QCD-like gauge theories using various phenomenological approaches/toy models, including linear sigma models, Polyakov-Nambu-Jona-Lasino models, and holographic models, see Refs. [108][109][110][111][112][113] for some recent representative studies.In many cases, these studies indicate relatively small (large) values of the phase transition strength (duration) parameters.For the Nnaturalness model, while it is not guaranteed, such values may still be potentially detectable by future space-based GW observatories, as we will discuss in Sec.V. We will not attempt to model the effective potential in this work, but will instead remain agnostic about the evolution of the phase transition.To illustrate the range of possibilities, we will consider several representative benchmark scenarios for the behavior of the phase transition and the parameters governing the spectrum, as we explain in the following.
Once a bubble nucleates the bubble wall experiences negative pressure from the potential difference between the true and false vacua causing it to accelerate.At the same time, the bubble wall faces pressure from the plasma in the symmetric phase which acts as friction on the expanding bubble wall.The largest frictional pressure a bubble wall faces is when the wall velocity v w approaches the Jouguet velocity v J [114].If the negative pressure from the potential difference between the true and false vacua (∆V ) is large enough (corresponding to large α −1 ) to overcome this maximum pressure from the plasma, the walls will exhibit ultra-relativistic velocities corresponding to a runaway scenario.On the other hand if the ∆V is not large enough (corresponding to a smaller α −1 ) to overcome the maximum frictional pressure from plasma, the wall reaches a terminal velocity with v w ∼ c s corresponding to a non-runaway scenario.Following the analysis in [115], we have checked that the boundary on the strength parameter between the runaway and non-runaway scenarios is given by α −1 ≈ 0.3 with larger α −1 corresponding to a runaway wall and smaller α −1 corresponding to a non-runaway wall with a terminal wall velocity that can be approximated by the speed of sound in the plasma, Another important parameter governing the GW spectrum is β/H, which is defined as where S 3 is the three-dimensional Euclidean bounce action, assuming the phase transition proceeds due to thermal fluctuations.The parameter β gives a measure of the duration of the phase transition.As is clear from Eq. ( 21), the calculation of β/H requires knowledge the tunneling action S 3 and hence the thermal potential during the phase transition, which, as mentioned above, is obscure due to the QCD strong dynamics.We will thus consider several benchmark choices for the phase transition duration parameter in the broad range 4 ].The lower bound is imposed to ensure efficient bubble percolation [116].
We note that the duration parameter β/H is expected to be inversely correlated with the strength parameter α −1 , see, e.g., Ref. [106] for discussion.
For our runaway scenario the energy of the phase transition is converted into accelerating the bubble wall implying that the dominant source of GWs is bubble collisions.For our nonrunaway scenario with a terminal wall velocity, the expanding wall pushes on the plasma in the symmetric phase, creating a coherent motion of the plasma.Therefore, most of the latent heat released during the phase transition is converted to sound waves.Given the above considerations, we will study the following two scenarios, which are illustrative of the range possibilities for the properties of the phase transition: • Runaway scenario: (5,10), (1, 10 3 ).
For the non-runaway scenario with v w = 1/ √ 3, we have used the numerical fitting function for the efficiency factor κ SW from Ref. [117] The amplitude of the GW signal is governed by the strength parameter α tot in Eq. ( 5), which is given by the product of α −1 and the fraction of the total energy density stored in the i = −1 sector.It is useful to ask how large this parameter may be while maintaining consistency with the ∆N eff constraints discussed in the previous section.Focusing on the regions of parameter space in which the first exotic sector provides the dominant contribution to additional relativistic degrees of freedom, ∆N eff ≈ ∆N eff,−1 , we may then relate α tot to ∆N eff,−1 as follows: The first line follows from Eq. ( 5) given that ρ perc tot ≈ ρ perc SM for parameters consistent with ∆N eff bounds.In the second line we have first related ξ perc −1 to ξ rh −1 using Eq. ( 8), and then related ξ rh −1 to ∆N CMB eff,−1 using entropy conservation and Eq. ( 12), assuming ∆N eff ≈ ∆N eff,−1 .Eq. (24) demonstrates that the strength parameter α tot may potentially be large enough to enable a detectable GW signal while satisfying bounds on additional relativistic degrees of freedom.

V. RESULTS AND DISCUSSION
Using the results of the previous section, in Fig. 2 we show the GW spectrum for the runaway and non-runaway scenarios defined in Eqs. ( 22) and ( 23), respectively.To exhibit FIG. 2. The GW spectrum from the runaway (blue) and non-runaway (orange) scenarios defined in Eqs.(22) and (23), respectively.To assess how large the signal can be, we assume the exotic sector contributes ∆N CMB eff,−1 = 0.7, saturating the bound from Planck + Lensing + BAO + SH0ES.Also shown are the NANOGrav 15-yr results (teal violin), the PLIS curves for upcoming experiments SKA and LISA (solid gray) and proposed experiments THEIA, µAres, asteroid laser ranging, and BBO (dashed gray).The PLIS curves for SKA, LISA, and BBO are adopted from [118] but scaled to observation times of 20 yrs [77] for SKA, 3 yrs for LISA [119], and 4 yrs for BBO [79].
The µAres PLIS is taken from [81], scaled to SNR = 1.For the asteroid ranging proposal, we adopt the strain sensitivity given in [82] and calculate the PLIS curve using the procedure outlined in [119] for SNR = 1 and assumed experiment duration of 7 yrs.For THEIA we adopt the PLIS sensitivity calculated in [120] for SNR = 1 and a mission lifetime of 20 yrs.For Ultimate-DECIGO (UDECIGO) we have adopted the PLIS in [80].Black dashed lines represent foregrounds from galactic and extragalactic compact binaries (CB) [121,122] and the SMBHB best fit to the NANOGrav 15-yr measurement [73].
Assuming that the astrophysical foregrounds either can be resolved and subtracted (see Ref. [125] for SMBHB foreground resolution) or are somewhat weaker in strength than currently expected, Fig. 2 demonstrates that there are promising opportunities to probe Nnaturalness with future GW measurements.As emphasized several times, this depends sensitively on the precise nature of the first exotic sector QCD phase transition, about which there are significant theoretical uncertainties, as well as the fractional energy density contained in this sector, which is dictated by the Nnaturalness model parameters, as we will discuss in detail shortly.Runaway transitions can be probed by PTAs, astrometric measurements, and spaced-based interferometers, while non-runaway transitions could lead to a signal in space-based interferometers.It is also clear from Fig. 2 that even under the most optimistic assumptions (runaway transitions, large α −1 , small β/H, and maximal ∆N CMB eff,−1 ) it is unlikely that Nnaturalness can fully account for the stochastic gravitational background reported in the NANOGrav 15-yr dataset (a similar point was made recently for the NANOGrav 12.5-yr dataset for generic stable secluded sectors [94]).Furthermore, we observe that if α −1 is too small, or β/H is too large, as may be suggested by detailed studies of FOPTs in toy models of QCD-like theories (see discussion in previous section), the GW signal from Nnaturalness may lie outside the reach of proposed experiments.The different scenarios considered here, Eqs.(22,23), serve to illustrate the range of possibilities.
Next, we map out the regions of the Nnaturalness parameter space that can potentially be probed by future GW experiments.Specifically, we determine the values of m ϕ and r that yield a GW signal intersecting (or tangential to) the PLIS curves shown in Fig. 2.
In Fig. 3 we show this reach for two runaway scenarios, (α −1 , β/H) = (5, 10) (top) and (α −1 , β/H) = (1, 10 3 ) (bottom).For each scenario, we show both the full parameter space for a reheaton mass lighter than 300 GeV (left), as well as a zoomed-in region of parameter space near m ϕ ≈ m h (right). 8The figures also show the predictions for ∆N CMB The behavior of the GW sensitivity curves in Figs. 3 and 4 can be understood by recalling that the strength parameter α tot is approximately linearly related to ∆N eff,−1 in the parameter regions where the total energy density is dominated by the SM bath, see Eq. ( 24).Thus, the GW sensitivity curves largely overlap with isocontours of ∆N eff,−1 , as can be seen by comparing with Fig. 1.The reach of GW experiments is strongest in the regions of parameter space where the reheaton has a relatively sizable branching ratio into the first exotic sector.
In the allowed regions of parameter space, consistent with CMB constraints on ∆N eff , this may occur in the regions where m ϕ ≲ 2m H −1 and the reheaton decays dominantly to the SM.These two requirements combine to sculpt two regions where GW experiments can be sensitive: 1) near the "Higgs funnel", m ϕ ≈ m h , and 2) above the threshold for reheaton decays to SM gauge bosons, as is observed in Figs. 3 and 4.

VI. CONCLUSIONS
Nnaturalness is a novel approach to the hierarchy problem.The key prediction of the framework is the existence of many decoupled hidden sectors containing small fractional energy densities, which can be probed through cosmological measurements such as ∆N eff .
In this work, have explored the potential to probe Nnaturalness through GW observations.
Considering the scalar reheaton model for concreteness, in certain parameter regions the first exotic sector, with the smallest positive squared Higgs mass, is predicted to have a sizable fractional energy density.QCD in this sector is expected to feature a cosmological first-order chiral symmetry breaking phase transition since all quarks are much lighter than the confinement scale, which then yields an associated stochastic GW signal.The resulting GW spectra are expected to peak in the nHz − mHz frequency range, with a strength scaling with the fraction of the reheaton energy density stored in the first exotic sector.
We have delineated the regions of parameter space where the first exotic sector has a substantial energy density, and thus a potentially detectable stochastic GW signature.The observational prospects for the GW signal depends sensitively on the detailed evolution of the exotic sector QCD phase transition, which involves strong coupling dynamics.We have remained agnostic about the phase transition properties, exploring several scenarios designed to encompass the spectrum of conceivable possibilities.Depending on these assumptions, as well as the eventual capabilities to discriminate various astrophysical foregrounds, we find that a GW signal from Nnaturalness can potentially be observable in several future exper- Our study reveals several interesting open questions.First, it would be valuable to further clarify the nature of the exotic sector QCD phase transition, using lattice studies as well as phenomenological models.Additional uncertainties in our predictions come from the modeling of the GW production from a first-order phase transition (see, e.g., Ref. [127] for a recent discussion).Progress on these issues will lead to a better understanding of the capabilities of GW observatories to probe Nnaturalness.
We have only considered the scalar reheaton model in this work, and it would be very interesting to also explore the potential gravitational wave signatures of the other reheaton models considered in Ref. [25].More accurate studies of the cosmological perturbations and their impact on the CMB and structure formation in the regions where the exotic sectors are populated would also be valuable (for a study considering the SM-like sectors, see Ref [28]).The exotic sector QCD FOPT phase transition could also be associated with other novel phenomena, such as the formation of dark quark nuggets [32,108] in the presence of a corresponding baryon asymmetry, or the production of primordial black holes [128][129][130].Future studies along these directions may point the way to even more new probes of Nnaturalness.Estimates are given at the epochs of reheating from reheaton decay (RH), exotic sector QCD FOPT at the point of percolation (perc) and reheating (rh), SM sector recombination (CMB), and today (0).We have assumed T RH = 100 GeV.
QCD phase transition, we estimate g perc

π − 1 with f π − 1
the corresponding pion decay constant, spontaneously breaking the approximate global chiral symmetry SU (6) L × SU (6) R → SU (6) V .This condensate thus also breaks the weakly gauged electroweak subgroup down to electromagnetism in the usual way.Of the 35 pions associated with this chiral symmetry breaking, three linear combinations form the true Nambu-Goldstone bosons eaten by the W ± and Z bosons to give them masses.The quark condensate triggers an effective tadpole for the Higgs field, inducing a VEV ⟨H −1 ⟩ ̸ = 0.
the default constraint on ∆N CMB eff .Comparable bounds can be placed on ∆N eff during the epoch of big bang nucleosynthesis (BBN), however, Nnaturalness generically predicts ∆N CMB eff > ∆N BBN eff .We evaluate ∆N CMB eff in the Nnaturalness model near the epoch of recombination at a SM temperature T CMB = 0.3 eV.Including the contributions from all sectors, this is we show several contours of ∆N CMB eff in the parameter space of m ϕ and r.The solid contours show the total ∆N CMB eff from all sectors.We see that for r ≲ 0.2 any mass of ϕ passes constraints from ∆N CMB eff .For larger values of r there are two viable regions.The first region is where 110 GeV ≲ m ϕ ≲ 140 GeV.Here the mixing between the reheaton and the SM Higgs grows much larger than the mixings between the reheaton and the Higgs particles from the other sectors.The large relative energy density in the SM means ∆N eff is small.Values of r ≳ 1 are possible in this region.The second region, which permits r ≳ 0.5 is where 160 GeV ≲ m ϕ ≲ 230 GeV.In this region the decay of the reheaton to a pair of SM W bosons goes on-shell which increases the relative energy density in the SM.As the mass of the reheaton increases the energy density in the SM-like sectors and exotic sectors grows which leads to the upper limit of this region.

�N = 10 4 TFIG. 1 .
FIG. 1. Contours of ∆N eff in the r vs. m ϕ space, evaluated at the time of the CMB.The number of sectors is N = 10 4 , the reheating temperature is T RH = 100 GeV, and the phase transition strength is α −1 = 1.The solid contours show the total ∆N eff contribution while the dashed contours show the contribution only from the first exotic sector.The light blue shaded region shows where the decay ϕ → H −1 H −1 is on-shell and the dark blue shaded region shows where the decay ϕ → H −2 H −2 is on-shell.

7
∆N CMB eff,−1 = 0.7 is only used for Fig.2.When results are shown in the Nnaturalness parameter space, as in Figs.3 and 4the ∆N eff limits are compared to ∆N CMB eff,tot .
the Planck + Lensing + BAO + SH0ES data.In the most optimistic scenario (Fig.3, top), there are several experiments and techniques (PTAs, space-based interferometry, astrometry) that can explore uncharted Nnaturalness parameter space, and, in particular, µAres, Ultimate-DECIGO, and THEIA even have the potential to compete in reach with future precision CMB measurements, e.g., CMB Stage IV[126] (∆N CMB eff ≲ 0.03).Similarly, in Fig.4we show the reach of future GW experiments for two non-runaway scenarios, (α −1 , β/H) = (0.3, 10 2 ) (top) and (α −1 , β/H) = (0.1, 10 3 ) (bottom).For these scenarios, the GW signal from sound waves is predicted to lie in the µHz range and thus can potentially be probed by future space-based interferometers such as µAres and asteroid ranging.It is worth noting that the phase transition parameters for the second scenario, (α −1 , β/H) = (0.1, 10 3 ) (bottom), are broadly consistent with results from studies modeling the phase transitions of QCD-like theories.

�� 3 FIG. 3 .
FIG. 3. Reach from various gravitational wave experiments in the r vs. m ϕ space over the full parameter region (left) and zoomed into the Higgs funnel (right).Results are shown for the runaway scenario with β/H = 10 and α −1 = 5 (top) and β/H = 10 3 and α −1 = 1 (bottom).Short dashed lines indicate contours of ∆N eff evaluated at T CMB .
iments, including PTAs (SKA), planned (BBO) and proposed (Ultimate-DECIGO, µAres, asteroid ranging) spaced-based interferometers, and astrometric measurements (THEIA).In some of the more optimistic phase transition scenarios, future GW observations may even complement tests of Nnaturalness from next generation CMB experiments such as CMB Stage IV.

TABLE I .
Typical values of g * ρ and g * s for the SM, first SM-like sector, and first exotic sector for allowed parameter regions featuring a relatively large fractional energy density in the first exotic sector.