Probing exotic phases via stochastic gravitational wave spectra

Stochastic backgrounds of gravitational waves (GWs) from the pre-BBN era offer a unique opportunity to probe the universe beyond what has already been achieved with the Cosmic Microwave Background (CMB). If the source is short in duration, the low frequency tail of the resulting GW spectrum follows a universal frequency scaling dependent on the equation of state of the universe when modes enter the horizon. We demonstrate that the distortion of the equation of state due to massive particles becoming non-relativistic can lead to an observable dip in the GW spectrum. To illustrate this effect, we consider a first order chiral symmetry breaking phase transition in the weak-confined Standard Model (WCSM). The model features a large number of pions and mostly elementary fermions with masses just below the critical temperature for the phase transition. These states lead to a 20% dip in the GW power. We find potential sensitivity to the distortions in the spectrum to future GW detectors such as LISA, DECIGO, BBO, and μAres.


I. INTRODUCTION
The prospects for gravitational waves (GWs) to probe exotic particle phenomena have gained increasing attention in light of their discovery by the LIGO scientific collaboration [1].The proof of principle that they provided has led to the proposal for a large number of new GW telescopes covering a wide range of frequencies, from (10-100) nHz with Pulsar Timing Arrays (PTAs) [2][3][4] to the mHz to 10 kHz window using ground and space based telescopes [5][6][7].The range between PTAs and telescopes is challenging to cover, but there exist proposals such as the µAres detector [8] and through the use of asteroids [9].At much higher frequencies, in the MHz-GHz range, novel laboratory techniques have already constrained gravitational wave backgrounds [10,11] and there are plans for future detectors capable of extending current sensitivity further [12,13].These observatories will significantly broaden the scope of new physics that can be explored through gravitational waves in the coming decades.
This development has the potential to radically alter our understanding of the universe in much the way detailed studies of the Cosmic Microwave Background (CMB) [14] have over the past few decades.In studying the CMB [15][16][17][18], we were able to infer with great precision the properties of the early universe fluid by studying the spectrum of background photons, as well as the anisotropies therein.It would be remarkable indeed if the same level of detailed study could be performed by GW experiments.For the purposes of understanding the early universe, it is particularly interesting to learn what can be inferred from GWs emitted prior to big bang nucleosynthesis (BBN) [19].The CMB can only directly probe the universe back to recombination at T ∼ eV, but observations of the primordial abundance of light elements probes the universe back to T ∼ MeV.Above this, we thus far have no real probes and the best hope for accessing this era is in GWs; the universe is essentially transparent to them since they only interact via gravity and this makes them ideal messengers for probing the pre-BBN era.This is, however, not without challenges because the universe was really tiny at these times and effects like the Sachs-Wolfe effect [20,21], that have proven essential to determining our picture of the universe from CMB photons, are suppressed.
Several early universe events, such as first order phase transitions [22][23][24], hybrid preheating after inflation [25], dynamics of topological defects [26], and soliton collisions [27], can produce GWs in the pre-BBN era.These are referred to as "stochastic" and would constitute a nearly isotropic background of gravitational radiation just like the CMB is a nearly isotropic background of electromagnetic radiation.In this work, we focus for concreteness on GWs produced via a first order phase transition, though none of the qualitative results of our analysis rely on much beyond assuming that the GWs are produced nearly instantaneously and isotropically on cosmological times.
A first order phase transition proceeds by nucleation of the new vacuum phase as bubbles that expand.Eventually, those bubbles collide violently enough to provide a first source of GWs.Colliding bubbles can also percolate and interact with the surrounding thermal plasma, causing bulk motion in the fluid which, in turn, provides two further sources of GWs: acoustic (sound wave) and magnetohydrodynamic disturbances [28].The spectrum of GW emission from these processes have been simulated numerically and depends on multiple parameters such as the duration of the phase transition, the amount of latent heat released, and the fraction of energy which gets converted into the bulk motion of the plasma (see, for example, Refs.[28][29][30] for more detailed discussions).The precise spectrum of gravitational waves from these three sources is an area of ongoing research.Despite these model dependent features, a pair of related recent papers [31,32] pointed out the interesting fact that the low frequency tail of GWs from any "instantaneous" event would have a universal spectrum which only depends on features of the universe when those modes enter the horizon.Intriguingly, the spectrum is sensitive to the equation of state of the plasma.Provided the universe is in a radiation dominated phase, as expected in a standard cosmology, we expect w ≈ 1/3.Nevertheless, there can be small deviations from w = 1/3 and a careful measurement of the spectrum of GWs would be sensitive to such deviations.The equation of state has proven to be a useful probe of early universe phases, as illustrated in e.g.[33,34].
There are several potential sources such a distortion, but all require a breaking of scaling invariance.Any scale-invariant theory will have w = 1/3 exactly as where J λ is the scale transformation current.One could imagine several small sources of breaking in the radiation-dominated era.The recent works metioned above, Refs.[31,32], considered interactions inducing a scale anomaly, leading to w − 1/3 ∝ β, where β is the beta function of the theory, and free-streaming of particles that decouple from the plasma.
In this work, we focus primarily on scale-invariance breaking due to the masses of particles.We demonstrate that when several degrees of freedom simultaneously become nonrelativistic, the resulting distortion of the equation of state is sufficiently significant to be visible in future GW observations.Since there are hundreds of bosonic degrees of freedom at high temperatures, the resulting distortion goes like g/g * , where g is the number of degrees of freedom that are becoming non-relativistic and g * is the total number of relativistic degrees of freedom.Unlike the signals considered previously, this distortion is transient.If the mass of the particles is not too far from the temperature at which the GWs are sourced, then a localized distortion of the long-wavelength tail of the GW spectrum could be observed.
To demonstrate how this could work and could help determine qualitative information about the non-standard early universe cosmology, we consider the weak-confined standard model (WCSM) presented in Ref. [35,36].In that model, a scalar field coupling to the SU (2) L gauge kinetic operator leads to an effective shift in the scalar field coupling.This shift occurs when the scalar field has a non-zero expectation value, a mechanism first presented in Ref. [37].If the shift is toward a stronger SU (2) L coupling, then it is possible for the weak force to become strongly coupled at T > T c , the critical temperature for the electroweak crossover of standard cosmology.In a strongly coupled SU (2) L phase, B + L-violating transitions would occur readily, opening the door for a possible baryogenesis mechanism.For such a mechanism to work, the standard cosmology phase would need to be restored at T < T c so that the electroweak crossover does not washout the generated asymmetry.
More to the point of this work, the transition to strong coupling could be a first order phase transition, leading naturally to a source of stochastic gravitational waves.What makes the WCSM a particularly interesting source is that the nature of the phase of stronglycoupled weak force is in question.It is unclear whether this would be a chiral symmetrybreaking confined phase, as with standard model Quantum Chromodynamics (QCD), or an infrared conformal phase.In the latter case, the only distortions away from w = 1/3 would be due to the scale anomaly of the remaining interactions (QCD, electromagnetic, and Yukawa).On the other hand, the former would break scale invariance more strongly.A number of species would have masses near and a bit below the critical temperature, leading to modifications of the equation of state in the region below the phase transition.We show that in a chiral symmetry-breaking phase, pions becoming non-relativistic below the phase transition lead to more than 10% deviation from w = 1/3, correspondingly leading to an order 10% dip in the characteristic strain or 20% dip in the GW power at a frequency corresponding to modes that enter the horizon just as the pions are becoming non-relativistic.
For a WCSM phase transition above the would-be electroweak crossover, such a dip could yield potentially observable gravitational wave signatures.
In this work, we develop a calculation of the distortion of the gravitational waves due to particle masses in the WCSM.This requires modifying the results of Refs.[31,32] to account for a time-dependent equation of state.We solve the resulting gravitational wave equations numerically for benchmark models.We find potential sensitivity to the distortions in the spectrum to future GW detectors such as the Laser Interferometer Space Antenna (LISA) [38], Deci-hertz Interferometer Gravitational wave Observatory (DECIGO) [5], Big Bang Observer (BBO) [6], and the µAres detector.
The remainder of this paper is structured as follows.In Section II, we develop the formalism required to study time-dependent distortions of the equation of state and demonstrate the effect of a small number of non-relativistic species on the gravitational wave spectrum.
We then review the relevant details of the WCSM model in Section III.The results of our study are presented in Section IV.Finally, in Section V, we discuss the implications of our results for future measurements.

