Disentangling the primordial nature of stochastic gravitational wave backgrounds with CMB spectral distortions

The recent detection of a stochastic gravitational wave background (SGWB) at nanohertz frequencies by pulsar timing arrays (PTAs) has sparked a flurry of interest. Beyond the standard interpretation that the progenitor is a network of supermassive black hole binaries, many exotic models have also been proposed, some of which can potentially offer a better fit to the data. We explore how the various connections between gravitational waves and CMB spectral distortions can be leveraged to help determine whether a SGWB was generated primordially or astrophysically. To this end, we present updated $k$-space window functions which can be used for distortion parameter estimation on enhancements to the primordial scalar power spectrum. These same enhancements can also source gravitational waves (GWs) directly at second order in perturbation theory, so-called scalar-induced GWs (SIGWs), and indirectly through the formation of primordial black holes (PBHs). We perform a mapping of scalar power spectrum constraints into limits on the GW parameter space of SIGWs for $\delta$-function features. We highlight that broader features in the scalar spectrum can explain the PTA results while simultaneously producing a spectral distortion (SD) within reach of future experiments. We additionally update PBH constraints from $\mu$- and $y$-type spectral distortions. Refined treatments of the distortion window functions widen existing SD constraints, and we find that a future CMB spectrometer could play a pivotal role in unraveling the origin of GWs imprinted at or below CMB anisotropy scales.


INTRODUCTION
The recent detection of a stochastic gravitational wave background (SGWB) by numerous pulsar timing array collaborations (Agazie et al. 2023c;Antoniadis et al. 2023a;Reardon et al. 2023;Xu et al. 2023) at nanohertz frequencies represents a major step forward in our understanding of the underlying tensor perturbations that permeate our universe.Unlike the initial detection of gravitational waves by LIGO (Abbott et al. 2016), the progenitor of the SGWB detected in pulsar timing residuals is highly uncertain.The standard astrophysical interpretation is that such a background may be sourced by a population of super-massive black hole binaries (SMBHBs) residing at the centres of galaxies (Agazie et al. 2023a).The current paradigm of structure formation indicates that the growth of galaxies is mediated both via accretion as well as a series of hierarchical mergers (White & Rees 1978).It is well known that most galaxies host a super-massive black hole at their centre (Kormendy & Ho 2013;Akiyama et al. 2019), and in the process of merging, dynamical friction will drag the supermassive black holes towards the centre of the new galaxy.This process is known to produce SMBHBs (Begelman et al. 1980), where the gradual inspiral of the two black holes will emit gravitational waves in the nanohertz frequency regime.
The expected population of SMBHBs is uncertain, rendering the amplitude of the resultant SGWB difficult to estimate a priori.However, the timing residual power spectral density is expected to follow a power law with spectral index  BHB = 13/3 ≈ 4.33 (Phinney 2001).In the NANOGrav analysis (Agazie et al. 2023c), it was shown that a general power law fit to the observed background produces maximum likelihood values of  GWB ≃ 6.4 × 10 −15 , and  ≃ 3.2.This result is roughly 3 away from the theory prediction, introducing a mild tension between the data and this model.It should be noted that if one models dark matter fluctuations as a Fourier-domain Gaussian process (as opposed to the piecewise-constant representation consid-ered primarily by NANOGrav) these posteriors can shift significantly.Additionally, a SGWB sourced by SMBHBs is expected to exhibit some small level of anisotropy, but so far none has been detected, leaving the SMBHB interpretation on uncertain ground (Agazie et al. 2023b).
Following the announcement of this discovery, a number of papers have appeared which discuss the constraints and interpretations of this SGWB in ways which modify the small-scale ( ≳ 1 Mpc −1 ) primordial power spectrum.Two well-studied examples are by sourcing the signal through SIGW (Cai et al. 2023;Huang et al. 2023;Wang et al. 2023;Zhu et al. 2023;Yi et al. 2023a,b;Firouzjahi & Talebian 2023;You et al. 2023;Balaji et al. 2023;Yuan et al. 2023;Choudhury et al. 2023;Unal et al. 2023;Jin et al. 2023) or primordial black holes (Franciolini et al. 2023;Depta et al. 2023;Inomata et al. 2023;Wang et al. 2023;Bhaumik et al. 2023;Mansoori et al. 2023;Huang et al. 2023).Enhancements to the power spectrum on these scales can be constrained strongly by distortions to the frequency spectrum of the cosmic microwave background (CMB) (e.g., Chluba et al. 2012b).Dissipation of small-scale curvature perturbations (Silk damping) mix photons with different temperatures, inducing spectral distortions even if the power spectrum remains nearly scale invariant at arbitrarily high -modes.Naturally, if the scalar spectrum is enhanced on these small scales, significant distortions can be generated which allows us to derive constraints on a variety of models.
Crucially, for Gaussian scalar perturbations, COBE/FIRAS has already ruled out a primordial origin of the supermassive BHs residing in the centers of galaxies (Kohri et al. 2014;Nakama et al. 2018), which can be directly deduced from the limits first presented in Chluba et al. (2012b).The dissipation of small-scale perturbations during Big Bang Nucleosynthesis (BBN) further limits the amplitude of perturbations (Jeong et al. 2014).Various observational consequences of primordial black holes (PBHs) also present much weaker constraints at even smaller scales.
While it will undoubtedly take some time before a precise model stands above the rest, one major question can still to be addressed: are these gravitational waves of a primordial origin?One possible way to tackle this question is by looking for -mode polarization in the CMB.The LiteBIRD collaboration will eventually be able to directly constrain very large scale tensor modes (10 −19 Hz ≲  ≲ 10 −15 Hz) by searching for this polarization signal (Allys et al. 2023) from space.Additional ground based efforts such as the Simons Observatory (Ade et al. 2019) and CMB-S4 (Abazajian et al. 2019;Abazajian et al. 2022) will also probe -mode polarization signals, promising stringent limits in the coming ≃ 5 − 10 years.
There is, however, a complementary way to determine if a given tensor spectrum is primordial, using CMB spectral distortions.Similar to the scalar perturbations, distortions are also generated through the direct dissipation of tensor modes present at very early times (5 × 10 4 ≲  ≲ 10 6 ) (Ota et al. 2014;Chluba et al. 2015b;Kite et al. 2021a).In addition, models which produce the SGWB via small scale features in the power spectrum may also be probed by spectral distortions if these enhancements take place at 1 Mpc −1 ≲  ≲ 10 4 Mpc −1 (e.g., Chluba et al. 2012a,b;Khatri et al. 2012b;Pajer & Zaldarriaga 2013;Chluba & Grin 2013;Chluba et al. 2015a).
Within ΛCDM, a clear target of  = 2 × 10 −8 for the -distortion exists (Chluba et al. 2012a;Cabass et al. 2016;Chluba 2016), and any departure from this value would point towards new physics occurring in the pre-recombination era.This would provide significant evidence for not only a primordial source, but could also single out a small class of preferred exotic physics explanations.In this work, we further explore the synergies between GWs and CMB spectral distortions, highlighting how to properly leverage this invaluable information in a multi-messenger approach to unraveling the origin of the gravitational wave background.
In Section 2, we review the physical mechanism of spectral distortion generation from enhanced small scale power, as well as the dissipation of primordial tensor modes.We extend results in the literature by considering -type distortions in addition to the -distortions previously studied.We also provide a set of window functions for use in future computations, and showcase how broad features in the primordial power spectrum (PPS) generally strengthen constraints derived from CMB distortions.Following this, in Sections 3 and 4 we discuss the possible distortion signatures of scalar induced gravitational waves and primordial black hole models, which may also produce a significant SGWB.We conclude in Section 5.In what follows, we mainly take the NANOGrav results as a case study, keeping in mind that other PTA collaborations have reported qualitatively similar results.

