Gravitational-wave cosmology across 29 decades in frequency

Quantum fluctuations of the gravitational field in the early Universe, amplified by inflation, produce a primordial gravitational-wave background across a broad frequency band. We derive constraints on the spectrum of this gravitational radiation, and hence on theories of the early Universe, by combining experiments that cover 29 orders of magnitude in frequency. These include Planck observations of cosmic microwave background temperature and polarization power spectra and lensing, together with baryon acoustic oscillations and big bang nucleosynthesis measurements, as well as new pulsar timing array and ground-based interferometer limits. While individual experiments constrain the gravitational-wave energy density in specific frequency bands, the combination of experiments allows us to constrain cosmological parameters, including the inflationary spectral index, $n_t$, and the tensor-to-scalar ratio, $r$. Results from individual experiments include the most stringent nanohertz limit of the primordial background to date from the Parkes Pulsar Timing Array, $\Omega_{\rm gw}(f)<2.3\times10^{-10}$. Observations of the cosmic microwave background alone limit the gravitational-wave spectral index at 95\% confidence to $n_t\lesssim5$ for a tensor-to-scalar ratio of $r = 0.11$. However, the combination of all the above experiments limits $n_t<0.36$. Future Advanced LIGO observations are expected to further constrain $n_t<0.34$ by 2020. When cosmic microwave background experiments detect a non-zero $r$, our results will imply even more stringent constraints on $n_t$ and hence theories of the early Universe.