II. GENERAL FORMALISM A. Gravitational Wave Equation
Our starting point is to consider a Friedmann-Robertson-Lemaître-Walker universe with metric Following Refs.[31,32], we work in conformal time τ = dt/a(t) and expand the FLRW metric in small perturbations to obtain the linearized gravitational wave equation in comomving momentum space: where ′ denotes derivatives with respect to conformal time, H = a ′ /a is the conformal Hubble rate, k is the comoving wavenumber, and h is the metric perturbation amplitude.
Suppose that at a time τ * an instantaneous event (meaning that the duration of the event ∆τ * ≪ 1/H * , where H * ≡ H(τ ⋆ )) generates gravitational waves.Then one must add a source on the right-hand side of (3), and the resulting equation can be solved using Green's function methods, which is equivalent to solving the source-free equation with initial conditions where J * is related to the dimensionless anisotropic stress inducing the gravitational waves, projecting onto the GW polarization tensor ϵ.
When k ≪ H * , the modes are initially frozen outside the horizon.As they enter the horizon, they become sensitive the equation of state determining H.When k ≫ H, the modes lose sensitivity to the equation of state and simply oscillate.Thus, if a number of degrees of freedom are becoming non-relativistic at the time when a given mode is entering the horizon, the spectrum for that mode gets distorted.
We demonstrate that future experiments should have sensitivity to this distortion.To see this, we must first generalize the calculation of [31,32].In that work, the authors considered the equation of state to be approximately constant.For the effect we are interested in, this approximation is not accurate as the modification to the equation of state is happening on a Hubble time, which is the same as the time scale for the modes to enter the horizon.
Defining the equation of state w = p/ρ and plugging the first Friedmann equation H 2 = 8 π G ρ a 2 /3 into the acceleration Friedmann equation, we get an equation for a in terms of w: This equation can be solved to determine H as As w(τ ) can be determined from the matter content of the universe, we can solve the wave equation ( 3), at least numerically.
It will be helpful to change variables in equation ( 3), defining In terms of these variables, equation (3) becomes with primes now denoting derivatives with respect to τ .In terms of these variables, for a purely conformal equation of state w = 1/3, H = 1/τ , and the equation reduces to that of a harmonic oscillator with angular frequency 1.
At these late times, the solution goes like where ϕ is a phase shift induced due to the change in equation of state and A quantifies the change in amplitude relative to the case of no change in equation of state.The ensembleaveraged power is approximately a period-averaged power at this time, so that we can write independent of the phase shift.We can therefore examine A 2 (k) as an observable measure of the effect of the distortion of the equation of state.