SYNERGIES WITH SPECTRAL DISTORTIONS
Distortions in the frequency spectrum of the cosmic microwave background provide a sensitive probe of perturbations to the thermal history of the universe up to redshifts of roughly  th ≃ 2 × 10 6 .The inefficiency of number changing processes (Bremsstrahlung and double Compton scattering) below this redshift ensures that a blackbody spectrum cannot be recovered if departures from a thermal spectrum are present.Instead, Compton scattering drives the spectrum to a Bose-Einstein shape, characterized by a small chemical potential () at redshifts of   ≲  ≲  th (where   ≃ 5 × 10 4 ).At still lower redshifts ( ≲   ), Compton scattering freezes out and any nonthermal energy injection is instead driven towards a -type distortion.As the freeze out of Compton scattering is not instantaneous, energy release between 10 4 ≲  ≲ 3 × 10 5 will also generate residual type distortions.These distortions are not simply a superposition of  and  types, but also contain non-linear scattering residuals which can be used to probe the exact time dependence of energy injection scenarios (e.g., Chluba & Sunyaev 2012;Khatri & Sunyaev 2012;Chluba 2013a) and can be parameterized using the principal components of the distortion (Chluba & Jeong 2014).
One such mechanism for sourcing spectral distortions is through the dissipation of primordial perturbations in the pre-recombination era (Sunyaev & Zeldovich 1970b;Daly 1991;Hu et al. 1994;Chluba et al. 2012a).When present, energy stored within these perturbations is transferred to the thermal bath through electron scattering and free-streaming effects which mix blackbodies with slightly different temperatures.Assuming that the spectrum of scalar fluctuations remains nearly scale-invariant, the dissipation of power from small scale acoustic modes (50 Mpc −1 ≲  ≲ 10 4 Mpc −1 ) induces a distortion with amplitude ≃ 2 × 10 −8 .This is one of the most important primordial distortion signals expected within the standard concordance model of cosmology (ΛCDM) (Chluba 2016), and presents a tantalizing target in reach of next-generation spectral distortion experiments (Kogut et al. 2011(Kogut et al. , 2016;;Chluba et al. 2021), although mitigating low-frequency foreground contamination is a significant challenge (Abitbol et al. 2017).Limits on the power spectrum at these scales are complementary to the tight constraints one can derive from CMB temperature anisotropies (Planck Collaboration et al. 2020) at scales of 10 In addition to this standard model signal, many extensions exist which invoke an enhancement of power on sufficiently small scales.In particular, several inflationary scenarios can achieve small scale enhancement by postulating interactions between multiple fields (Silk & Turner 1987;Polarski & Starobinsky 1992;Braglia et al. 2020), or from adjustments to the inflaton potential (Ballesteros & Taoso 2018;Garcia-Bellido & Ruiz Morales 2017).Depending on the amplitude and wavenumber of these enhancements, it may be possible to form primordial black holes (PBHs) (Carr et al. 2021;Green & Kavanagh 2021), source gravitational waves (at second order in perturbation theory, see Domènech (2021) for a review), and induce significant spectral distortions (Chluba et al. 2012b;Nakama et al. 2018).The detection of these spectral distortions could provide an important multi-messenger signature of a stochastic gravitational wave background present at high redshifts, providing us with a pathway to disentangle astrophysical sources from primordial ones.
At present, the state of the art measurement of CMB spectral distortions comes from the COBE/FIRAS satellite, which reported limits on the distortion parameters of  ≤ 9 × 10 −5 and  ≤ 1.5 × 10 −5 (Fixsen et al. 1996) at 2.While these bounds can be strengthened by roughly a factor of 2 (Gervasi et al. 2008;Bianchini & Fabbian 2022), the next generation of ground based (TMS Rubiño Martín et al. 2020) and balloon borne (BISOU Maffei et al. 2021) experiments will likely further strengthen the bounds on / by one order of magnitude.Satellite efforts such as PIXIE (Kogut et al. 2011(Kogut et al. , 2016) ) or the Voyage 2050 program (Chluba et al. 2021), however, can increase the constraining power of distortions by several orders of magnitude.With current technology, a PIXIE-type mission is expected to reach a sensitivity of  ≲ 10 −8 and  ≲ 2 × 10 −8 , even if significant uncertainties with respect to foreground removal remain (Sathyanarayana Rao et al. 2015;Desjacques et al. 2015;Mashian et al. 2016;Abitbol et al. 2017;Rotti & Chluba 2021;Zelko & Finkbeiner 2021).Spectral distortions are therefore uniquely situated to provide us with a deep view behind the surface of last scattering due to the impressive increase in sensitivity possible with future experiments.

Spectral distortions from enhanced small-scale power
The dissipation of small scale modes is primarily caused by the diffusion of photons which gives rise to the well known damping tail of CMB anisotropies (Silk 1968;Planck Collaboration et al. 2020) at ℓ ≳ 500.In the context of spectral distortions, photon diffusion is equivalent to the mixing of blackbodies with slightly different temperatures which subsequently induces a distortion.This effect has been known for some time (Sunyaev & Zeldovich 1970a;Daly 1991;Hu et al. 1994), though more recent developments have expanded our understanding by providing complementary approaches and accurate numerical treatments (Chluba et al. 2012a;Pajer & Zaldarriaga 2013;Chluba & Sunyaev 2012;Khatri et al. 2012a,b).For the remainder of this subsection, we consider distortions generated by adiabatic fluctuations, following closely Chluba et al. (2012b); Chluba et al. (2015a), while a treatment relevant to isocurvature modes can be found in Chluba & Grin (2013).
The effective heating rate ( ac ) for this process is given by with Θ N being the Nth moment of the photon temperature or polarization transfer function (the superscript P refers to polarization), v the baryon velocity (transfer function),   = 2 2  −3 P  the (dimensionful) primordial power spectrum,  =  T  e  is the rate of Thomson scattering,  the Hubble rate and  the scale factor normalized to be unity at the present time.For primordial distortions, the tight coupling regime ensures that the quadrupole anisotropy dominates the signal, and so we neglect multipoles with ℓ ≥ 2 (Chluba et al. 2012a) as well as additional polarization corrections (Chluba et al. 2015b), which only affect the results at the percent level.The background CMB energy density (at redshift ) scales as   ∝ (1 + ) 4 .For adiabatic initial conditions and in the tight coupling regime (Hu & Sugiyama 1995), it is possible to derive a more functionally useful form of this rate, namely ≈ 0.9 is a dimensionless coefficient related to the neutrino loading and receives various -dependent corrections in the case of isocurvature perturbations (Chluba & Grin 2013).The damping scale  D is determined through (Kosowsky & Turner 1995;Dodelson 2003) and encodes the fact that perturbations with  >  D () have already dissipated their energy into the plasma at a given redshift .
Furthermore, one can find that a given  mode dissipates the majority of its energy at redshift  diss ≈ 4.5×10 5 (/10 3 Mpc −1 ) 2/3 .From here it is straightforward to show that -distortions will be efficiently generated by power in modes of wavenumber 50 Mpc −1 ≲  ≲ 10 4 Mpc −1 , allowing us to constrain regions of the primordial power spectrum complementary to CMB anisotropy measurements.
For a specified form of the visibility function, one can package the distortions in a simple way that highlights the dependence on the form of the PPS by inserting Eq. ( 2) into Eqs.( 3) and ( 4) and computing -space window functions, Here,  min ≈ 1 Mpc −1 is a cutoff scale introduced due to the fact that modes with  ≲  min are tightly constrained by CMB anisotropy, and because the efficiency of energy injection drops for larger scale modes, which dissipate in the post-recombination era through heat conduction/velocity terms rather than shear viscosity (Chluba et al. 2012b).The exact form of the window function depends on our choice of distortion visibility function, of which we have five options with varying degrees of accuracy.We briefly mention each method here but refer the reader to Chluba (2016) for a detailed discussion.
Method A: The simplest approximation one can make is to assume that the  −  and full thermalization transitions take place instantaneously, While the least accurate, this approximation can simplify the computation significantly and hence has been used in simple estimates for anisotropic distortions signals (e.g., Pajer & Zaldarriaga 2012;Ganc & Komatsu 2012).Method B: As a marginal improvement, one can include the fact that even for  >  th , small -distortions can be generated by introducing the thermalization efficiency function, J bb () ≈  − (/ th ) 5/2 .This leads to Given the wide use of this description in the literature, approximations to the -space window functions have also been made (Chluba & Grin 2013;Chluba et al. 2015a) for this Method These can be useful for gaining analytic intuition for e.g.-function features.Note that k = /Mpc −1 .Method C: We can further treat the  −  transition redshift in a more consistent way by simply fitting the  and -distortions to their values computed numerically using a Green's function approach Method D: The previous approach suffers from energy leakage into the residual distortions, and as a result only allows a determination of  and  to within roughly 10% − 20%.To enforce energy conservation, one can choose PCA Method: Each of the previous methods relied on peeling back layers of approximations for the treatment of the microphysics of thermalization.None of these methods, however, was able to model the complicated dynamics of the (non-/non-) residual distortion.In order to capture and utilize this information, one can resort to a principal component analysis (PCA) in which the residual part of a distortion is projected onto a basis of its eigenmodes and subsequently used to constrain energy release histories Here, the distortion signature in a narrow frequency bin  can be decomposed into its contributions from the temperature shift, -, and -type spectra, plus a residual contribution.The   are a set of distortion parameters akin to  and  which can be used to improve the accuracy of the visibility functions.Unfortunately, these improved visibility functions lack a compact analytic expression, though the interested reader can consult Chluba & Jeong (2014), Chluba (2016) and Lucca et al. (2020) for further details.Additional parametrizations based on boosts of the -distortion are discussed in Chluba et al. (2022) but will not be considered here.
For the purposes of distortions from small scale power, the window functions  / () are most relevant.Fig. 1 shows a comparison of the various methods for both the  and -distortion.With the exception of Method A, the variations between procedures are most evident in the low  tail of   and the high  tail of   , due to different treatments of the residual era.In addition to being the most accurate method, the PCA approach provides the most generous window function as it takes a more comprehensive distortion history into account.Below we will mainly focus on comparing Method B and the PCA, highlighting how the constraints on the small-scale power are underestimated using Method B.

Features in the small-scale spectrum
At scales relevent to CMB anisotropy measurements, the spectrum of primordial scalar perturbations is well measured to be adiabatic, nearly scale-invariant, and follow Gaussian statistics (Planck Collaboration et al. 2020).The Planck collaboration has characterized the power spectrum at scales 10 −3 Mpc −1 ≲  ≲ 10 −1 Mpc −1 to be with pivot  p = 0.05 Mpc −1 , spectral index  s = 0.9641, running  run = −0.0045,and amplitude   = 2.1 × 10 −9 .These correspond to the maximum likelihood values found in Tables 4 and 5 of Planck Collaboration et al. ( 2020) (TT,TE,EE+lowE+lensing) when adding a running of the scalar spectral index.Using Eq. ( 6), it is straightforward to show that if this spectrum extends to scales of  ≃ 10 4 Mpc −1 , a spectral distortion with amplitude  ≈ 2 × 10 −8 will be generated.2This distortion presents a target that next generation space missions such as PIXIE 3 (Kogut et al. 2011(Kogut et al. , 2016) ) or Voyage2050 (Chluba et al. 2021) are capable of observing with high significance.Perhaps more exciting, however, is the possibility that these space missions could see a -distortion with a different amplitude, as this would be a smoking gun signal of a departure from the standard ΛCDM cosmology.If one finds limits such that  ≲ 2 × 10 −8 , it would signify a strong red-tilt to the small scale power spectrum, while  ≳ 2 × 10 −8 would imply a boosted feature in the spectrum, or some other exotic energy injection mechanism such as decaying or annihilating dark matter (Bolliet et al. 2021;Liu et al. 2023), cosmic strings (Cyr et al. 2023a,b), primordial black holes (Nakama et al. 2018), or other possibilities (e.g., Chluba 2013b; Chluba & Jeong 2014).
It is therefore possible to use spectral distortions as a tool to constrain departures from the nearly scale invariant power spectrum on small scales.Since distortion signatures are an integrated effect, constraints derived depend sensitively on the shape of the feature in the PPS.The NANOGrav collaboration has hinted that their recent detection could have been generated by some of these features at second order in perturbation theory (scalar induced gravitational waves), and that these models offer a better fit when compared against the expected astrophysical sources (Afzal et al. 2023).While the precise dynamics of GW emission from scalar perturbations poses a highly non-linear and complicated problem (Pen & Turok 2016), in the perturbative regime one can convert the distortion PPS constraints into GW limits.In following Afzal et al. (2023), we thus consider spectral distortion constraints on three different shapes: -function, Gaussian peak, and box-like features in the PPS.Additional cases are discussed in Chluba et al. (2012b).
-function: Here we assume that the feature in the power spectrum is a single -function located at  * , with amplitude Depending on the position of the peak, sizeable distortions can be generated with amplitudes given by  =      ( * ) and  =      ( * ).In Fig. 2, we show constraints on the amplitude of the -function from both COBE/FIRAS (blue) as well as a PIXIE-type satellite (green).We additionally show contours of intermediate distortion constraints.The precise primordial value of  that PIXIE will be able to constrain depends on how reliably one can subtract off the contribution coming from low redshift effects of reionization and clusters (which source a -distortion with amplitude ≃ 10 −6 , Refregier et al. 2000;Hill et al. 2015;Thiele et al. 2022).Through upcoming direct measurements of the SZ effect of galaxy clusters and X-ray observations we expect to be able to model this contribution to the level of 1% − 10%, which we use here as benchmarks.
The left and right plots differ by the window function that was applied in computing the distortion.The left plot makes use of the full PCA procedure and provides stronger constraints in the range of 3 Mpc −1 ≲  ≲ 30 Mpc −1 compared to the right plot.The right plot uses the Method B window function [and its analytic approximations given in Eqs. ( 11) and ( 12)] which are most commonly employed in the literature.Finally, we show the nearly scale-invariant power spectrum inferred by Planck with the red (solid and dashed) line.While it looks like PIXIE cannot see the extrapolated small-scale PPS signal, we remind the reader that integrating this spectrum over the window function produces  ≃ 2 × 10 −8 , which is in reach of these next-generation experiments, assuming that the foreground challenges presented in Abitbol et al. (2017) can be overcome.
In Fig. 3 we show an expanded look at these constraints in the context of a broader experimental reach.On large scales precise measurements of CMB anisotropy provide stringent limits on departures ) of a -function feature in the scalar power spectrum using the full PCA window function (left), and the most commonly used analytic expressions (Method B) for  / ( ) (right).In addition to the nominal COBE/FIRAS and PIXIE-type distortion limits/forecasts, we also showcase a variety of intermediate limits on a next-generation primordial  distortion measurement.Low redshift effects from reionization and clusters produce a -distortion with amplitude ≃ 10 −6 , meaning that primordial  constraints will be dependent on how reliably one can marginalize out that standard model component.Constraints on broader features can be estimated by integrating the corresponding power which consequently tightens the limits, as we show in Figs. 4 and 5.
from a nearly scale invariant spectrum (Planck Collaboration et al. 2020).Strong features in the power spectrum can lead to the abundant production of primordial black holes which are constrained by a variety of astrophysical and cosmological phenomena (Carr & Kuhnel 2022) as indicated by the grey contour.Constraints from scalar induced gravitational waves (SIGWs) can be derived by assuming a non-detection of gravitational wave backgrounds from upcoming telescopes such as SKA, LISA, and the Einstein Telescope (Gow et al. 2021;Moore et al. 2015;Bartolo et al. 2016;Maggiore et al. 2020).Big bang nucleosynthesis (BBN) also drives constraints on the right edge of the distortion contours (Jeong et al. 2014;Nakama et al. 2014;Inomata et al. 2016).We differentiate between current and future constraints by filled and unfilled contours respectively.
As we will discuss in the next section, the recent detection of a stochastic gravitational wave background may be sourced by a strong scalar feature.The detection contours depend sensitively on the shape of the feature, so here we choose to only show the 1 and 2 contours for a -function feature as discussed in Antoniadis et al. (2023b).In general,  features are unphysical (Cole et al. 2022), implying that some width in -space is necessary.Wide features in the power spectrum also typically boost the predicted distortion signature, which means that a PIXIE-type experiment could potentially see a  distortion if the detected stochastic background has a SIGW origin.For example, the broad models illustrated in Figs.18 and 19 of Antoniadis et al. (2023b) lead to an integrated signal of  ≃ 10 −8 − 10 −6 due to their enhanced power at  ≲ 10 5 Mpc −1 , thus becoming visible to future CMB spectrometers, and providing a litmus test for these scenarios.
Gaussian Peak: For our second PPS shape, we consider departures from scale-invariance by the addition of a Gaussian peak in ln , modelled by This model is characterized by three parameters, and we chose to illustrate constraints on the amplitude  Gauss  by benchmarking values of the width Δ, and the peak position  * in Fig. 4. In contrast to the -function feature where we showcased constraints from  and , here we simply focus on  constraints for illustrative purposes.The left hand plot of this figure shows that the constraints become more stringent as the Gaussian widens, as even peak positions of  * ≳ 10 5 Mpc −1 will have significant tails in the active regions of   ().As expected, the constraints approach that of a -function in the Δ → 0 limit.
The right hand plot shows a broader range of widths for peak positions spread throughout active areas of   ().For very narrow widths, the constraints become roughly constant as one approaches the  feature, with an amplitude dependent on precisely where  * lies in the window function.For large widths, significant leakage of the signal in the high and low  modes outside of the window function lead to a weakening of the constraints, coupled with an overall reduction of the effective amplitude at the peak controlled by  Gauss  /Δ.

