Causal gravitational waves as a probe of free streaming particles and the expansion of the Universe

The low frequency part of the gravitational wave spectrum generated by local physics, such as a phase transition or parametric resonance, is largely fixed by causality, offering a clean window into the early Universe. In this work, this low frequency end of the spectrum is analyzed with an emphasis on a physical understanding, such as the suppressed production of gravitational waves due to the excitation of an over-damped harmonic oscillator and their enhancement due to being frozen out while outside the horizon. Due to the difference between sub-horizon and super-horizon physics, it is inevitable that there will be a distinct spectral feature that could allow for the direct measurement of the conformal Hubble rate at which the phase transition occurred. As an example, free-streaming particles (such as the gravity waves themselves) present during the phase transition affect the production of super-horizon modes. This leads to a steeper decrease in the spectrum at low frequencies as compared to the well-known causal $k^3$ super-horizon scaling of stochastic gravity waves. If a sizable fraction of the energy density is in free-streaming particles, they even lead to the appearance of oscillatory features in the spectrum. If the universe was not radiation dominated when the waves were generated, a similar feature also occurs at the transition between sub-horizon to super-horizon causality. These features are used to show surprising consequences, such as the fact that a period of matter domination following the production of gravity waves actually increases their power spectrum at low frequencies.


Introduction
With the first direct detection of Gravitational Waves (GWs) [1], we have gained access to a new probe of the universe. One especially exciting aspect of GWs is that, because they are so weakly interacting, they can carry information about the early universe in a very clean and unpolluted form. The future of GWs is bright as there are many proposed future detectors such as LISA [2], BBO [3], MAGIS [4] and DECIGO [5]. Using stochastic GW background searches, these new GW detectors can probe the universe at a time 10 28 times earlier than our current best direct probe, the Cosmic Microwave Background (CMB).
There are many exciting possibilities that stochastic GW backgrounds will allows us to explore. Stochastic GW backgrounds can be produced by a myriad of early universe phenomena such as inflation [6][7][8][9], reheating/preheating [10][11][12][13][14][15], phase transitions [16] and topological defects [17,18] such as cosmic strings or domain walls, or by scalar perturbations at second order in perturbation theory [19][20][21][22][23][24]. The detection of such a background would teach us details about which of these phenomena were active and what were the relevant parameters governing their dynamics.
Aside from learning about the actual source of GWs, their discovery also teaches us about the dynamical evolution of the early universe. There exists a vast literature discussing these effects and their detectability in various scenarios [16,[25][26][27][28][29][30][31][32][33][34][35][36][37][38][39][40][41][42] as well as reviews on the subject [17,43]. Intuition for this dependence can be most easily developed focusing on a situation that has been largely overlooked, the low frequency spectrum of a stochastic GW background produced by a causal source. This is the case for GWs generated by a phase transition or parametric resonance for example. The low frequency modes that correspond to length and time scales much longer than the source duration are largely independent of the details governing the generation of GWs and are instead fixed by causality. The largest region in causal contact with itself at the time GWs were generated has a size at most of the order of the Hubble radius, but is generically smaller. These low frequency modes can only be affected after entering horizon, making the spectrum sensitive to the change of the Hubble rate as a function of time, or in other words, to the expansion history of the universe.
In this paper we consider the low frequency part of the power spectrum of a stochastic background of GWs. These modes are typically referred to as the causal part of the spectrum in order to emphasize the effect of causality, and under standard assumptions have a k 3 scaling [44]. Nonetheless, as we will see later, there are effects other than causality that determine their dynamics. We will refer to these modes as the "causality-limited" part of the spectrum. The shape of the spectrum for modes which are causality-limited but were subhorizon at generation (sub-horizon causality) is completely independent from the evolution of the universe 1 . Causality arguments affect modes which were super-horizon at generation in a two-fold way. The first factor is that the production of these modes is highly suppressed. The reason for this is that they have frequencies ω H (H being the Hubble rate) at the time of production, and therefore behave like an over-damped harmonic oscillator, making such modes very hard to excite. As a result, this effect decreases the power at low frequencies.
The second factor is that once excited, their amplitude stays constant while they are superhorizon, rather than decreasing, leading to an overall increase in the power at low frequencies.
The competition between these two effects determines the shape of the low frequency tail of GWs.
More mathematically, the physical intuition that we will develop is that the initial conditions for causality-limited modes immediately after generation are h ij = 0 whileḣ ij = 0, in much the same way that the initial condition of a hammer hitting a string creates no displacement but a large velocity. If the GW is sub-horizon, ω H, the mode obtains an amplitude h ∼ḣ/ω. If the mode is super-horizon, friction is important and the mode obtains a smaller than expected amplitude h ∼ḣ/H. After this initial displacement, over-damped modes are frozen in place until Hubble crossing at which point they become the standard under-damped modes seen in GW detectors.
There are several interesting observations that result from an analysis along the previous lines. The first is that generically there is a difference between sub-horizon and super-horizon causality. As a result, there will generically be a feature at horizon crossing. Finding such a feature would give a model-independent direct measurement of the conformal Hubble rate at the time when GWs were created. Upon specifying an expansion history, it is possible to deduce the Hubble rate, and hence the temperature, at which the GWs were generated. We will discuss two scenarios in which GW modes that are super-horizon or sub-horizon (but at a wavelength larger than the correlation length of the source) at the time of generation behave differently.
In the first, we study how super-horizon GWs depend on the equation of state parameter w. We show that while sub-horizon causality forces the spectrum of GWs to fall as k 3 , that super-horizon modes fall as k (1+15w)/(1+3w) . Thus, as long as w = 1/3, the scaling of GWs changes between modes that are sub-horizon at the time of generation versus super-horizon. The change in behavior between sub-horizon modes and super-horizon modes during matter domination (MD) was previously studied in Ref. [45] and followed up in Refs. [46][47][48]. We improve upon these results by developing simple physical intuition for these effects and by generalizing them to arbitrary equations of state.
For the second scenario we look at the effect of free-streaming particles on the production and propagation of GWs. This was originally studied in Ref. [49] for the case where some species of particles started free-streaming after GWs were produced, and it showed that this leads to a constant suppression in the amplitude as modes entered the horizon. We focus on the contrasting case where free-streaming particles were already present during the creation of the GWs. We find that sub-horizon modes have the standard k 3 scaling. Super-horizon modes have a scale dependent suppression, and for a sufficiently large fraction of free-streaming particles there is a a surprising oscillatory feature on top of a suppressed scaling of k 4 . These effects come about because in this scenario super-horizon GWs are not frozen out but instead slowly roll in the potential of the free-streaming particles, suppressing the power in GWs. This effect is always present, even if at a small level, as the high frequency part of the GW spectrum itself acts as free-streaming particles.
Another interesting result we find is that an intermediate period of matter domination after the GW generation actually increases the power at low frequencies. There are well motivated scenarios where various effects due to reheating generate GWs [10,11]. In the typical scenario where reheating occurs at high scales, only the causal part of the GW spectrum is accessible as the peak frequencies are too high. An additional feature of many of these models is that they generically imply matter domination after the production of GWs. This effect was previously considered to be undesirable as it dilutes the peak power of the GW. However, for the low frequencies that are experimentally visible, the effect of a period of matter domination is to increase the visibility of these scenarios, making it a desirable feature.
In Sec. 2, we discuss the physical intuition behind causality-limited GWs and compare our physically motivated approximations against exact results. In Sec. 3, we study the effects of free streaming particles on causal GWs. In Sec. 4, we study how the frequency spectrum of gravitation waves changes as the universe transitions between matter, radiation, and subhorizon modes. Finally, we conclude in Sec. 5.
2 The spectrum of causality-limited gravity waves

The physical intuition
Throughout this work we will be exploring "causality-limited" GWs, that we define as follows. Consider a source that is active for a short amount of time 1/β. The causality-limited part of the GW spectrum consists of the waves whose period and wavelength are much longer than the source's temporal and spatial correlations respectively. In other words, λ 1/β, where λ is the wavelength and 1/β is the duration of the process generating the GWs. The standard example of such causality-limited GWs is the low frequency part of the GW spectrum generated by a cosmological phase transition.
Using conformal time τ and conformal Hubble rate H = a /a, the equation of motion for a comoving mode k of the graviton h ij is (we follow a notation close to the one of [44]) where Π ij is the energy momentum tensor of the sector generating the GWs. We will project both h ij and J ij onto the two independent (+, ×) polarisations, and assume that the respective amplitudes J and h of the two polarisations are equal. We will also assume that the GWs are produced on a time scale fast compared to Hubble so that we can approximate H as a constant over its production and the time over which Π is non-zero to be small. The time at which the phase transition occurs will be denoted by τ . The solution to Eq. (2.1) can be found using Green's functions to be h(k, τ ) = dτ e H(τ −τ )/2 We will be interested in the initial conditions right afterwards τ = τ + assuming that J occurs fast so that we can take J(k, τ ) = J (k)δ(τ − τ ), and we will drop the k dependence from here on. Using this approximation and Eq. (2.2), we find that Afterwards, we have the second order differential equation and initial conditions The observable relevant for GW detectors such as LISA is d Ω gw /d log k. This is obtained using We will also take with P h (k, τ ) the dimensionful power spectrum of GWs. Fourier transforming gives where ρ c = 3H 2 /(8πG) is the critical energy density. In the second equality we use the simplification valid at late times that h = k h, and where an additional k 3 comes from the phase space factor. We can then obtain the k scaling of d Ω gw /d log k by taking k 5 and multiplying by the k scaling of h 2 .
The observable of interest depends on h(k, τ ), which comes from solving Eq. (2.4). From that we can see that there are two contributions to the k scaling of hh , the k dependence of the source J (k) and propagation effects. Here is where causality plays a central role: as long as we are interested in wavelengths much longer than the typical spatial correlation of the sources, the Fourier transformed source J ≈ constant (see e.g. [44]), and hence the k dependence coming from the source is trivial for long wavelength modes. 2 Assuming correlation lengths for the source, k −1 source , that are small compared to the horizon size, we can separate the solutions of the causality-limited modes into two regimes, k source k H (sub-horizon) and k H (super-horizon). Note that in the case of the phase transitions, the existence of a causal sub-horizon regime depends on whether all the sources of GWs (including in particular the sound wave contribution) have a short duration compared to the Hubble time (see e.g. [50][51][52][53] for recent studies on this topic).

Sub-horizon regime k
H : These modes are under-damped by definition. Thus, the standard approximation is that of a frictionless solution rescaled by a(τ )/a(τ ) to take care of Hubble friction as can be seen by using the WKB approximation. This is easily shown to be From this, we immediately see that sub-horizon causality automatically gives regardless of what the equation of state of the universe is.
Super-horizon regime k H : There are two competing factors controlling the dynamics of the k H modes. The first is that one is attempting to excite an over-damped harmonic oscillator. This effect suppresses the power in these modes. The second is that 2 More explicitly, if for separations x larger than a length scale λ the two-point function PΠ(x) of the source vanishes, then its Fourier transform super-horizon modes are frozen in place until they enter the horizon, increasing the relative power in these modes. The competition between these two effects is what controls the behavior of k H modes. Very quickly (roughly after a single e-fold) the initial conditions due to Hubble friction, and remain frozen while the mode is super-horizon. The O(1) number out front depends on the exact equation of state of the universe. This scaling can be confirmed by solving ∂ 2 τ h ij + 2H ∂ τ h ij = 0 with a non-zero initial velocity. Notice that the amplitude of the oscillation is suppressed compared to what one would expect based on the amplitudes for sub-horizon modes: ∼ J /k. This suppression is the effect of attempting to excite an over-damped harmonic oscillator.
At this point, the mode is frozen out and its amplitude is constant horizon entry, which occurs at a conformal time τ k given by After the mode enters the horizon, the GW can be approximated as completely under-damped and redshifts away as We have been cavalier about the phase information just putting in sin kτ , as the phase is unimportant for computing the power due to averaging over many oscillations.

Radiation domination
It is well known that in radiation domination (RD) there is a simple solution for the GWs. During RD, a ∼ τ and H = a /a = 1/τ . We can solve Eq. (2.4) exactly in this case and find Thus, the full solution is given by the solution in a friction-less universe times a red-shifting factor of a(τ )/a(τ ). We now show that the physical intuition given in Sec. 2.1 replicates this exact solution. The sub-horizon modes given in Eq. (2.8) have the same scalings as the exact solution. The super-horizon modes given in Eq. (2.13) simplify using which is again the same scaling as the exact solution, confirming our physical intuition. From these results we can calculate how dΩ gw /d log k scales with k: This reproduces the famous fact that causality in a RD universe leads to the spectrum of GWs falling off as k 3 in the small k limit.

General equation of state
We now repeat our exercise for a general equation of state p = wρ. In this case, we have with n = 2/(1 + 3w), so that for radiation (matter) domination we have n = 1 (n = 2). In this case the exact solution of Eq.
where j n (x) and y n (x) are spherical Bessel functions of first and second kind respectively. While this is the exact solution, it is useful to consider two limits. The first is the limit of sub-horizon modes: The second is the limit of super-horizon modes at production after they enter the horizon, Again, we will show that the physical intuition given in Sec.
After some algebra, we find that reproducing the parametric behavior of the exact solution for the super-horizon modes. The scaling of dΩ gw /d log k can be directly obtained from the scaling of the solutions. For sub-horizon modes k H As expected, as long as the wavelength of the mode is sufficiently larger than the distances over which the sources are spatially correlated, sub-horizon causality forces these modes to fall as k 3 . For the super-horizon modes, we have We see that there is more power at low k for any equation of state with w < 1/3. In particular, for the case of MD, n = 2, we find that From the discussion above we find that there is a distinct kink in the spectrum at horizon crossing as the modes transition from sub-horizon to super-horizon causality. This kink in the spectrum occurs for all cosmologies, except for exact RD where the spectrum scales as k 3 for long-wavelength modes. Identification of this generic feature in a signal allows for a model independent measurement of the value of conformal Hubble at which the GWs were generated.
We can use the exact analytic results of Eq. (2.18) to obtain dΩ gw /d log k. To do this, we take the late time limit of Eq. (2.18) and find From the equation above we quickly see that where φ is a constant phase factor. Finally, we arrive at an exact expression for the energy in GWs (after averaging over the oscillations) where we have neglected unimportant proportionality constants. In Fig. 1 we use these exact analytic results to show how sub-horizon scaling becomes super-horizon scaling for various equations of state.

The effect of free-streaming particles
In this section, we show how free-streaming relativistic particles induce not only a dampening of GWs at horizon crossing but also change the shape of the causality-limited part of the GW spectrum. We specialize to the case of a radiation-dominated universe with a fraction f fs = ρ fs /ρ total in free-streaming relativistic particles since for other equations of state the fraction f fs changes significantly with the expansion of the universe. As was derived in Ref. [49], GWs are affected by the presence of particles whose free-streaming length is larger than the Hubble radius. In this case, the propagation of the free-streaming particles along geodesics is affected by and affects the propagation of GWs.
It is well known that during BBN and CMB a large fraction of the energy density in radiation was free-streaming, since neutrino interactions freeze out at temperatures below approximately 2 MeV. Little is known about our Universe prior to this era, and therefore the existence of other relativistic free-streaming species at earlier times remains an interesting possibility in the early Universe. In fact, since GWs themselves are free-streaming, the high frequency part of the GW spectrum itself is an irreducible contribution to the population of free streaming particles. We will see that if such particles were present during the generation of GWs, they would lead to striking features in the causality-limited part of the spectrum.
The equation of motion for GWs in this scenario is [49] where u = kτ . As indicated before, we are interested in causality-limited GWs. Thus, we are solving these equations subject to the boundary condition There are two important effects that free-streaming particles have on GWs: post-production dynamics and horizon entry. The effects of free-streaming particles on horizon entry is the standard effect emphasized in the case of neutrinos. Note that the LHS of Eq. (3.1) is larger than the RHS by a factor of u 2 and hence the largest effect is expected to be at early times when u = kτ 1. From this we see that sub-horizon modes are completely unaffected by free-streaming particles. Conversely, super-horizon modes are hit by a uniform frequency independent suppression when they enter the horizon, e.g. GW amplitudes are suppressed by a factor of 0.8 when f = 0.4, as originally shown in Ref. [49]. Amusingly, we will see that for GWs generated by phase transitions, this frequency-independent suppression is in fact a subdominant effect.
The novel effect here is that free-streaming particles also cause a suppression in the production of GWs. What we mean by this is that the transition of super-horizon GWs from a large velocity with a small amplitude to a large amplitude with no velocity discussed in Eq. (2.11) is significantly affected by the presence of free-streaming particles. The suppressed production can be computed analytically using the techniques of the previous sections.
We first consider the limit u = kτ 1, where the modes are very super-horizon. In this limit, Eq. (3.1) can be simplified to The only difference between this equation and Eq. (3.1) is that we have taken the kτ 1 limit of the RHS, eliminating the integral. This approximation allows us to study the effect of free-streaming radiation on the evolution of super-horizon GWs. The approximation breaks down near horizon entry, however the effect of horizon entry has been studied in the past and approximately amounts to a k-independent reduction of the amplitude, proportional to f fs for f fs 1, and so it can be treated separately. The suppression of GW production by free-streaming particles can be directly seen from the RHS of Eq. (3.3), where the free-streaming particles act like a Hubble scale mass. Thus, in contrast to the normal scenario where wave amplitudes freeze out and remain constant while outside of the horizon, the effective Hubble scale mass is constantly reducing the amplitude of the GW due to a slow roll type effect. Note that this effect is absent if the wave starts out with negligible velocity like in previously studied scenarios, in which cases the RHS vanishes.
Eq. (3.3) can be solved exactly with our initial conditions to give where j α and y α are spherical Bessel functions and we have taken the large time limit in the second line. From this, we find that the suppression factor is simply which follows directly from Eq. (3.4) (and holds even when α is imaginary). While Eq. (3.5) is rather un-illuminating in and of itself, it can be simplified in two interesting limits 4f fs 1 and 4f fs > 1. In the first limit, we find that In the other limit, 4f fs > 1, we find where C 1 , C 2 , C 3 are unenlightening functions of f fs . Interestingly, we find that in this limit there is an overall suppression of the amplitude that is independent f fs and that there is an additional oscillatory feature on this suppressed amplitude. The appearance of an oscillation as f fs increases comes from when the super-horizon modes go from being over-damped to under-damped even while super-horizon. In this limit, the mass coming from free-streaming particles overcomes Hubble friction and induces oscillations. To obtain the shape of the GW spectrum as a function of k, we numerically solve Eq. (3.1). The results for the suppression can be found in Fig. 2, which shows the numerical results as well as the analytical approximation, Eq. (3.5). As can be clearly seen, the analytic results and the numerical agree very well for f fs = 0.01, 0.1 and 0.4. This agreement holds despite the fact that the transition region, where u = kτ ∼ 1, is not accurately captured by our approximation, and hence, the traditional effect of horizon entry is not entirely taken into account.
We can re-express our results for the amplitude of the GWs in terms of the observable  for large free-streaming fractions, where C 1,2,3 are constants. This shows that the presence of free-streaming radiation during the generation of GWs leads to a spectrum that decreases faster at low frequency, with the spectral index growing linearly in f fs from 3 until it saturates at 4 for f fs ≥ 5 32 . As mentioned before, the high k (sub-horizon) modes of the GWs themselves necessarily play the role of free-streaming particles. The derivation of Eq. (3.1) assumes that the freestreaming particles follow geodesics, which GWs do as well. Thus, the only possible difference between the high frequency GW modes and free-streaming particles is that the high frequency GW modes were produced by the phase transition itself. Imagine that the GWs are produced at a time after the phase transition, then Eq. (3.1) becomes As long as H 1, we recover Eq. (3.3) and we see that the high frequency GW modes behave exactly like any other free-streaming particles.

Implications of current N eff limits
Results from CMB and BBN measurements limit the amount of non Standard Model radiation present in the Universe at those times. These results are presented as limits on the effective number of neutrinos, where ρ γ and ρ ν are respectively the energy density in photons and neutrinos and ρ X represents any new contribution to the energy density in radiation.
A combined analysis of BBN and CMB measurements leads to a bound ∆N eff 0.3 [54]. This implies that at neutrino decoupling a new type of radiation could contribute to at most 4.4% of the radiation energy. However, due to the decoupling of massive degrees of freedom, the relative temperature between the standard model plasma and any decoupled species changes as the Universe evolves. In particular, a ∆N eff 0.3 at BBN translates to f fs 0.09 for temperatures above the weak scale.
This shows that if below the weak scale there are no modifications to ΛCDM, the maximum allowed f fs is 0.09, and if a non-zero N eff is measured in future CMB experiments, the nature of this radiation at high energies could be tested using GWs. 3 Nevertheless larger f fs are possible with simple modifications of early cosmology, such as late equilibration of the radiation with the SM with subsequent decay or through an entropy injection into the standard model plasma, which dilutes the relative energy density in the new species. The latter case would also lead to a temporary change in the equation of state that changes the GW signal as discussed in the next section.

Effect of non-standard expansion histories
In this section we show how the causality-limited spectrum of GWs can be used to probe the expansion rate of the Universe after the GWs were generated. There are two main effects of a non-standard expansion history in the spectrum of GWs: an overall re-scaling of the spectrum due to the change in expansion history; and a change in the power-spectrum of causality-limited modes that enter the horizon during an era of non-radiation domination.

An intermediate period of matter domination
In this subsection we consider the effect of an intermediate period of matter domination (MD). We will compare two scenarios. In the first scenario, GWs are produced at a temperature T and the universe is radiation dominated (RD) until CMB. In the second scenario, GWs are again produced at a temperature T but there is an intermediate period of matter domination between the temperatures T r→m , where the universe goes from RD to MD, and T m→r , where it goes from MD to RD. A summary of the results in this section can be seen in Fig. 3 where we compare the numerical solutions for the spectrum in both cases. We find that if one is interested in modes that enter the horizon before T r→m , then the extra red-shifting induced by the period of MD simply decreases the amplitude of these modes. If one considers instead the modes that enter the horizon after T m→r then their amplitude increases due to the red-shifting of the frequencies moving power from high frequencies to lower frequencies. In what follows, we give an intuitive explanation of the results while in App. A we give a more analytical derivation.
Matter domination gives an entropy dump that increases the red-shift between the source and present day observers yielding two distinct effects on the spectrum. The first is that the extra expansion dilutes the energy density, decreasing the overall power. The second is that the frequencies themselves are red-shifted, an effect that increases the power at low frequency as it shifts the entire spectrum to lower frequencies. The competition of these two effects determines the behavior of the spectrum. The net effect of MD can be understood intuitively and is shown visually in Fig. 4. We first discuss the effect of the dilution of the energy density from the additional red-shifting. For modes that enter the horizon before T r→m , the overall power is redshifted by an additional (∆a) −4 due to the additional expansion. Here we take ∆a = (a /a 0 ) 1 (a /a 0 ) 2 = T r→m T m→r to be the difference between scenarios 1 (pure RD) and 2 (intermediate MD) in the amount of red-shifting between the phase transition and us. For modes that re-enter the horizon after T m→r , the modes were frozen out during matter domination and so their overall power is not modified. Intermediate frequencies interpolate between the two results with a f scaling. We now discuss the effect of the red-shifting of momenta. The observable GW frequency f corresponds to the physical momentum k/a of the GW, which in scenario 2 redshifts more by a factor ∆a. The f 3 fall off of the power spectrum during RD means that when comparing the two spectra at the same physical momenta, the power is enhanced by a factor of (∆a) 3 for all modes entering during RD. In totality, for modes that re-enter the horizon before T r→m , these two effects combine to give a total suppression of 1/∆a while for modes that re-enter the horizon after T m→r , they combine to give a total enhancement of (∆a) 3 . Thus we find that an intermediate period of MD increases the power at low frequencies but suppresses it at high frequencies.
These estimates can be made more concrete as follows. In the first scenario where the universe was always radiation dominated, the spectrum of GWs is where A gw,1 is the amplitude of the GW spectrum at its peak and this scaling holds up to the peak frequency f ,1 . Frequencies higher than this are sensitive to the properties of how the signal was generated. As discussed before, there are two main effects of a period of intermediate matter domination. Firstly, the amplitude of the peak signal is decreased by (∆a) −4 , Secondly, the frequencies are shifted by 1/∆a so that the new peak frequency is Finally, the universe transitions to and from MD at the critical frequencies From this starting point, the entire spectrum can be found by scaling as f 3 as long as the Universe is in RD and as f 1 as long as it is in MD. The final spectrum is then We see all of the effects that we had mentioned before. At high frequencies, there is a suppression of 1/∆a = (T r→m /T m→r ) −1/3 while at low frequencies there is an enhancement of (∆a) 3 ∼ T r→m /T m→r . A limit of the previous scenario is when the period of matter domination extends to the point in time when the GWs were produced. As shown in Fig. 5, in this case the power in modes that are sub-horizon when the spectrum is generated scales as k 3 while for those outside it goes as k. Thus even though this also has a transition from a k 3 scaling to a k scaling, the origin for the k 3 scaling is different between the two scenarios, sub-horizon physics instead of radiation domination. It is then an interesting question to see if it is possible to differentiate between the two. As the transition from sub-horizon to MD and the one from RD to MD both lead to identical scaling away from the transition region, we cannot differentiate between the two based on their spectrum alone. However, the two transitions are not identical. Thus, in principle, if one sees the transition between a k 3 scaling to a k 1 scaling, one could use its shape to determine whether the k 3 scaling was due to radiation domination or modes being sub-horizon.
We illustrate the difference between these two scenarios in Fig. 6. From this, one sees that the two different transitions are in fact distinguishable, but that their deviation is rather small. Without extra information from shorter wavelengths that are sensitive to the generation mechanism for the waves, one would need an extremely precise measurement of the transition region in order to determine which kind of transition it is.

Using the GW spectrum to measure w(t)
In the early universe, it is highly likely that w was not constant and was a changing function of time. In this section, we consider how to measure w(t). We will demonstrate that w(t) can be to a very good approximation read off directly from the slope of dΩ gw /d log k.
In an expanding universe with a constant w, the slope of dΩ gw /d log k for super-horizon  Figure 6. Comparison between the evolution of the tilt of GW spectra in two different scenarios: transition from sub-horizon modes to MD super-horizon modes (green) and from RD super -horizon modes to MD super-horizon modes (purple). Sub-horizon and RD super-horizon modes both scale as dΩ gw /d log k ∼ k 3 while MD super-horizon modes scale as dΩ gw /d log k ∼ k so that both scenarios have the exponent transition from 3 to 1. As is shown in the picture, the transition region is different, showing that the two scenarios can in principle be differentiated. modes is given by Eq. (2.24) and repeated below for convenience When one considers a generic w(τ ), there is no longer an exact solution. However, the intuition developed before still applies. As mentioned before, we can approximate the GW as having a suppressed production coupled with being frozen out until k = H. These combined to give Eq. (2.13), repeated below for the sake of convenience As long as w (τ ) H, then we can proceed as before and obtain where, for each momentum k, w(τ ) should be taken at the time τ when k re-enters the horizon, k = H(τ ). We see that in this approximation, one can simply read off w(τ ) straight from the slope of dΩ gw /d log k. More explicitly d log (dΩ gw /d log k) d log k = 1 + 15w(τ ) 1 + 3w(τ ) . Following this expression, one could in principle extract w(τ ) by just measuring the tilt of the low-frequency tail of the GW spectrum. The approximations made to obtain this simple expression are not valid during the transition region between matter and radiation domination. To see if they work in practice, we compared this procedure to several w(τ ) numerically. The results are shown in Fig. 7. We see that this procedure is surprisingly accurate.

Conclusion
In this paper, we have presented a physical picture of how to understand causality-limited gravitational waves. This picture allows one to easily estimate and predict the shape and behavior of low frequency gravity waves. As an example, we showed that generally sub-horizon and super-horizon modes behave very differently, allowing one to make model independent measurements of the conformal Hubble rate at which the gravity waves were created. Additionally, we showed how an intermediate period of MD has the effect of suppressing power at high frequencies and enhancing power at low frequencies.
Perhaps the most surprising result is that we showed that free-streaming particles change the scaling of the GW power spectrum. Unlike the usual case where particles start freestreaming after GWs are produced, if there are free-streaming particles present when GWs are produced, they lead to a suppression in the production of super-horizon gravitational waves, bringing the scaling from k 3 down to a maximum of k 4 . In the limit of a large number of free-streaming particles, oscillatory features are present demonstrating that interesting effects unrelated to causality can sculpt the spectrum of gravity waves.
There is much still left to be done with causality-limited gravitational waves. This range of the spectrum is special because its shape can be calculated from first principles. However, because causality sends the power to zero at low frequencies, a dedicated analysis would be needed to see to what extent gravity wave detectors such as LISA would be able to detect changes on this falling spectrum. On a separate note, since preheating and reheating generated GWs are typically followed by a period of matter domination, they are more visible than previously expected. It would be interesting to see if these models naturally generate signals detectable by existing or future gravitational wave detectors. where τ r is some fixed reference time at which the EOS was still determined by n k . Fixing the reference time τ r allows us to compare different modes that enter horizon during an era of fixed n.
We will now use the result in Eq. (A.6) to compare two scenarios to show the effects of non-standard expansion histories. As before in the first scenario, GWs are produced at a temperature T at a scale factor a and the universe is radiation dominated until the CMB. In the second scenario, GWs are produced at the same temperature as before, T , but there is an intermediate period of MD between the temperatures T r→m (where the universe becomes matter dominated) and T m→r (where the universe becomes radiation dominated). For simplicity we will take the scale factor at production to be the same a , which implies that the scale factor today will differ due to the difference in expansion history.
As discussed in Section 4.1, one of the major effects of an intermediate period of MD is that the frequency spectrum of GWs is shifted by the additional expansion due to the matter dominated phase. Since we normalized the scale factor at generation a to be the same in both scenarios, the relation between the scale factor and plasma temperature is identical in both cases up until the temperature reaches T r→m , when matter domination starts in the second scenario. The expansion history of both scenarios differ from that point on until the temperature gets to T m→r , after that point the scale factors corresponding to the same normalization at late times correspond to distinct frequencies in the two cases. In order to compare the spectra of the two scenarios and their relation to experimental searches, it is more useful to compute the differential spectrum in terms of frequency, f = k/(2πa): Using the fact that scale factors after the end of MD satisfy a m→r,1 /a 1 = a m→r,2 /a 2 , we find dρ 2 d log f (f ) = a m→r,2 a m→r,1 where f r→m = a r→m H r→m /(2πa 2 ) and f m→r = a m→r,2 H m→r /(2πa 2 ), which reproduces the result discussed in the main text.