B. Time-Dependent Equation of State
Each species in equilibrium with temperature T in the universe contributes an energy density and a pressure where + is for fermions and − is for bosons, g i is the number of degrees of freedom corresponding to the species and x = m i /T for a species of mass m i .The total energy density and pressure are obtained by simply summing over the species, From this, we can determine the equation of state parameter w as a function of temperature via Then, we can calculate δw as a function of temperature using To illustrate the change in the equation of state, we first consider a toy model.We consider a scenario in which a thermal bath contains 2 different bosonic species.One of them has To write w as a function of τ , we need the relation between conformal time and temperature, obtained by combining the first Friedmann equation with the conservation of entropy, We can then obtain the solution The entropy is given by the usual Euler relation Eqs. ( 18) and ( 19) allow us to express the conformal time as a function of temperature.
This relation can be inverted to get the temperature as a function of conformal time, which allows us to express δw as a function of τ .

C. Gravitational Wave Model
The phase for the gravitational wave is arbitrary and depends on the random initial conditions of the universe.To estimate the sensitivity, we therefore consider the ensembleaveraged power in the gravitational waves, In Fourier space, we define the power spectrum by From this, we can define the differential power by As discussed before, for a first order phase transition there are three contributions to the gravitational wave spectrum, from collisions of the bubbles themselves and from acoustic and magnetohydrodynamic distrubances induced by their expansion in the thermal plasma.
The three sources are expected to combine linearly, To determine which one of these sources dominate we must consider the system undergoing a first order phase transition.In a strongly interacting theory such as the WCSM, one can reasonably expect that the bubble wall interacts sufficiently with the thermal plasma to cause bulk motion of the latter.In this scenario, numerous numerical studies (see for example [29,39]) have found that sound waves are the dominant source of GWs compared to the collisions of bubbles themselves.[40] The peak frequency is given by [29,30] f sw = 1.9 where β is the inverse duration of the phase transition, H * and g * are respectively the Hubble parameter and number of relativistic degrees of freedom at a temperature T * .The latter is usually taken to be the nucleation temperature T n when the bubble nucleation rate is one per Hubble volume since that is when the production of GWs from all three sources is most significant.The dimensionless gravitational wave spectrum today obtained numerically in [39] is well fitted by [29] where α is the ratio of the vacuum energy released during the phase transition to the energy of the radiation bath, and κ is the fraction of the energy of the vacuum energy that gets converted into the bulk motion of the fluid.Aside from α and κ, the density is parameterized in terms of other properties of the phase transition such as the ratio of the inverse duration of the phase transition to the conformal Hubble at the phase transition H * , the bubble wall velocity v w , and the number of relativistic degrees of freedom at T * .
We consider only long wavelength waves, with f < f sw .The f 3 scaling is a universal property of short duration phase transition.Rather than directly calculate this power spectrum in our model, we perform a simple rescaling by where we define the averaging procedure here by averaging over a period in τ of 2π at τ → ∞, when the gravitational waves are well within the horizon and lose their sensitivity to the equation of state.We denote by ⟨ h2 ⟩ 0 the averaged gravitational wave in a purely radiation-dominated universe.Taking the now-arbitrary normalization J * = 1, we find The goal for much of the remainder of this paper is then to solve the equation ( 9), with the initial conditions h(k τ * ) = 0 and h′ (k τ * ) = τ * .
Before going into that calculation, we note on the observable spectrum in the present day universe.The gravitational waves are a bit colder than the radiation, as they decouple immediately after the phase transition.The spectrum gets scaled by a factor of a(τ * ) −4   and we need to convert from co-moving wave number to physical frequency.We work in a convention with a(τ * ) = 1, in which case f = k/(2 π a(τ 0 )).