Box Feature:
We additionally consider a box-like feature, also consisting of three parameters and characterized by where Θ are Heaviside step functions (Not to be confused with the photon transfer functions) and  upper/lower correspond to the upper and lower cutoffs of the box respectively.The amplitude of the distortion can be determined by though it is more instructive to transform to parameters which characterize the central box position ( * = ( lower +  upper )/2), and total width ( w =  upper −  lower ) when discussing constraints on the amplitude.Similar to the Gaussian case, we show these constraints for benchmarked values of  w and  * in the left and right plots of Fig. 5 respectively.For the left hand plot, the constraints gradually strengthen for wider boxes, as expected.The block constraints on the left hand side of the plot arise due to the fact that there  lower ≲ 1Mpc −1 and can be strongly excluded by CMB anisotropy constraints.The right hand plot illustrates this same feature with blocked constraints on the high  w end of the figure .A note worth mentioning is that while the  and Gaussian features are properly normalized to their respective amplitude parameters (in log space), the box feature is not.This leads to weaker constraints for thin boxes ( w ≲ 50 Mpc −1 ) around 10 2 Mpc −1 ≲  ≲ 10 4 Mpc −1 when compared with the  limits.For properly normalized features, the  feature always produces the most conservative bounds.We make this somewhat inconsistent choice of normalization to match those models considered by the NANOGrav collaboration (Afzal et al. 2023) In this case, taking the limit  w → 0 will once again reproduce the most conservative  feature constraints as expected.