Gravitational-wave astronomy is now a reality.The LIGO Scientific Collaboration has recently announced the first direct detection of gravitational waves coming from the merger of a binary black hole [1].Other experiments worldwide are ready to measure gravitational radiation across a wide range of frequencies.From the cosmic microwave background (CMB) to ground-based GW interferometers, these experiments cover more than 21 orders of magnitude in frequency-29 with complementary but indirect bounds from big bang nucleosynthesis (BBN), CMB temperature and polarization power spectra and lensing, and baryon acoustic oscillations (BAO) measurements.Each of these experiments is sensitive to a primordial stochastic GW background, originating from quantum fluctuations in the early Universe, and amplified by an inflationary phase [2][3][4][5].Standard inflationary models predict a primordial GW background whose amplitude is proportional to the energy scale of inflation [6].Observations of primordial GWs therefore provide unique insights into poorly understood processes arXiv:1511.05994v2[astro-ph.CO] 28 Feb 2016 in the very early Universe and its evolution from 10 −32 s after the Big Bang through to today.
In standard inflationary theories, the GW energy spectrum is expected to be nearly scale-invariant-above a certain frequency, the GW energy density decreases monotonically with increasing frequency [6].The gravitational field has quantum mechanical fluctuations which are dynamic at wavelengths smaller than the cosmological horizon, H −1 = 3c 2 /(8πGρ), and static due to causality at wavelengths larger than the horizon.During inflation, modes are redshifted and pulled outside the horizon where their power is frozen in with an amplitude that corresponds to the size of the cosmological horizon, and hence, to the energy density of the Universe at that time.As inflation progresses the energy scale of the Universe decreases, and the cosmological horizon grows.This is a consequence of the null energy condition, which posits that the energy density of the Universe cannot increase as a function of time.Modes that freeze out at larger physical wavelength have less power in them.Therefore, the slowly and monotonically decreasing energy density of the Universe during inflation is responsible for the monotonically decreasing shape of the primordial power spectrum of all fields.Spectra that decrease with increasing frequency are referred to as "red" spectra, and those that grow with increasing frequency are "blue".
A red spectrum, combined with observational constraints on the amplitude of GWs from the CMB, imply that GW detectors such as Pulsar Timing Arrays (PTAs) and ground-based interferometers such as the Laser Interferometer Gravitational-wave Observatory (LIGO) [7] and Virgo [8] are not sufficiently sensitive to detect primordial GWs predicted by the simplest model of inflation [e.g., 9].Detection at frequencies at or above PTAs may require extremely ambitious detectors such as the Big Bang Observer [10] or DECIGO [11].However, some non-standard models for the early Universe predict blue GW spectra, which could be detected by PTAs and/or LIGO (see below).
A blue spectrum can be generated from inflation depending on what happens when GW modes exit the horizon, either by non-standard evolution of the Universe during inflation or if there is non-standard power in these modes when they exit.This idea gained recent popularity in the wake of some early interpretations of the BICEP2 observations [12], where a flat GW spectrum was unable to simultaneously explain both the lower-frequency Planck observations [13] and the higher-frequency BI-CEP2 results [e.g., [14][15][16].
Standard models of inflation suggest that the slope of the GW spectrum should be approximately equal to the slope of the power spectrum of density perturbations.This prediction can be modified by having more than just a simple scalar field driving inflation.These non-minimal models can predict either red GW spectra whose spectral index varies from that of standard inflation [17] or blue spectra [e.g., 18,19].The latter modification is so dra-matic that the system violates the null-energy condition; a desirable, but by no means compulsory, property of the stress-energy tensor.Alternatively, blue spectra can be generated if the propagation speed of primordial GWs varies during inflation [20], or by introducing new interactions between the scalar field and gravity, where these interactions are low-energy remnants of some (unknown) modification of general relativity at much higher energy scales such as the Planck scale.Couplings of this form do not change any of the standard predictions of general relativity, but the theories that predict them allow us to treat the (unknown) high-energy theory of gravity in an effective low-energy limit for some energy scales.The simplest of such effective field theories produce a blue spectrum [21].
It is also possible to abandon inflation altogether and replace it with a scenario which preserves the observed spectrum of density perturbations.Two classic examples are string-gas [22] and ekpyrotic cosmologies [23].In the former, an ensemble of fundamental strings have thermodynamic properties that produce a high-temperature, quasi-static state, which produces a blue GW spectrum, whose size is comparable in magnitude to the standard red spectrum [24,25].Ekpyrosis posits that the primordial spectrum of perturbations is a result of a pre-big bang contracting phase.Such a phase has an increasing energy density and would create a blue power GW spectrum [23,26].
Importantly, a blue primordial GW spectrum may yield a primordial background immediately below present day limits, which may be detectable in the near future.While CMB experiments are likely to make direct measurements of the tensor-to-scalar ratio, r, they will poorly constrain the tensor index, n t .In this paper, we show how the combination of constraints on the primordial GW background from CMB, PTA, BBN, BAO and ground-based interferometer GW experiments can place stringent constraints on n t , yielding insights into the physics of the early Universe not accessible by any other means.

I. GRAVITATIONAL WAVE EXPERIMENTS
Current results from experiments trying to measure the primordial GW background do little to constrain the possible tilt of the spectrum.However, combining GW experiments over all frequencies allows us to constrain cosmological parameters from non-standard inflationary cosmologies [27][28][29][30].Combined CMB observations from the Planck satellite and the BICEP2 experiment constrain the stochastic GW background at frequencies of ∼ 10 −20 -10 −16 Hz, while PTAs are sensitive to GW frequencies of ∼ 10 −9 -10 −7 Hz and ground-based interferometers are sensitive at ∼ 10-10 3 Hz [31].As we show below, constraints on the total energy-density of GWs from BBN, gravitational lensing, CMB power spectra and BAO are sensitive to GWs as high as 10 9 Hz.There-fore, even a small blue tilt in the GW spectrum may be detectable in the GW frequency band covered from the CMB to LIGO/Virgo and provide more stringent constraints on the overall shape of the GW background [30].
A first step in our effort to apply experimental constraints to the GW energy-density spectrum is to assume that it can be well-approximated by a power law: where the pivot frequency, f CMB , is taken to be the standard value f CMB = (c/2π)0.05Mpc −1 [e.g., 32].It is conventional to re-express the amplitude of the primordial GW spectrum in terms of the tensor-to-scalar ratio, r ≡ A t /A s , where A s is the amplitude of the primordial power spectrum of density perturbations, and both are evaluated at the pivot scale.
Equation (1) is the simplest approximation one can make about the primordial GW spectrum.Most early-Universe theories predict only a small deviation from pure power-law behavior.The next level of complexity is to replace n t with n t +α t ln(f /f CMB )/2 in Eq. ( 1), where α t is known as the running of the spectral index.For example, single-field, slow-roll inflationary models predict α t (1 − n s ) 2 [e.g., see 33], where n s = 0.9645 ± 0.0049 is the measured value of the scalar spectral index [13].Therefore, within this class of theories, we expect a correction to the total GW power-law index of 10 −2 .Over 29 decades in frequency, this can have a marginal effect on the results presented here, a point we discuss in more detail below.
Due to the expansion of the Universe, the primordial GW spectrum that we observe today has evolved since it was created.This evolution is expressed in terms of a transfer function, T (f ), which encodes information about how GWs change as a function of frequency [34].The energy density of GWs today is given by The PTA/LIGO communities commonly present the GW spectrum in terms of the energy density in GWs as a fraction of the closure energy density per logarithmic frequency interval [31,35], where ρ c ≡ 3c 2 H 2 0 /(8πG), H 0 = 100h km s −1 Mpc −1 is the Hubble expansion rate, and h = 0.67 is the dimensionless Hubble parameter [13].Indirect constraints on the GW background are typically "integral bounds" on Assuming a standard expansion history that includes non-relativistic matter and radiation, the GW spectrum today is given by [34,[36][37][38]: where f eq is the frequency of the mode whose corresponding wavelength is equal to the size of the Universe at the time of matter-radiation equality with frequency Here, Ω m and Ω r are respectively the total matter and radiation energy-density evaluated today, and Equation ( 4) is key in our analysis since it allows us to combine constraints on r and n t from the CMB with constraints to Ω gw (f ) and n t from PTAs and LIGO.We use cosmological parameters obtained from the latest Planck satellite data release [13].
In the following sections, we combine observational constraints spanning 29 orders of magnitude in frequency to derive stringent constraints on backgrounds with a non-zero spectral index n t .Figure 1 highlights the key idea: we show the current best upper limits on Ω gw (f ), and a series of curves given by Eq. ( 4) that are constrained by these limits.Starting from the lowest frequency limits, we summarize current upper limits before combining them to derive joint constraints.
Note that in evaluating the observed GW spectrum we have neglected the effects of neutrino free-streaming [39] and phase transitions occurring in the early universe [36].As shown in Refs.[36,39,40], the free-streaming of neutrinos and phase transitions during the very early Universe lead to a suppression of Ω gw by a factor of ∼ 1/2 -1/3 for PTA and LIGO frequencies.However, because of the large lever-arm between the frequencies probed by these detection methods and the CMB, including this suppression in the analysis changes our constraints on n t by only a few percent.We therefore neglect these effects in our analysis.
Three recent papers have presented combined constraints on n t and r using some combination of CMB, LIGO and PTA data [30,40,41] 1 .Huang and Wang [41] presented their analysis soon after the original BICEP2 results were reported [43,44], and as such, focussed on the fact that those data preferred a slightly positive blue tilt for the tensor power spectrum, which resulted from the inconsistency between the original BICEP2 data and Planck observations.On the other hand, Liu et al. [40] presented constraints from only CMB and PTA data, but focussed on what a positive detection could do for our understanding of the early-Universe equation of state, cosmic phase transitions and relativistic free-streaming.
Our analysis improves on those of Refs.[30,40,41] in a number of significant ways.First, we include the indirect GW constraints in a self-consistent way, which allows us to compare integral and non-integral constraints with varying spectral indices (see Section I D).Second, we present a new analysis of PPTA data [45] that give the best limit on Ω gw (f ) in the PTA band by a factor of four over previous published results.Finally, we provide our own analysis of the raw PPTA time-of-arrival data to allow for varying spectral indices, instead of assuming a constant n t for PTA observations taken from older PPTA analyses as is done in [30,40].

A. CMB Intensity and Polarization
Primordial GWs imprint a characteristic signal onto the intensity and polarization of the CMB that can be measured by ground-based and space-borne observatories.A joint analysis [13,44] of Planck satellite and BICEP2/Keck array data found that r < 0.12 at 95% confidence level (CL) under the assumption that n t = 0.The solid black curve in Fig. 1 labeled "CMB" shows the estimated sensitivity of the Planck satellite.The CMB sensitivity curve is calculated by determining the value of the spectral density Ω gw (f ) that yields a marginally detectable signal given a model of the Planck satellite noise properties [46].
Observations of the CMB intensity and polarization are analyzed by re-expressing the real-space data in a spherical harmonic expansion.The intensity measurements can be expanded in the standard (scalar) spherical harmonics, whereas the polarization data must be expanded in spin-weighted spherical harmonics [47].We can further separate various physical processes by dividing the polarization data into a curl-free (E-mode) and curl (B-mode) basis.In order to compare these data to a theoretical model, the measured spherical harmonic coefficients are further analyzed to estimate their statistical correlations.The presence of a primordial GW spectrum fundamentally alters the expected correlations leading to an enhanced correlation for the intensity of the CMB on the largest angular scales as well as a non-zero correlation for the B-mode polarization [48,49].
The expected effect of a non-zero primordial GW spectrum on the CMB is calculated by solving the Boltzmann equation for the various components of matter that fill the Universe.The Boltzmann equation for the photons encodes all of the information about correlations in the intensity and polarization of the CMB.In particular, for each spherical harmonic multipole, the expected correlations can be expressed as an integral over cosmic time and frequency [50].Therefore, the total expected CMB signal due to primordial GWs can be expressed as an integral over its spectrum, Ω gw (f ).
The CMB sensitivity curve shown in Fig. 1 is calculated by setting the total CMB signal-to-noise ratio equal to two (corresponding to a 95% CL bound).The squared signal-to-noise ratio is calculated from a sum in quadrature of CMB B-mode multipoles divided by the estimated polarization noise for each multipole from the Planck satellite's 143 GHz detector [46,51].Since, as dis-cussed above, each B-mode multipole is an integral over the GW spectrum, we express the integral over frequency as a sum so that we can evaluate the contribution at each frequency interval.We relate the primordial amplitude to the present-day spectral density using Eq. ( 21) of [36].The noise is calculated under the hypothesis of no primordial GWs, although weak gravitational lensing also induces a non-zero B-mode correlation, which we treat as an additional source of noise.Finally, the limit is converted into a "power-law integrated" curve using the formalism from Thrane and Romano [52].Any model intersecting this curve is ruled out at 95% CL.
Current CMB constraints to r and n t come from measurements made by Planck [53], the BICEP2/Keck array [44], and SPTpol [54].Constraints from these datasets were determined using the Boltzmann solver CAMB and a modified version of the Monte Carlo stepper cosmoMC [55][56][57].

B. Pulsar Timing Arrays
The incoherent superposition of primordial GWs is expected to imprint on the arrival time of pulses from the most stable millisecond pulsars.A number of PTAs around the world are engaged in the hunt for GWs, including the Parkes Pulsar Timing Array [PPTA; 58], the North American Nanohertz Observatory for Gravitational Waves [NANOGrav; 59], and the European Pulsar Timing Array [EPTA ; 60].Here, we use recent data from the PPTA [45] to provide the strongest constraints to date on Ω gw (f ) from a primordial background in the PTA band.
The PPTA monitors 24 pulsars with the 64-m Parkes radio telescope in a bid to directly detect GWs, and currently has the most stringent upper limits on the GW background from supermassive black hole binaries [45].We derive our limit on the primordial GW background by performing a similar Bayesian analysis to that in [45], with the exception that we utilize the Bayesian pulsar timing data analysis suite PAL22 , and allow for an arbitrary strain spectral index.
The GW spectrum in the PTA band can be approximated as a power law, with where A gw is the amplitude of the characteristic strain at a reference frequency of f yr ≡ yr −1 .The star in Fig. 1, labelled "PTA", is the 95% CL upper limit assuming a spectral index of n t = 0. temperature (GeV) FIG. 1. Experimental constraints on Ωgw(f ); the black star is the current PPTA upper limit and all black curves and data points are current 95% confidence upper limits.The grey curve and triangle are respectively the predicted aLIGO sensitivity and PPTA sensitivity with five more years of data.The indirect GW limits are from CMB temperature and polarization power spectra, lensing, BAOs, and BBN.Models predicting a power-law spectrum that intersect with an observational constraint are ruled out at > 95% confidence.We show five predictions for the GW background, each with r = 0.11, and with nt = 0.68 (orange curve), nt = 0.54 (blue), nt = 0.36 (red), nt = 0.34 (magenta), and the consistency relation, nt = −r/8 (green), corresponding to minimal inflation.
the PPTA limit are the upper limits from the EPTA [61] and NANOGrav [62].Both the EPTA and NANOGrav present limits on the GW energy density from inflationary relics assuming n t = 0; our new limit for n t = 0 (cf.our limit on Ω gw (f ) for n t = 0.5 which only differs in the second decimal place) is a factor of 4.1 better than the previous best limit from [62].
The grey triangle below the star in Fig. 1 is a predicted GW upper limit derived by simulating an additional five years of PPTA data.We took the maximum likelihood red noise parameters in the existing data sets, estimated the white noise level using the most recent data that represents current observation quality, and assumed a twoweek observing cadence to derive the 95% CL upper limit of Ω 95% gw (f ) 5 × 10 −11 .However, the PPTA limit will be superseded before 2020 with limits placed from collating datasets from the three existing PTAs as part of the International Pulsar Timing Array [IPTA; 63].From Fig. 1, it becomes clear that PTAs may not play a significant role in constraining inflationary models where the GW spectrum is described by Eq. (3) when aLIGO reaches design sensitivity, given the significant improvements in the latter experiment.However, PTAs can still play an important role for cosmological models with a varying spectral index; that is, with a non-negligible running of the spectral index α t .
Giblin and Thrane [64] recently proposed a "rule of thumb" for the maximum GW energy density for cos-mological backgrounds based only on arguments of the energy budget of the Universe at early times.They presented optimistic, realistic and pessimistic upper limits for Ω gw (f ), with the optimistic limit representing the largest value of Ω gw (f ) possible given a reasonable set of conditions.The new PPTA limit reported here is the first time a GW limit in either the PTA or LIGO band has gone under this optimistic threshold, thus marking the first time the detection of cosmological GWs could actually have been possible according to arguments in [64].Conventional models of early-Universe particle physics do not predict such a large GW background in the PTA frequency band.The temperature of the Universe at the time when such GWs are produced is ∼ 1 GeV (see top axis of Fig. 1), a temperature at which physics of the early Universe is relatively well known.We note that the possibility of first-order phase transitions that generate a strong GW background in the PTA frequency band is not completely ruled out [e.g., [65][66][67][68] Of course it is possible that there is unknown physics that influences gravity without coupling strongly to the standard model of particle physics that could produce a strong GW background in the PTA frequency band.Combined, two-dimensional posterior distribution for the tensor-to-scalar ratio, r, and the blue tilt of the GW spectrum, nt, using CMB, PPTA, indirect, and LIGO observations.The contours are the 95 and 99% limits.The green, dashed curve shows the consistency relation, nt = −r/8, while the red and blue triangles correspond respectively to the red and blue curves in Fig 1 .For clarity, the right panel is a zoomed-in version of the left panel, with additional posterior distributions shown.See Section II A for a description of each posterior distribution.

C. Ground-based interferometers
LIGO [7] and Virgo [8] are long-baseline, ground-based GW interferometers with best sensitivity at frequencies of 10 2 -10 3 Hz.Data collected during the initial phases of these instruments have been used to place upper limits on a stochastic background of GWs from astrophysical and cosmological sources [9,69,70].We utilize data from the initial LIGO and Virgo observatories.These data were collected in 2009 -2010 as part of the fifth LIGO Science run.Two limits were originally obtained using these combined observations: a lower-frequency limit from combined LIGO-Virgo observations that assumed a flat, i.e., n t = 0, spectrum [69], and a higher-frequency limit from an analysis of the two co-located LIGO detectors at Hanford, which assumed n t = 3 [70].
We implement a new way to analyze LIGO/Virgo limits on the primordial background that allows for a varying spectral index.The analysis goes beyond Meerburg et al. [30], and Huang and Wang [41], which both assume n t = 0 for their LIGO/Virgo constraints.We combine published data from Refs.[69] and [70] to generate a power-law integrated curve [52], shown in Fig. 1.Any power-law model intersecting a power-law integrated curve is ruled out at 95% CL.Then, utilising the formalism from Mandic et al. [71], we obtain constraints on Ω gw (f ) for arbitrary spectral indices.The limits on Ω gw (f ) are converted into constraints on n t and r.
At the time of writing, the LIGO experiment has begun taking data for the first observing run of the advanced detector era, with the Virgo experiment to follow in 2016.At design sensitivity, advanced detectors are forecast to achieve nearly four orders of magnitude of im-provement in Ω gw (f ); see the curve marked "aLIGO" in Fig. 1, which is the projected sensitivity given two LIGO detectors operating for one year at design sensitivity.

D. Indirect constraints
Indirect constraints on GW backgrounds have been obtained using a variety of data including CMB temperature and polarization power spectra, lensing, BAOs, and BBN [e.g., [72][73][74].Indirect bounds are "integral bounds," which apply to Ω gw and not to Ω gw (f ); see Eq. 3. Recently, Pagano et al. [75] combined the latest Planck observations of CMB temperature and polarization power spectra and lensing with BAO and BBN measurements (specifically observations of the primordial Deuterium abundance) to put an integral constraint on the primordial GW background of Ω gw < 3.8 × 10 −6 .
While there is a long history in the literature of plotting Ω gw integral bounds alongside Ω gw (f ), they are not directly comparable.However, the two quantities can be related if we assume that Ω gw (f ) is described by a powerlaw spectrum with a known cut-off frequency, which we choose to be f max = 1 GHz, corresponding to an energy scale typical of inflation, T = 10 17 GeV.Given this plausible assumption, we plot the indirect constraints in Fig. 1 as power-law integrated curves using the formalism from Ref. [52].Any power-law model intersecting a power-law integrated curve is ruled out at 95% CL.
Inspecting Fig. 1, it is apparent that the current best constraints on n t come from observations of the CMB combined with indirect bounds.
The strength of the indirect bounds depends in part on our choice of f max = 1 GHz, however changing the cut-off frequency by several orders of magnitude would not the change qualitative picture.For example, in alternative theories of inflation it is possible to posit an energy scale as low as 10 6 GeV, corresponding to a cut-off frequency of f max = 10 mHz.This choice of cut-off frequency shifts the minimum of the indirect bound curve from ∼ 100 mHz to ∼ 1 µHz, while the minimum value of Ω gw (f ) increases by a factor of ∼ 2. When aLIGO reaches design sensitivity, it will surpass indirect constraints on primordial backgrounds with non-running spectral indices.

II. COMBINED CONSTRAINTS ON THE PRIMORDIAL TILT A. Combined experimental constraints
Here we combine the current limits on Ω gw (f ) from the individual experiments mentioned above to constrain the tensor-to-scalar ratio, r, and the tensor index, n t .In Fig. 2, we plot these two-dimensional posterior distributions for r and n t .In both panels we plot two theory points and a theory curve.The green, dashed curve corresponds to the consistency relation from standard inflationary models (n t = −r/8), and the red and blue triangles have the same values of n t and r as the red and blue curves in Fig. 1.
Figure 2 combines the constraints from each experiment in a heuristic manner.The left panel shows all of the CMB constraints from direct detection experiments starting with Planck (grey shaded region), adding BI-CEP2 (green), and finally adding SPTpol (red) to get the overall constraints on the CMB from direct GW observations.Also plotted in the left panel is the PPTA posterior (blue).The PPTA search algorithm described in Section I B derives a posterior distribution in terms of n t and Ω gw (f ), which is converted to n t and r using Eqs.( 4) and (5).Finally, in the left-hand panel of Fig. 2, we plot the combined CMB + PPTA posterior (black).This distribution represents the state-of-the-art constraints one can derive from CMB and PTA experiments alone.At a reference value of r = 0.11, this limits n t < 0.68 with 95% confidence.In general, the 95% CL upper limit on n t as a function of r derived from these constraints is well-approximated by the simple relation where A = −0.13 and B = 0.68.Equation ( 7) allows one to extrapolate the constraints on n t to arbitrary small values of r.With the addition of each experimental constraint we get tighter limits on A and B; the 95% CL upper limits for each experiment are collated in Table I.
In the right panel of Fig. 2 we retain the CMB (red, labelled 'Planck + BICEP2 + SPTpol') and CMB + PPTA (black) posterior distributions from the left-hand panel.As with the PPTA analysis, the LIGO/Virgo data analysis algorithm described in Section I C constrains n t and Ω gw , which we convert to n t and r using Eqs.( 4) and ( 5).This distribution is shown in yellow, and the combined CMB + PPTA + LIGO/Virgo constraints are shown in pink.This pink contour represents the current best constraint from the direct GW experiments covering 21 decades in frequency.At a reference value of r = 0.11, the CMB + PPTA + LIGO/Virgo constraints yields an upper limit of n t < 0.54.For smaller values of r, the theoretical curves at the CMB frequency take on lower values of Ω gw , which implies higher-frequency experiments play even more of a role in constraining n t than comparatively lower-frequency experiments.For example, at r = 0.01, CMB + PPTA constraints imply n t < 0.82, while CMB + LIGO/Virgo constraints imply n t <0.60.The 95% CL from the CMB + LIGO constraint is well approximated by Eq. ( 7) where values for A and B can be found in Table I.
Also in the right panel of Fig. 2, we add the indirect constraints described in Section I D to the direct constraints (turquoise).This contour represents our total knowledge of n t and r using all experimental constraints.In this case, at r = 0.11 we find n t < 0.36 at 95% CL, and the upper limit as a function of r is well approximated by Eq. (7), where the values for A and B can be found in Table I.
From Fig. 2, it is clear that only direct observations of the CMB constrain n t < 0. This makes sense in the context of Fig. 1: given that the lever-arm for the GW theory curves, i.e., Eq. ( 4), are hinged at f CMB , a negative spectral index is not constrained by experiments that are only sensitive to values of Ω gw (f ) higher than at the CMB.
Finally, in Fig. 2 we show the projected constraints that one can expect by the year 2020 (dark blue, labelled ". . .+ aLIGO + PPTA(2020)") assuming five more years of PPTA observations and aLIGO at design sensitivity (see Sections I B and I C).These contours show that the constraint for the spectral index improves to n t 0.34 at r = 0.11, and is well-approximated by Eq. ( 7) with A and B given in Table I.As is evident from Fig. 1, the constraint at high n t will be dominated by aLIGO.Similar constraints in the PTA band are not expected until the era of the Square Kilometre Array and Five hundred meter Aperture Spherical Telescope (FAST) [e.g., 40,76].7).The value of B is therefore the 95% upper limit of nt at a reference value of r = 0.11.

B. Comparison with theory
In the previous subsection, we present stringent constraints on the blue tilt of the primordial GW background from experiments spanning 29 decades in frequency.These results can be used to comment on early Universe models.Those models whose spectral indices are near zero -or of comparable magnitude to standard inflationary models -are consistent with the data.String-gas cosmologies and modified inflationary scenarios with non-minimal couplings to gravity seem to be the least constrained, since these models predict relatively small values of n t and are unconstrained for even relatively large tensor-to-scalar ratios.
Ekpyrosis has a tendency to predict large values of the spectral tilt including n t ≈ 2 [23] that come from modes freezing out of the horizon during the contracting phase of the Universe.More modern incarnations of ekpyrosis produce blue tilts, but with relatively low values of tensor-to-scalar ratio, r, [e.g., 26].Here, our derivation of fitting formulae for n t as a function of r (see Table I) allow specific ekpyrotic predictions to be tested to arbitrary small values of r.
Our results will have important implications following the detection of non-zero tensor-to-scalar ratio by a future CMB experiment (e.g., see Ref. [77] 3 and references therein).Such a detection, together with the data from PTAs and ground-based interferometers will put very tight limits on n t , with larger values of r being the most constraining.For example, a confirmed detection of r ≈ 0.1 would put very tight bounds on n t with a strong preference for small and positive values.Such tight constraints are truly a result of CMB bounds on the low frequency end and PTA, LIGO and indirect bounds on the upper end.
When a detection is made in any of the frequency bands studied herein, it becomes even more pertinent to analyze all experimental data, consistently taking into account the spectral running, α t .Indeed, in the case of a detection, upper limits in each frequency band can be used to simultaneously constrain n t and α t ; three or more experiments are required to constrain both parameters.In a future work, we will present three-dimensional posterior constraints that include r, n t and α t , and also incorporate predictions for future CMB experiments.Distinguishing primordial backgrounds from astrophysical foregrounds may be a daunting task, though multiwavelength measurements could prove useful toward this end.

III. CONCLUSION
By combining limits from many different GW experiments probing 29 decades in frequency, we present new constraints on cosmological parameters n t and r, which are intimately related to the evolution of the early Universe.This interdisciplinary research also makes significant advances in PTA, LIGO/VIRGO and indirect GW limit analysis techniques.Specifically, we present new PPTA data that provides the most stringent limit on the primordial gravitational-wave background, Ω gw (f ) < 2.3 × 10 −10 ; more than a factor of four tighter than the previous best limit from [62].Moreover, we develop and implement a method to give the best limits on the primordial background from ground-based interferometers-a method we anticipate will become standard in future LIGO/Virgo primordial background analyses.Furthermore, we provide a new interpretation of indirect GW constraints from CMB temperature and polarization measurements, lensing, BBN and BAO observations that allow for a varying primordial spectral index, allowing us to directly compare these "integral" constraints on Ω gw with the usual frequency-dependent Ω gw (f ) constraints.Our technique for comparing direct and indirect limits can be widely adopted within the GW community to avoid the confusion created from 'applesto-oranges' comparisons.
While Refs.[30,40,41] present constraints on n t and r using combinations of CMB, LIGO and PTA data, the focus of their work was significantly different.Indeed, the work of Meerburg et al. [30] and Huang and Wang [41] were originally in response to the now defunct BICEP2 results [43,44], while Liu et al. [40] presented constraints from only CMB and PTA data, but focussed on what a positive detection could do for our understanding of the early-Universe equation of state, cosmic phase transitions and relativistic free-streaming.
A direct comparison between our results and that of Meerburg et al. [30] is not possible for a number of reasons.Notably, they use a linear prior on r, which, together with the use of the original BICEP2 results, ends in constraints that are not bounded below.Figure 1, together with Eqs. ( 4) and ( 5), show that r should be unbounded below given that there are no lower limits on the amplitude of Ω gw (f ) from any experiments.Moreover, Meerburg et al. [30] use an unconventional pivot scale for the theoretical GW spectrum (see, e.g., Ref. [32] for a discussion of the optimal pivot scale).
Our results are significantly more constraining than those of Liu et al. [40] (see their Fig.7), most notably due to the inclusion of indirect GW constraints.Our analysis quantifies how large the spectral index of the primordial spectrum, n t , can grow as a function of the tensor-to-scalar ratio, r; see Fig. 2 and Table I for a summary of the results.
Various theories of the early Universe predict a blue primordial gravitational-wave spectrum [e.g., 19,21,23,26], and indeed, some versions of ekpyrosis predict large values of n t which we can now rule out by our analysis -see Section II B. Observations of the CMB alone only limit the inflationary GW spectrum to n t 5 at a reference value of the tensor-to-scalar ratio of r = 0.11.Current observations by the PPTA and initial LIGO/Virgo reduce this limit to n t < 0.54 with 95% confidence, and including limits from indirect GW observations reduces this to n t < 0.36.We predict that observations by aLIGO at design sensitivity (circa 2020) will reduce this constraint to n t < 0.34.All upper limits on n t are applicable at a reference value of r = 0.11, but can be extrapolated to other values of r using Eqn.(7) and Table I.
Of course, it is a future direct detection of r that will have the most important implications.Such a detection will allow us to slice through the parameter space presented in Fig. 2, providing significant constraints on parameters governing theories of the early Universe.
FIG. 2.Combined, two-dimensional posterior distribution for the tensor-to-scalar ratio, r, and the blue tilt of the GW spectrum, nt, using CMB, PPTA, indirect, and LIGO observations.The contours are the 95 and 99% limits.The green, dashed curve shows the consistency relation, nt = −r/8, while the red and blue triangles correspond respectively to the red and blue curves in Fig 1.For clarity, the right panel is a zoomed-in version of the left panel, with additional posterior distributions shown.See Section II A for a description of each posterior distribution.

TABLE I .
95% CL upper limits on A and B as in Eq. (