III. PARTICLE PHYSICS MODELS
The  scalar is coupled to the kinetic term for the gauge bosons, where Λ is the scale of the operator.If the scalar Φ acquires an expectation value in the early universe, the value of the weak coupling g ′ is shifted.Since the weak coupling constant is asymptotically free, it will run to strong coupling if the nominal confinement scale Λ W is larger than the elecotroweak phase transition temperature T c ≲ v.In order to have Big Bang Nucleosynthesis as observed, the weak coupling should return to the SM phase at a temperature T SM ≳ 10 MeV.Further details of the WCSM can be found in Ref. [35]. The Out of 9 gauge bosons, 5 gauge bosons become massive and five pions are eaten.Out of the remaining 60 pions, 58 pions acquire masses through gauge interactions and Yukawa interactions.Four of these pion masses are tachyonic, but we expect these pions to be light after minimizing their potential.
In addition to the massive pions, there are fermionic states in the WCSM which contribute to modifying the equation of state.The left-handed fermions of the SM become mostly composite, while the right-handed states are mostly elementary.The composite fermions have masses around the confinement scale Λ W ∼ 4πf π .Only the elementary fermions, with masses ≪ Λ W , contribute significantly to modifying the equation of state at relevant temperatures.
The mass spectrum of the pions and the fermions is sensitive to dimensionless operator coefficients which encode loop corrections and are expected to be O(1) numbers.The values of these coefficients can be taken to be ±1 for the purpose of calculating the numerical spectrum of the WCSM given in Table I.
With this numerical spectrum, we can repeat the procedure outlined in Sec.II B to calculate δw for the WCSM.The behaviour of δw over a range of values of T for a benchmark value of f π = 80 TeV is shown in Fig. 2. For the gravitational wave spectrum, we illustrate three benchmark values for f π = (0.08, 8, 80) TeV, which sets the critical temperature of the chiral phase transition T C ≃ f π .
For the gravitational wave parameters in (25), the precise value of β/H * is obtainable only through a numerical lattice study of the WCSM.That study is beyond the scope of this work and we therefore adopt three benchmarks, β/H * = 1, 10, 100, noting that since we are assuming that the phase transition occurs quasi-instantaneously, β/H * = 1 is a lower limit.Furthermore, we adopt two benchmark (BI and BII) values for α and κ.BII assumes values for the latter that are fully consistent with the numerical simulation of [39] while BII is an extrapolation of the latter assuming that the GW spectrum scaling with respect α and κ is unchanged.For both benchmarks, we assume the wall velocity to be v w = 0.83 c, the upper limit considered [39].With these, we show the GW spectrum of an assumed first order WCSM phase transition and compare them to the reach of future experiments like LISA, DECIGO, and BBO.We also show the expected GW spectrum from galactic binaries, the dominant background for our model.The dotted lines in the spectrum are used to describe modes entering the horizon when the validity of the equation of state is Thor (GeV)  Thor (GeV)  Thor (GeV) A re s μ 10 1 10 7 10 3 10 5 Thor (GeV) The stochastic background of GWs from galactic binaries is shown in gray.The dotted lines in the spectrum denote frequencies with model-dependent effects.The left panels show the spectra for BI and the right panels the spectra for BII as described in the text.