Dissipation of tensor perturbations
Scalar perturbations are not the only type of fluctuation that can source spectral distortions at early times, both tensors (gravitational waves) and vectors offer contributions that can be significant in certain regimes.While a spectrum of vector modes sourced before the distortion window ( ≳ 2 × 10 6 ) will decay quite rapidly, interesting constraints can be put on tensor modes by considering their dissipation through CMB polarization fluctuations.This process was originally studied in Ota et al. (2014); Chluba et al. (2015b), while an update regarding the complementarity of distortion based constraints and other more traditional gravitational wave observatories can be found in Kite et al. (2021a) and Campeti et al. (2021).
Tensor perturbations source spectral distortions mainly through free-streaming effects, in contrast to scalar perturbations which are damped through free-streaming as well as direct interactions with the electron-photon fluid.Indeed, it is these direct interactions (through Thomson scattering) that introduced the scalar damping scale  D in Eq. ( 2), efficiently converting small-scale ( ≲  D ()) fluctuations into an effective heating term for the plasma.
The lack of an efficient damping mechanism for gravitational waves introduces two qualitative difference when compared against scalars.First, the conversion from gravitational waves to an effective heating term is suppressed relative to the scalars by roughly 5 orders of magnitude.In contrast, the second effect is that tensor perturbations can actively contribute to spectral distortions over a much wider range of scales, 1 Mpc −1 ≲  ≲ 10 6 Mpc −1 , with a power law decay of the window function in the UV (compared to the exponential suppression of scalars for  ≳ 10 4 Mpc −1 ).
For a given form of the primordial tensor power spectrum, the procedure to compute the associated -distortion is very closely related to the formalism presented in Eq. ( 6) for scalars.If the tensor spectrum is sourced at arbitrarily high redshifts (from inflation, for example), the average value of the  distortion can be determined by In analogy with the scalar sector,  T is the dimensionful tensor power spectrum, and  T  is the -dependent window function for tensors.This window function is computed numerically following the details of Chluba et al. (2015b).
There are also many scenarios in cosmology that produce a gravitational wave background significantly later than inflation, such as first order phase transitions (Kosowsky et al. 1992a,b;Kosowsky & Turner 1993), cosmic strings (Vilenkin 1981;Vachaspati & Vilenkin 1985), metastable domain walls (Hiramatsu et al. 2010;Kawasaki & Saikawa 2011), and scalar induced gravitational waves (Ananda et al. 2007;Baumann et al. 2007), to name a few of interest to the NANOGrav signal.In these cases, the tensor modes will be generated over a specific redshift range on sub-Hubble scales.To account for this, Kite et al. (2021a) defined a generalization of Eq. ( 25), namely Here, W T  (, ) is known as the tensor window primitive, and en-codes time dependent information relevant to the damping of these more general gravitational wave spectra.In the limit where a tensor spectrum appears instantaneously,  T (, ) =  T (), the window primitive is related to the more familiar window function via We show an illustration of the window functions 4 in Fig. 6 for a primordial tensor spectrum ( src ≳ 10 10 , much before the  era), as well as for backgrounds created at later times.The amplitude of the window function is typically 5 orders of magnitude below the scalar counterpart in the plateau region, generally leading to a weaker distortion signature except for the case of a strongly blue-tilted tensor spectrum.The amplitude of the window function decreases for tensor backgrounds created during the  era as there is less time overall for them to dissipate their energy.For  src ≲ 5 × 10 4 , no  distortion is possible and so the window function closes.In principle a window function could also be constructed for -distortions, though cluster marginalization makes detection of a small  signal very difficult once the strong suppression relative to scalar dissipation is accounted for and so we do not consider them here.
For a more in-depth analysis of the formalism, including an explicit form for the window primitive and various mappings between  T and the more familiar gravitational wave observable Ω GW ℎ 2 , the reader is referred to Kite et al. (2021a,b).As indicated in Fig. 10, spectral distortions from the direct dissipation of tensor modes are capable of probing a low frequency region of the GW parameter space not accessible by other experimental efforts.While the constraints from COBE/FIRAS are rather weak, next generation experiments will provide deep insights into tensor backgrounds over a wide range of scales.Any sufficiently strong primordial ( src ≳ 5 × 10 4 ) SGWBs will undoubtedly source spectral distortions, providing a valuable multi-messenger signal that can be useful for model discrimination.

SCALAR INDUCED GRAVITATIONAL WAVES
At linear order in perturbation theory, the scalar-vector-tensor decomposition of fluctuations provides us with a powerful tool to study various phenomena.At second order, however, it is well known that scalar perturbations can source tensor modes (Ananda et al. 2007;Baumann et al. 2007).As this is a second order effect, scalar induced gravitational waves (SIGWs) suffer from significant suppression if the amplitude of the scalar power spectrum remains near its scaleinvariant value at arbitrarily small scales, Ω GW ∝  2  .Therefore, theories which introduce enhancements to the small scale power will also source a spectrum of gravitational waves which can have interesting cosmological consequences.
Perhaps the most straightforward consequence is the production of a stochastic background that could be detected by various different gravitational wave observatories.In light of the recent detection, the NANOGrav (Afzal et al. 2023) collaboration has indicated that a fit to the data with a SIGW component is more convincing than considering SMBHBs alone.Additionally, the EPTA collaboration (Antoniadis et al. 2023b) has determined 1 and 2 "detection" contours for a SIGW progenitor, which we reproduced as the orange region in Fig. 3.These contours assume a -function feature in the scalar spectrum with an amplitude large enough that significant production of PBHs may also occur.Future constraints (assuming a non-detection) can also be forecasted for SKA, LISA, and the Einstein Telescope (see Gow et al. (2021) and references therein).The detection of this stochastic background has inspired a flurry of work on the topic (Cai et al. 2023;Huang et al. 2023;Wang et al. 2023;Zhu et al. 2023;Yi et al. 2023a,b;Firouzjahi & Talebian 2023;You et al. 2023;Balaji et al. 2023;Yuan et al. 2023;Choudhury et al. 2023;Unal et al. 2023;Jin et al. 2023), showcasing a broad interest in the community to further explore the parameter space of SIGWs.
Here, we would like to highlight a synergy between CMB spectral distortions and the gravitational waves induced by scalar perturbations.We choose to take a slightly different approach, where instead of converting gravitational wave detections and non-detections into constraints on the PPS (as is usually done), we ask how CMB spectral distortions and other constraints on the PPS can be transformed into limits on the gravitational wave parameter space for SIGWs.
The SIGW calculation has been refined since the seminal work of Ananda et al. (2007), including the derivation of useful analytic forms for simple shapes in the PPS (Kohri & Terada 2018;Espinosa et al. 2018).These results have been cross-checked (Inomata & Nakama 2019) and an extensive literature now exists on the subject (see Domènech (2021) for a comprehensive review).Here we summarize the results relevant to mapping the CMB spectral distortion constraints onto the gravitational wave parameter space assuming a -function feature (i.e., the constraints presented in Figs. 2 and 3).
We assume that the gravitational waves are induced during radiation domination.During this epoch, tensor modes are primarily sourced close to horizon re-entry for any scalar feature.Modes relevant to CMB spectral distortions generically cross the horizon in the radiation era (recall that  eq ≃ 10 −2 Mpc −1 ).We also assume these scalar fluctuations obey Gaussian statistics.Under these simplifications, the time-averaged primordial tensor spectrum is related (at second order) to the scalar spectrum by the expression (Inomata & Nakama 2019) where  is conformal time (d = d),  = / and  = |k − q|/ relate the wavenumbers of the scalar and induced tensor modes, and  (, , ) is a highly oscillatory kernel encoding the time- dependence of the source function (see Kohri & Terada (2018); Domènech (2021) for some exact expressions).
During radiation domination and in the subhorizon limit, this expression can be approximated as The Θ indicate Heaviside step functions, and the gravitational wave energy density per logarithmic  interval is given by For a  feature such as the one considered in Eq. ( 20), the GW energy density takes a simple analytic form (Kohri & Terada 2018) where k = / s, * (we now add the subscript  s when discussing features related to the scalar spectrum).
The left panel of Fig. 7 illustrates the general shape of the induced GW spectrum, and consists of two interesting features: A resonant peak, occurring at  r = 2 s  s, * ( s = 1/ √ 3 in the radiation era), and an extended IR tail, which has been proposed as a fit to the NANOGrav data (Afzal et al. 2023;You et al. 2023).We make use of both features when performing our mapping of spectral distortion constraints from the PPS.The GW spectrum cuts off at  ≥ 2 s, * due to momentum conservation.
After the production of these GWs, their energy density is simply redshifted to the present day for comparison with observations, where  c and  c are the conformal time and temperature of horizon crossing for a  feature located at  s =  s, * , and Ω r,0 is the radiation density today.The usual effective spin/entropic degrees of freedom are encoded by  * and  *  .
To translate the PPS constraints into limits on the SIGWs we proceed as follows: For a given value of  s (=  s, * ), we read off the corresponding upper limit on    ( s, * ) coming from Planck (large scales) or PBHs, COBE/FIRAS and PIXIE (small scales) as seen in the left panel of Fig. 8.This amplitude is then inserted into Eq.( 30) to infer the maximum allowed Ω GW at a given -mode.The right panel of Fig. 8 illustrates the evolution of these constraints when including more stringent measurements on the small scale PPS from Frequency (Hz) 10 20 10 19 10 18 10 17 10 16 10 15 10 14 10 13 10 12 10 11 10 10 10 9 10 8 10 7 10 6 10 5 10 4 10 3 10 2  2021) and Ellis et al. (2023).Additionally, we include limits from FIRAS and future experiments on direct tensor dissipation as computed in Kite et al. (2021a).CMB spectral distortion limits on the PPS push the constraints at  ≲ 10 −11 Hz down, while PBH limits set the bound at higher frequencies.In the event that LiteBIRD detects a primordial SGWB, a PIXIE-type experiment would be capable of determining if they were produced by SIGWs from a  feature.The shaded orange region represents the 2 detection contours reported by NANOGrav (Agazie et al. 2023c).
COBE/FIRAS and a future spectral distortion mission with PIXIElike sensitivity.
For a PIXIE-type experiment, the dominant constraint over modes of 10 −2 Mpc −1 ≲  ≲ 10 4 Mpc −1 comes from the IR tail of a  feature at the PBH bound on the right hand side of the (scalar) PIXIE contour.Therefore, it is clear that the most efficient way to strengthen constraints on the GW parameter space is by improving limits on the PPS over a wider range of scales.
The current bounds from FIRAS can also be seen in Fig. 9.For small and intermediate -modes, the smooth IR tails of two separate scalar  features (located at  s, * ≃ 2 Mpc −1 and  s, * ≃ 4 × 10 4 Mpc −1 ) set the leading constraints.The regions marked "resonant" are instead set by the resonance peaks of a forest of  features and follow a 1 → 1 mapping onto the PPS constraints.
One subtlety regarding our mapping is related to the regularization of the resonant peak exhibited in Fig. 7.The height of this resonance is in principle unbounded, though the total integrated energy density remains finite.We choose to regulate this amplitude by demanding that the fractional energy neglected in our scheme be  ≲ 0.3%.The height of the resonance peaks, and thus our constraint contours, are dependent on the choice of , so a couple of comments are in order.First, in a realistic setup, the minimum width of this resonance would be bound by finite-time effects in analogy to the finite widths of atomic spectral lines, implying that the amplitude itself is not unbounded.Secondly, gravitational wave detectors are sensitive to energy deposition in -modes with a finite bin width.Our regularization condition on  implies that we flatten the peak over a width of / r ≃ 10 −4 .Provided that this width is smaller than the binning done by any given experiment, our constraint curves will not be altered.We leave a more complete analysis of this regularization to future work.
We show an extended look at the mapping of PPS constraints on the GW parameter space by the blue shaded region in Fig. 10.As expected, the limits we set here do not preclude a SIGW origin to the NANOGrav signal, whose preferred region of parameter space is indicated by the orange contour.Instead, we highlight that stringent constraints can be set on low frequency gravitational waves originating from  features in the primordial power spectrum due to limits on CMB spectral distortions.
LiteBIRD (Allys et al. 2023), CMB-S4 (Abazajian et al. 2019;Abazajian et al. 2022), and the Simons Observatory (Ade et al. 2019) are hunting for large-scale gravitational waves by searching for mode polarization in the CMB.If such a signal is detected, disentangling its possible origin will be a top priority.For instance, if the Lite-BIRD satellite detects a signature of primordial gravitational waves, a PIXIE-like spectral distortion experiment will possess the capability to determine whether the source of such gravitational waves is of a SIGW ( feature) origin.This can be deduced from Fig. 10, noting that the green dashed contour fully covers the LiteBIRD sensitivity curve.In addition, we can immediately conclude that COBE/FIRAS already rules out a significant SIGW contribution over most of the frequencies probed by LiteBIRD.Moving forward, multi-messenger signatures such as this will be crucial in determining the physical origin of various primordial processes.
The constraint contours derived in this section are robust for the case of a  feature in the PPS.They can, however, change quite dramatically for features with more general shapes.We have chosen to show the  constraints first as they offer an intuitive and illustrative mapping between constraints on the scalar spectrum, and those on the tensors.The obtained limits can, however, be viewed as conservative.
The mapping procedure becomes more complicated with general shapes, because both the spectral distortion limits on the PPS, and the gravitational waveform responses transform non-trivially.In general, wider shapes will induce a larger  distortion, which in turn will lead to even more stringent constraints on the amplitude of the induced gravitational waves.We plan to explore this further and expand our results in a future work to more generic log-normal, box, and broken power-law shapes.