IV. RESULTS
We have solved (9) numerically in the WCSM phase with initial conditions described in Section (II).As discussed above, the particle physics effects of the WCSM on this equation As such, we trace the high frequency end of our spectrum within an order of magnitude of the cutoff scale Λ χ = 4πf π for chiral perturbation theory, as a dashed line.Furthermore, any modes that would not quite reach a point where they lose sensitivity to the equation of state before the universe deconfines back to the SM phase are subject to modification by that deconfining phase transition and changes to the equation of state just below that phase transition.Therefore, for which 10π/k < τ SM , the conformal time of the transition to the SM phase, are also drawn as dashed lines.
In Figure 4, we plot the ratio of the time-averaged gravitational wave amplitude today for δw ̸ = 0 compared to the δw = 0 case for f π ≃ 80 TeV.Since this value of f π corresponds to the variation of the equation of state with respect to temperature (upper axis) shown in GeV and corresponds to the mode that enters the horizon when particles of mass around f π (mainly 49 pions which get their mass predominantly from the strong coupling or the top Yukawa) become non-relativistic.
We also show a zoomed in version of the f π = 80 TeV, α = 0.5, κ = 0.1, β/H * = 1 benchmark (Fig. 3f) in Fig. 5.To make the f 3 scaling at low frequencies apparent, we plot h 2 Ω gw /f 3 instead of h 2 Ω gw .We show the δw = 0 case, corresponding to no conformal symmetry breaking, and the WCSM with chiral symmetry breaking, where conformal symmetry is strongly broken.We see that LISA has some sensitivity to the two curves in a region where the conformal and non-conformal scenarios differ which lies at the upper end of the Thor (GeV)  region where the causality limited f 3 scaling is valid.There are also pronounced dips well into the region where the f 3 scaling is valid, but this lies in the µHz gap described in Ref. [9].Unfortunately, the sensitivity of the proposed detection scheme in this work still does not quite reach the sensitivity required to see the deviation from the expected f 3 scaling for this benchmark, but we hope that our model adds further motivation to thinking about future detectors capable of probing this frequency range, where there is already significant interest for astrophysical sources of GWs.3f with only β/H * = 1 benchmark case shown.To make the f 3 scaling at low frequencies apparent, we plot h 2 Ω gw /f 3 instead of h 2 Ω gw .The blue curve represents the δw = 0 case, corresponding to no conformal symmetry breaking, the dark blue curve represents WCSM with chiral symmetry breaking.We see that LISA and, especially, µAres have some sensitivity to the two curves in a region where they defer and the universal f 3 scaling is valid, making such an effect potentially observable.

V. DISCUSSION AND OUTLOOK
The signal we describe here is quite general and applies to any scenario with a long wavelength gravitational wave background and a significant number of particles becoming non-relativistic.For example, one could examine supersymmetric scenarios with a large spectrum of superpartners at similar masses.The signal also bears some resemblance to that studied by Ref. [41].One could imagine slight tweaks of that scenario to get similar dips.
One further intriguing possibility is to modify the QCD phase transition such that it is first order.One scenario in which this would occur is a supercooled electroweak phase transition [42], as the QCD transition would then occur with six light flavors of quarks.
After the QCD phase transition, the amount of radiation in the universe is rather low.As such, the process of the muons becoming non-relativistic is sufficient to induce a roughly 5% change in w, which is more challenging to observe.
While the prediction of a f 3 scaling of the power, up to modifications due to the equation of state, is universal, the normalization of the spectrum and peak frequency are dependent on the properties of the phase transition.Further calculations are needed to determine these properties of the phase transition.Lattice gauge theory is rather challenging for this model due to the addition of the Higgs field, but this would be the most robust way to study the phase transition.A lattice study would also elucidate the nature of the strongly interacting phase, which we assume be chiral symmetry breaking as in QCD.It is possible that the gravitational wave properties can be also estimated from a Nambu-Jona-Lasinio model [43,44].
Note that any dimensionful parameters in the model could in principle lead to such a deviation in the spectrum.In particular a super-renormalizable scalar cubic coupling or non-renormalizable operator will lead to a small shift in the power of the low frequency spectrum.
The promising signal type presented in this work will ultimately require further study by gravitational wave telescopes.

Figure 1 .FIG. 1 . 4
Figure 1.In this toy model, a peak change in w is found at T = 0.38m with δw = 0.03.

2 FIG. 2 .
FIG. 2. Deviations from w = 13 in the WCSM for a benchmark scenario where f π ≃ 80 TeV.The peak near T /f π ≃ 1 corresponds particles of mass around f π -mainly 49 pions which get their mass predominantly from the strong coupling or the top Yukawa -decoupling from the thermal bath.

FIG. 3 .
FIG.3.The GW spectrum of an assumed first order WCSM phase transition compared to the reach of future experiments like LISA (purple), DECIGO (orange), BBO (green), and µAres (brown).
describing the propagation of GWs are encoded in the small changes to the Hubble parameter due to deviations of the equation of state from w = 1/3 during the radiation domination era.The latter is determined assuming the validity of chiral perturbation theory, but there is some uncertainty, inherent to the strongly coupled nature of this theory, on what the critical temperature for the phase is, what the properties and nature of the phase transition are, and what distortions of the equation of state would occur due to relatively strong pion interactions close to the confinement scale.A robust answer to these questions can only come from a detailed numerical lattice study, which is beyond the scope of this work.

Figure 2 ,
Figure 2, we can compare the two figures and see that the largest change occurs at k ≲ 10 −9

FIG. 4 .
FIG.4.Average gravitational wave power relative to that in the absence of changes to the equation of state as a function of comoving wavenumber in (GeV).The vertical axis is the temperature of the thermal bath at the moment of horizon entry for each mode.The largest change -the dip at k ≲ 10 −7 GeV corresponds to the mode that enters the horizon when particles of mass O(f π )mainly 49 pions -decouple from the thermal bath.

FIG. 5 .
FIG. 5. A zoomed in plot Figure3fwith only β/H * = 1 benchmark case shown.To make Weak Confined Standard Model (WCSM) is a model in which the SU(2) L component of the electroweak force is strongly coupled and weak isospin is confined.To achieve this, a

TABLE I
. Numerical spectrum of the WCSM at f π = 80 TeV.Additionally, there are 4 light tachyonic degrees of freedom as well.