PRIMORDIAL BLACK HOLES
The presence of enhanced density perturbations also gives rise to the possibility of PBH formation.Indeed, there have already been numerous claims of a PBH interpretation (or by-product) to the NANOGrav signal (Franciolini et al. 2023;Depta et al. 2023;Inomata et al. 2023;Wang et al. 2023;Bhaumik et al. 2023;Mansoori et al. 2023;Huang et al. 2023), as well as some work scrutinizing the idea (Gouttenoire et al. 2023).
A population of primordial black holes can induce CMB spectral distortions in up to three unique ways.On the low mass end, PBHs in the range 10 −20  ⊙ ≲  PBH ≲ 10 −17  ⊙ will undergo their final moments of evaporation during the distortion era, directly injecting significant amounts of energy and entropy (Acharya & Khatri 2020; Chluba et al. 2020).Spinning PBHs undergo slightly different dynamics, but can be constrained through similar effects (Pani & Loeb 2013).Over a wide range of intermediate mass scales the accretion of background material onto the PBHs can generate distortions through the creation of primordial 'jet'-like features (Ricotti et al. 2008).While these distortions are too faint for COBE/FIRAS to see (Horowitz 2016;Ali-Haïmoud & Kamionkowski 2017), there is still hope for next generation instruments, with the challenge of disentangling the -type distortion from the SZ cluster contributions.Finally, the production of primordial black holes requires large features in , the scalar perturbation.As we have discussed above, if these features exist in a certain range of -modes, they can generate sizeable spectral distortions.We review this calculation (closely following Nakama et al. 2018) and present updated constraints using the more precise PCA scalar window functions as shown in Fig. 1.
For simplicity, let us once again posit that the power spectrum possesses a -function feature at some scale  * , namely P   =    [ln() − ln( * )].Assuming Gaussian statistics (for a non-Gaussian extension see Nakama et al. 2018) the probability density function for the curvature perturbation in a given patch of the sky is given by (Nakama et al. 2018;Carr & Kühnel 2019) where  is the dispersion of the perturbations smoothed over the horizon.For -function features, the amplitude and dispersion are related by    =  2 .Numerical simulations indicate that gravitational collapse of a Hubble patch occurs when the local curvature perturbation exceeds some critical value,  ≳  c .The precise value of  c is still a topic of some debate, but likely lies within the O (0.1 − 1) range.For our purposes, we take  c = 0.67 as found in Harada et al. (2017) to provide a more direct comparison with the results of Nakama et al. (2018) who used this same value.
The fraction of patches (i.e. the initial abundance) which collapse into black holes is determined by where erfc() is the complement of the error function.When a critical 10 2 10 4 10 6 10 8 10 10 10 12 10 14 10 16 10 18 M PBH (M ) density fluctuation re-enters the horizon, that region collapses to form a primordial black hole with some fraction (, which we take to be unity here) of the horizon mass,  pbh =  H ().More precisely, the location of the -function peak determines the (monochromatic) mass of the PBH distribution via (Nakama et al. 2017)  pbh ≃ 10 9  ⊙   * 10.75 scaled here to directly make the link to super-massive black holes apparent.The initial and late time abundance (  pbh = Ω pbh /Ω dm ) are related through Recalling that for  features the induced distortion is  ≃      ( * ), severe constraints can be mapped onto the PBH parameter space from nondetection by COBE/FIRAS.Previous work (Nakama et al. 2018) utilized less precise versions of the window functions, so we present updated constraints using the PCA method described above in Fig. 11.Most notable, the inclusion of information from the  and residual eras can extend previously derived constraints by almost two orders of magnitude.
The results in Fig. 11 show a strong overlap between -and excluded regions.This is primarily due to the fact that a large feature in the PPS is needed to create any primordial black holes (notice how contours extend down to  pbh → 0).This means for PBH constraints that even the tails of the  and  window functions have large constraining power.These tails overlap strongly, unlike the maximal plateaus.The importance of the window function tails also explains why possible gains with PIXIE, while significant in overall distortion sensitivity, are surprisingly modest for  pbh .
Since both  and  can independently constrain wide PBH mass ranges, a simultaneous detection of both distortion parameters (e.g. from PIXIE) would combine to narrow down a precise mass range for PBHs.A larger amplitude in the  distortion compared to the  would signify higher mass PBHs and vice versa.
In Fig. 11 we exclude scales above  rec ≈ 3 × 10 −2 Mpc −1 , since these modes enter the horizon after recombination, and the production of SDs is more complex than the window function description given above.It is also noteworthy that beyond scales of  ∼ 1 Mpc −1 there would also be direct limits on the PPS arising from Planck.
Finally one important extension not considered here is the formation of PBHs from non-gaussianities, thus moving beyond the simple PPS picture.This is considered in more detail in Nakama et al. (2018) and Ünal et al. (2021), where typically speaking the constraints are weaker for non-gaussian fluctuations due to a lowering of the collapse probability at a given PPS amplitude.In the results this impacts not just the mass ranges probed, but also limits the  pbh visibility in some cases (not all contours stretch to arbitrarily low  pbh ), which would give PIXIE a significant boost in constraining power when compared against COBE/FIRAS (unlike in Fig. 11).In addition, one could expect significant distortion anisotropies to form from the stochastic nature of the dissipation process.These scenarios may be probed with  − / correlations (see Rotti et al. (2022) for limits from Planck and Bianchini & Fabbian (2022) for discussions on COBE/FIRAS), which can now in principle be computed using a full spectro-spatial descriptions of the CMB anisotropies (Kite et al. 2022).

DISCUSSION AND CONCLUSIONS
In this work, we have further elucidated the link between CMB spectral distortions, and exotic models which can source a significant stochastic gravitational wave background.After a brief discussion on the generation of spectral distortions from scalar perturbations, we reviewed a simple window function approach to determine the relevant distortion parameters in various approximate schemes (Methods A/B/C/D/PCA).The different window functions can be seen in Fig. 1, and spectral distortion constraints (assuming a -function feature) for methods B and PCA are shown in Fig. 2. Method B is most widely used in the literature, but misses out on constraining power from some dissipation occurring during the residual era (5 × 10 4 ≲  ≲ 3 × 10 5 ).As a result we recommend using the PCA window function in future works to generate more robust and stringent constraints.
Spectral distortion limits on scalar dissipation are an integral constraint, which means that they transform (and in general, become stronger) when the shape of a primordial feature transforms away from a  function.To appreciate how the constraints evolve, Figs. 4 and 5 show limits for Gaussian and Box-shaped features, which have been of recent interest 5 .Perhaps most insightful is the widening of constraints for broad Gaussians, which could dramatically strengthen the GW exclusion plots for SIGWs.
Following this, the formalism from Kite et al. (2021a) was presented to compute the  distortion from direct dissipation of primordial tensor modes.In analogy to the scalar case, a window function approach was derived and generalized to include GW backgrounds sourced at  ≲ 10 8 when their observational signatures could be altered.A summary of the tensor window functions can be found in Fig. 6, and the corresponding limits seen in Fig. 10 highlight the complementary nature of SD constraints in bridging scales between experiments.
From here we moved on to computing the distortion signatures of a non-exhaustive list of exotic models also capable of producing 5 Additional cases can be found in (Chluba et al. 2012b) a sizeable SGWB.We first looked at the case of scalar-induced gravitational waves from -like features in the primordial power spectrum.Contrary to what is conventionally done in the literature, we performed a mapping which took the constraints on the PPS and translated them into conservative limits on the GW parameter space for these  features.A discussion on the subtleties and substructure of these limits is presented in and around Fig. 7, while the results of this mapping and the specific impact of spectral distortion constraints are highlighted in Figs. 8 and 9.
A more complete look at the GW landscape for SIGWs (from  features) was presented in Fig. 10.For  ≳ 10 −10 Hz, the excluded region is driven by a forest of resonant peaks and has a 1 → 1 mapping to PBH constraints on the PPS.Interestingly, in the event that LiteBIRD makes a detection of primordial tensor modes, a next generation spectral distortion experiment would be capable of determining if SIGWs were the source.The overall shape of these limits will change significantly when considering scalar features other than a -function.In general, we expect that for wider features the GW constraints will get stronger, as these shapes will also necessarily widen and deepen the CMB distortion constraints on the PPS.We plan to investigate this in future work.
We additionally rederived the constraints on primordial black holes coming from large enhancements to the power spectrum (assuming Gaussian fluctuations).As an improvement to the previous literature, we utilized the full PCA window function for both the  and  distortions, slightly strengthening the constraints (see Fig. 11).The production of a distribution of PBHs can lead to a significant GW background, particularly in models with large amounts of clustering.CMB spectral distortion signatures of such models are another possible avenue of future research, including the study of their SD anisotropies (e.g., Kite et al. 2022).
The opening of the gravitational wave window is a major step forward in our understanding of cosmology, but comes with a number of new challenges that must be addressed.In particular, the ability to confidently disentangle astrophysical and primordial signals will be crucial for exploration of the parameter space of exotic models.CMB spectral distortions provide a powerful multi-messenger probe, complementary to other more direct measurements, that is capable of discriminating such signals.

Figure 1 .
Figure1.The -space window functions for  (left) and  (right) distortions.We compare the four different analytic methods discussed in the text with the numerically computed PCA projection.Method  treats  th as a hard cutoff which explains the high  discrepancy for   .The other variations between methods come from different treatments of the residual era (5 × 10 4 ≲  ≲ 3 × 10 5 ), with the PCA being the most faithful representation of the full numerical solution to the thermalization procedure, as computed in CosmoTherm.

Figure 2 .
Figure2.Constraints on the amplitude (   =    ) of a -function feature in the scalar power spectrum using the full PCA window function (left), and the most commonly used analytic expressions (Method B) for  / ( ) (right).In addition to the nominal COBE/FIRAS and PIXIE-type distortion limits/forecasts, we also showcase a variety of intermediate limits on a next-generation primordial  distortion measurement.Low redshift effects from reionization and clusters produce a -distortion with amplitude ≃ 10 −6 , meaning that primordial  constraints will be dependent on how reliably one can marginalize out that standard model component.Constraints on broader features can be estimated by integrating the corresponding power which consequently tightens the limits, as we show in Figs.4 and 5.

Figure 3 .
Figure3.An extended view of the primordial scalar power spectrum.The largest angular scales are tightly constrained by CMB anisotropy, with spectral distortions providing competitive bounds at scales of 1 Mpc −1 ≲  ≲ 10 4 Mpc −1(Chluba et al. 2012b).At still larger scales, gravitational wave experiments can set constraints by non-detection of a stochastic background induced through SIGWs (the PTA, SKA, and LISA contours).In addition, large scalar fluctuations can seed abundant PBH production, which are constrained by various astrophysical and cosmological signatures (seeCarr et al. (2010);Carr & Kuhnel (2022) and references therein for specific details).Disruptions to BBN also lead to a set of marginal constraints(Jeong et al. 2014).Current and future constraints on the spectrum are indicated by filled and unfilled contours respectively.The 1 and 2 detection contours reported by EPTA for a ( -function) SIGW signal are labelled in orange(Antoniadis et al. 2023b).

Figure 4 .
Figure 4. COBE/FIRAS constraints on the amplitude of a Gaussian-type feature in the primordial power spectrum.Left: The peak position is varied for different benchmark widths.In the limit Δ → 0 the constraints approach that of the  feature.Right: The peak positions are benchmarked and width varied.

Figure 5 .
Figure 5. Amplitude constraints from COBE/FIRAS for the box type feature in terms of the central position of the box  * , and the total (linear) width,  w .Left: Peak positions are varied for different width benchmarks.The rectangular constraints on the left side of the plot are generated due to stringent constraints coming from CMB anisotropy.Right: Varying the widths for fixed peak positions.The rectangular contours on the right come from the same source as on the left.Constraints appear weaker when compared to the -function feature for thin boxes due to an inconsistent choice of normalization, followingAfzal et al. (2023).

Figure 6 .
Figure6.The window function for tensor perturbations.For a long-lived gravitational wave background (blue line,  ≳ 10 10 ), the window function is maximized as tensors will dissipate over the entire  regime.For backgrounds sourced at later times, the window function degrades as the tensors have less overall time to dissipate their energy.

4Figure 7 .
Figure 7.The spectrum of induced gravitational waves from a -function feature in the PPS.The GW resonance occurs at kr = 2 s ≈ 1.15 for tensors sourced in the radiation era.

Figure 8 .Figure 9 .
Figure8.Left: Upper limits on the primordial power spectrum for three different choices of constraints, Planck + PBH (black), including FIRAS (blue), and including PIXIE (green).Right: The constrained region for the gravitational wave parameter space of SIGWs, assuming a -function feature in the primordial power spectrum.The colour coding of the contours match the upper limits traced out on the left hand plot.

Figure 10 .
Figure10.An extended look at the gravitational wave parameter space for SIGWs from a -function feature in the PPS.The sensitivities of various GW observatories have been compiled fromCampeti et al. (2021) andEllis et al. (2023).Additionally, we include limits from FIRAS and future experiments on direct tensor dissipation as computed inKite et al. (2021a).CMB spectral distortion limits on the PPS push the constraints at  ≲ 10 −11 Hz down, while PBH limits set the bound at higher frequencies.In the event that LiteBIRD detects a primordial SGWB, a PIXIE-type experiment would be capable of determining if they were produced by SIGWs from a  feature.The shaded orange region represents the 2 detection contours reported by NANOGrav(Agazie et al. 2023c).

Figure 11 .
Figure 11.Updated spectral distortion constraints on the  pbh - pbh phase space.The region already excluded by COBE/FIRAS is shown in blue, and future PIXIE-type forecasts are shown in green.Dashed and solid lines show the  and  contours respectively.Note that the COBE/FIRAS  contour extends all the way to the black line on the right, which corresponds to scales entering horizon post-recombination.Both PIXIE lines similarly extend all the way to the right.