Footprints of the QCD Crossover on Cosmological Gravitational Waves at Pulsar Timing Arrays

Pulsar Timing Arrays (PTAs) have reported evidence for a stochastic gravitational wave (GW) background at nHz frequencies, possibly originating in the early Universe. We show that the spectral shape of the low-frequency (causality) tail of GW signals sourced at temperatures around $T\gtrsim 1$ GeV is distinctively affected by confinement of strong interactions (QCD), due to the corresponding sharp decrease in the number of relativistic species. Bayesian analyses in the NANOGrav 15 years and the previous International PTA datasets reveal a significant improvement in the fit with respect to cubic power-law spectra, previously employed for the causality tail. This suggests that the inclusion of Standard Model effects on GWs can have a potentially decisive impact on model selection.


I. Introduction
A stochastic background of Gravitational Waves (SGWB) may be the only direct probe into the early stages of cosmological evolution, where it can be produced by physics beyond the Standard Model (SM).The recently-reported evidence for a nHz SGWB in the NANOGrav 15 years [1, 2] (NG15), European PTA Data Release 2 (EPTA-DR2) [3,4], Parkes [5,6] and Chinese [7] PTAs datasets, whose uncorrelated component was already detected with previous data [8][9][10][11], has attracted ample attention in the astrophysics, cosmology and particle physics communities.
One of the most urgent endeavors is determining whether the signal is of cosmological or astrophysical origin, in the latter case sourced by supermassive black hole binaries (SMBHBs), see e.g.[12,13] for an overview.This is difficult for multiple reasons.First, while the NG15 analysis suggests that the astrophysical model may face challenges [14], the current understanding of SMBHBs is not sufficiently precise to draw conclusions [11,15,16]; Second, the spectrum of a cosmological SGWB generically depends on the microphysical nature of the source and often requires case-by-case numerical simulations.
Our starting point is a fortuitous coincidence of scales: the nHz frequencies probed by PTAs coincide with those of GWs that re-enter the Hubble horizon at the epoch of the QCD crossover phase transition, i.e. at temperatures T ∼ 100 MeV (see e.g.[38]).While the crossover is not expected to source GWs, the rapid drop of relativistic degrees of freedom in the thermal bath significantly changes the equation of state (EoS) of the Universe, that is precisely determined by means of lattice QCD [39].
A causality-limited GW source, active before the QCD crossover, produces a model-independent low-frequency GW signal, which we refer to as Causality Tail (CT), that is affected by the SM-induced change in the EoS (as similarly pointed out for other GW signals in [40][41][42]), due to the different evolution of GW modes whose wavelength is larger than the Hubble radius when the source is active, see also [43][44][45]).
This Letter derives the spectral shape of CT signals at nHz frequencies, that can be readily used by PTA collaborations, GW and BSM communities for model comparison (improving upon simple power-law CTs, currently adopted by NG15 [26] and EPTA-DR2 [16], see also [21,22]).Importantly, we show that our novel inclusion of QCD-induced features significantly impacts the interpretation of current PTA data, by performing a Bayesian search for a CT signal in the International PTA data release 2 (IPTA-DR2) [11,46] and NG15 [2].

II. Standard Model features in the causality tail of primordial GW backgrounds
For cosmological SGWBs, a powerful property is ensured in a broad class of primordial sources where GWs are generated locally, independently in each spatial patch, and in a limited amount of time.We denote by f ct the frequency of GWs entering the Hubble radius when emission shuts off.The remarkable property of the CT is that GWs with frequencies f < f ct evolve independently of the source, because the corresponding wavelengths are larger than the source's correlation length.The evolution of each tensor mode h k (t) in this regime is sensitive only to the expansion of the Universe and GW propagation [41,.Cosmological PTs are a typical example [59,67,77,78].In this case, bubble collisions, sound waves, and plasma turbulence act as causality-limited sources, each with its own finite correlation length.
The GW energy fraction is customarily defined as where we introduced the critical energy density today ρ cr = ρ γ,0 /Ω γ,0 in Eq. (1) (where Ω γ,0 h 2 = 2.47 × 10 −5 is the SM radiation abundance today and h ≡ H 0 /(100 km/s/Mpc) is the reduced Hubble constant).Ω gw (f ) exhibits a spectral shape, that is typically peaked at some frequency f ⋆ > f ct .The CT of the spectrum behaves as Ω gw (f ≲ f ct ) ∝ f 3 in a universe filled by a perfect relativistic fluid with EoS w(t) = p/ρ = 1 3 , see e.g.[51].
This tilt of the CT can be modified as the GW energy density ρ gw (f ) in the CT is determined by the Universe's expansion history and its effect on evolution of tensor modes [43].After the GW source shuts off, the emitted super-Hubble modes freeze due to Hubble friction, remain practically constant until they re-enter the Hubble sphere, and then proceed with under-damped oscillations diluted as 1/a.If the equation of state parameter is constant at this epoch, the CT scales as On the other hand, during the QCD crossover the EoS is not constant, since heavy hadrons form and both the effective number of degrees of freedom in energy g * (T ) = ρ/(π 2 T 4 /30) and entropy densities g * ,s (T ) = s/(2π 2 T 3 /45) decrease rapidly.From the thermodynamical relation sT = p + ρ one then infers the temperature-dependent EoS w(T ) = 4  3 (g * ,s (T )/g * (T )) − 1.The values of g * (T ), g * ,s (T ), and the EoS shown in Fig. 1, are precisely determined with lattice QCD techniques [39].The temporary decrease of w(T ) is due to the pressureless contribution of QCD matter to the thermal bath, followed by its rapid depletion.The CT is thus distinctively modified, when the corresponding GWs with wavenumber k re-enter the Hubble horizon The effect of the temperature dependence of w(T ) on the CT can be only approximately captured (see e.g.[41] in the context of inflationary GWs)by accounting in (1) for the fact that SM radiation is slightly reheated during the crossover, such that: where a(a 0 ) is the scale factor (today), while GWs are decoupled from the thermal bath.
We perform a precise computation by solving the equations of motion for h k (t) (see App. B), accounting for the temperature dependence of w and g * , and obtain the spectrum shown in Fig. 1.A clear deviation from the commonly employed f 3 approximation is found, accidentally located in the nHz range where PTA experiments are most sensitive.
Beside the SM effects presented above, a nHz CT signal probes the cosmic expansion history back to ∼ T QCD .In particular, the possible presence of a fraction f fs ≡ ρ fs /ρ tot of free-streaming species is easily accounted for, by multiplying the Ω gw spectrum shown in Fig. 1 by a factor f 16ffs/5 , for f fs ≲ 15% [43].
The causality tail (CT) of a cgw is commonly modeled as a power-law with n T = 3 (α = 1 2 , γ = 2).The main novelty of our work is that the actual CT is significantly different, when f ct lies above the frequency bins employed in PTA analyses, i.e. 10 nHz where S(f ) is tabulated in App.B and normalised to S(f yr ) = 1, while A ct is the amplitude at f ct .Fig. 1 shows Ω gw (f ) ∝ f 2 S 2 (f ) (i.e. with f fs = 0).A crucial difference between astrophysical and cgw backgrounds is that only the latter contribute to the energy budget in the early Universe, and can affect cosmological observables as any other relativistic freestreaming component beyond the SM.cgws contribute to the effective number of neutrino species as N eff ≡ 3.044 + ∆N cgw eff , with ∆N cgw eff = ρ cgw /ρ ν,1 and ρ ν,1 is the energy density of a single neutrino species.Specifically, the total (integrated) GW abundance is Ω cgw h 2 ≃ 1.6 • 10 −6 (∆N cgw eff /0.28) [59].Measurements of the Cosmic Microwave Background (CMB) [86] and Baryon Acoustic Oscillations (BAO) constrain ∆N eff ≤ 0.28 at 95% C.L.For the peaked sources of interest, Ω cgw ≃ Ω cgw (f ⋆ ) and the constraint on cgw backgrounds reads , (95% C.L.), (7) for signals that can be approximated with power-laws up to f ⋆ .While this assumption is often not valid (see e.g.[87,88] for PTs), the ∆N cgw eff bound can be applied model-independently to the CT, thereby giving A ct ≤ 5 × 10 −14 (f yr /f ct ).This often captures the approximate strength of the constraint, since typically f ct ≲ (0.1 − 1)f ⋆ and the spectrum flattens close to the peak.By comparison with Eq. ( 5), it is evident that CMB constraints can affect signal interpretation when f ct ≳ f yr , if α ≥ −1 (see also [89]).Finally, let us remark that GWs with frequency f ≳ f ct provide an unavoidable contribution to f fs ∼ 0.3∆N cgw eff , see also App.B.

IV. Bayesian Analysis
We now present the results of a Bayesian search for the CT of a cgw signal in IPTA-DR2 [11] and in the recent NG15 [2] dataset.We aim to quantify the impact of the QCD-induced deviations from a power-law with n T = 3 on signal interpretation.Following the IPTA and NG15 collaborations, we limit our analysis to the first 13 and 14 frequency bins of each datasets respectively, to avoid pulsar-intrinsic excess noise at high frequencies.
We consider two cgw models: the CT and the powerlaw n T = 3, both with only one free parameter, A ct and A cgw respectively.We impose ∆N eff constraints in (below) Eq. ( 7) to A cgw (A ct ).For the CT signal, we fix f fs to the unavoidable contribution from GWs of frequency f = f ct (see App. B for details).We set f ct = f yr (such that the first bins of both datasets are deep in the CT) and similarly f ⋆ = f yr for the n T = 3 model, and comment on other choices below.We also compare with common power-law processes (non-cgw) with free exponent γ, as well as with γ = 13  3 (i.e.n T = 2 3 , α = − 2 3 ).For the NG15 analysis, we also compare with the SMBHB expectation obtained by the NG15 collaboration [26], via population synthesis studies assuming circular orbits and GW-only energy loss.App.E lists all parameters and prior choices.
We compare models by computing the Bayes factors B ij , quantifying the preference of model i with respect to model j.These are reported in Tab.I.In this analysis, we consider only auto-correlation terms, rather than the full HD function, because of computational time limitations.We expect only minor effects in our NG15 results when including HD correlation [2].
One important result is the comparison between the CT and power-law n T = 3 signals.As highlighted in bold in Tab.I, the CT provides a better fit than f 3 to both datasets, with a Bayes factor indicating "substantial evidence".
The comparison with a possible astrophysical SGWB depends on the dataset, and is subject to uncertainties on the SMBHB signal.For NG15, no substantial evidence in favor of nor against the CT is observed in comparing with the SMBHB expectation of [26], while the n T = 2 3 and the free power law models are both favored over the CT.The preference for these two power law models is stronger in the older IPTA-DR2.
) for the comparison of the CT with alternative models.The CT spectrum (see Fig. 1) is substantially preferred with respect to Ω gw ∼ f 3 .The ∆N eff upper bound has been imposed on the amplitude of both CT and f 3 spectra, fixing f ct = f yr for the CT and f ⋆ = f yr for n T = 3.
The GW spectra and CP time delays for maximumposterior values of parameters obtained by our analyses are shown in Fig. 2 (here we include HD correlations in the NG15 results, using the faster fitting procedure described in [90]).Notice that for A ct = A cgw the CT allows for a larger amplitude of Ω gw and of CP delays in the first bins (roughly by a factor of 1.5).
The curves in Fig. 2 correspond to negligible values of f fs from GWs of frequency f = f ct .Additional contributions from high-frequency and other freestreaming species are model-dependent, but would affect our results only if ∆N cgw eff ≳ 0.1.In this case, an excess in upcoming cosmological surveys [91][92][93] is expected.
For a qualitative comparison with simple power-laws, one can approximate the CT as a power-law with slope evaluated through a given number of frequency bins.The result of this procedure is shown by the vertical gray lines in Fig. 3, obtained by fitting the CT with a power-law through the lowest 4th to 8th frequency bins of NG15, where evidence for GWs is reported (and within the 9 frequency bins employed by EPTA-DR2).The 1, 2 and 3σ posteriors on power-law parameters reported by IPTA-DR2 and NG15 are also shown, together with the results of EPTA-DR2 (using the DR2new HD dataset).We expect the CT to provide as good a fit to the NG15 (EPTA-DR2) data as a power-law model with parameters inside the 3σ (1σ) region of the posteriors, whereas the f 3 signal lies at the border of the 3σ (2σ) region.We stress that this is only approximate, since the CT deviate significantly from a power-law at the relevant frequencies.
Let us comment on how different f ct values affect our model comparison.The CMB bound on power-law cgws becomes stronger for f ct > f yr , as shown in Fig. 3 for peak frequency f ⋆ = 1, 3, 10f yr .Thus, we find that the IPTA-DR2 dataset decisively disfavors f ct ≳ 32 nHz, while for the new datasets the CT signal is excluded at 3σ (affected) only if f ct ≳ 100 nHz (≳ 50 nHz).The future reach of CMB-S4 [93] observations is also plotted.

V. Discussion and Implication for Particle Physics
We have pointed out that SM physics predicts a specific shape for the CT of a GW signal, clearly distinguishable from Ω gw ∼ f 3 already with present datasets.This applies to any GW source that is active before QCD confinement (i.e.T ≳ GeV, see App.C), and provides a much-needed signature to help determining whether the current PTA excess is of cosmological or astrophysical origin.If the SMBHB interpretation (or a different cgw signal) is preferred, then our work should prove useful in the future to disentangle additional contributions from cosmological sources (see e.g.[94,95]).
We have quantified that, despite the CMB bound on ∆N eff [89], a CT spectrum can explain the PTA signal and is substantially favored over the commonly employed n T = 3 approximation in both IPTA-DR2 and NG15 datasets, while providing a similar fit to NG15 as simple SMBHB models.
The detection of such a signal would have dramatic implications for particle physics, most likely pointing to the breaking of a (global or gauge) symmetry in a dark sector (DS), for instance via a PT (or causing the annihilation of domain walls for discrete symmetries) at temperatures T ≳ GeV.Some macroscopic properties of the corresponding GW source can be estimated rather model-independently (see App. (e.g. the sub-Hubble size of the source, see App.C), the posteriors (2σ) in Fig. 3 imply the lower bound , for 15 nHz ≲ f ct ≲ f yr .The required fraction is even larger for larger f ct .This may be problematic for sources with ϵ ⋆ ≪ 1, as the DS would then dominate over the SM background (α ⋆ ≃ 1), and thus suggests that the source should be relativistic and Hubble-sized.We discuss further the PT scenario and its challenges in App.C (see also [23,25,96]).
Such a large abundance should decay to the visible sector (or be in the SM, if e.g.QCD confinement proceeds via a first order PT in the early Universe, although this scenario requires additional ingredients, see e.g.[97]), because if it remains in dark relativistic species it contributes as [23] ∆N DS eff ≃ 0.28 where we normalised to the current 2σ constraints from CMB.Interestingly, the required couplings between the SM and the dark sector may be additionally probed at colliders, laboratory experiments and astrophysical environments.For non-CT cgw signals, NG15 gives a 2σ bound A cgw ≳ 10 −14.6 to explain the current excess, thus requiring α ⋆ ≳ 0.07ϵ , well within the reach of upcoming CMB and LSS surveys [91][92][93]98].
Finally, in our work we have conservatively employed the SM prediction for the EoS of the Universe.We notice that the measurement of the CT at PTAs offers a probe of the whole cosmological expansion history up to T QCD , as encapsulated by Eq. ( 3).Modifications to the CT could occur in BSM scenarios where w(T ), g * (T ) or f fs are varied, or if the Universe undergoes a phase of matter domination below the QCD crossover (we present a search for this scenario in App.D).Additionally, PTAs could test the EoS of hot QCD matter in the early Universe (see e.g.[97,99]), in the presence of a CT signal.The near-future detection of a SGWB at PTAs could disclose a unique window around the epoch of the QCD crossover.As highlighted in this paper, this coincidence of scales provides a robust feature to assist in the discrimination between the astrophysical and cosmological origin of the SGWB.

A. Relation between frequency today and temperature at horizon crossing
We derive here the present frequencies corresponding to modes crossing the horizon at a given time, in order to clarify a potential confusion with the typical peak frequency of primodial SGWBs and to highlight the coincidence of frequencies probed by PTAs and modes corresponding to the QCD crossover.
The frequency (redshifted to its present day value) of modes crossing the horizon at a given epoch can be simply derived as follows.We impose that, at the temperature of horizon crossing one has k = aH, Using the conservation of entropy g * ,s (T )T 3 a 3 (T ) = constant, the ratio of scale factors between GW production at a temperature T , and the present time at T 0 , is T 0 T ≃ 8.0 × 10 −14 g * ,s (T ) 100 where we have used T 0 ≃ 2.4×10 −13 GeV for the photon temperature today [104] and g * ,s (T 0 ) ≃ 3.91 for the Standard Model degrees of freedom with three light neutrino species [105].Additionally, employing the Friedmann equation GeV is the reduced Planck mass, the energy density in a radiation dominated era ρ = (π 2 /30)g * (T )T 4 and g * ,s (T ) ≃ g * (T ), one obtains Eq. ( 3).
If one considers the emission of GWs from a cosmological phase transition (or from e.g.domain wall annihilation), the characteristic frequency f ⋆ of the SGWB, i.e. the frequency at which one expects the SGWB to peak, is located around the inverse typical time or length scale of the system, which is the inverse of the duration of the PT or of the bubble size, depending on the details of the source.Following Ref. [59], we set k ⋆ /a ⋆ ≃ 2π β, where β is the inverse of the nucleation timescale.Therefore, one obtains the characteristic frequency today as Notice that this frequency is larger by a factor of 2π than the frequency in Eq. (3) corresponding to horizon crossing at a given T , because they refer to two different physical quantities.Eq. (A3) gives the typical peak f ⋆ of the SGWB and identifies the wavenumber k = 2π/λ = 2πf associated to the excitation of a GW over a period f −1 ≃ β −1 , whereas Eq. ( 3) determines the wavenumber k = a(T )H(T ) crossing the Hubble radius at a given time, and allows to fix the boundary f ct of the CT.

B. Tabulated causality tail
This Appendix provides more details about the CT that we compute and use for our analysis.For completeness, we provide a summary of the dynamics for low frequency GW modes excited by a local, causality-limited source in the early Universe (see [43] for details).The equations of motion in conformal time τ for a GW mode h k (τ ) sourced by sector with energy density ρ s and dimensionless anisotropic stress Π read (we understand the space-time tensor indices and the polarisation (+, ×) indices for both the GW and source) where H = a ′ /a is the conformal Hubble time and ′ denotes conformal time derivatives.At small k, the source takes a simple and universal form when it lasts for a finite amount of time (up to a time τ ⋆ ) and has a finite correlation length (≲ 2π/k ⋆ , where k ⋆ ≡ H(τ ⋆ ) is the mode crossing the Hubble radius at τ ⋆ ).This corresponds to GW periods longer than the duration of the source, so that the time dependence of the source is basically a Dirac delta function, and to wavenumbers where the power spectrum of the source is independent of k (and is given by white noise, meaning that sources in different Hubble patches are uncorrelated).Therefore, in this limit J(τ, k) → J ⋆ δ(τ − τ ⋆ ), so that an initially vanishing GW mode gets a velocity h ′ k = J ⋆ shortly after τ ⋆ , and thereafter evolves according to We solve Eq. (B2) for a wide range of momenta k < k ⋆ by accounting for the expansion history of ΛCDM with the SM content around the QCD crossover.We notice that a positive detection of these features on the CT allows to confirm that the cosmological expansion history up to T ∼ T QCD has followed the predictions of ΛCDM, and possible deviations would hint at a different evolution of the scale factor.We can then compute the GW power spectrum P h ∼ ⟨(a(τ )h k (τ )) 2 ⟩ after horizon crossing, and Ω ct ∝ f 5 P h (f ) [43].
To include the effect of free-streaming (high-frequency) GWs on the CT, one starts with the modification on Ω gw h 2 = S 1 (f )f 16fFS/5 .Here S 1 is a function that describes the CT for f fs = 0. Then we transform to the characteristic strain h c (f ), as follows which can then be rewritten as where S is the function tabulated in Tab.II.When f ct = f yr , then A ct is the amplitude of the signal at f = yr −1 , as in our paper, whereas generally it is just the amplitude at f ct .With this notation, the ∆N eff constraint on Ω gw h 2 at f ct reads which implies Now let us relate f fs to A ct .The definition of ∆N eff is where ρ ν,1 is the energy density in one neutrino species.At CMB, this is related to the total energy density in radiation by from entropy conservation, therefore we have where f fs = ρ gw /ρ rad and we assume that the background is always dominated by SM radiation between GW production and matter-radiation equality.Using the above relation between ∆N eff and Ω gw h 2 , we then have where we have normalized g ⋆ to its value at temperatures corresponding to f ct = f yr .Plugging this into (B4), we obtain the GW-altered causality-tail characteristic strain as a function of two parameters: f ct and A ct .In our runs, we shall fix f ct and we are thus left with only one free parameter.
Table II contains the tabulated result, parameterised through the function S(f ) defined in Eq. ( 6), which expresses the characteristic strain as a function of frequency in the PTA range.This function can be directly used as an input for Bayesian parameter estimation.We also include the tabulated curve for the corresponding Ω ct (f ) ∝ f 2 S(f ) 2 shown in Fig. 1.

C. Implications for Particle Physics
We collect here various important considerations from cosmology and particle physics related to our work, and in particular some consequences of the microscopical interpretation of the GW primordial signal.We focus especially on the SGWB generated by the prototypical example of a causality-limited signal, a first-order phase transition.
a. Peak frequency for a causality-limited GW signal around the QCD phase transition.In this work we are studying the impact of the QCD era1 on modes that are excited by a local GW source when they are still outside the Hubble horizon, i.e. k QCD ≪ a(T QCD )H(T QCD ), where we denote by the subscript "QCD" quantities evaluated when the temperature of the universe was T = T QCD ≃ 150 MeV.The mode k QCD corresponds to frequencies f QCD ≡ k QCD /(2π) ≪ aH/(2π).Therefore, for consistency, we have to require the emission of frequencies around f QCD ≈ 3 nHz to take place when the corresponding modes were super-horizon, i.e. when the temperature of the universe was larger than T QCD .Requiring f QCD to be emitted at temperature larger than the QCD transition, forces the peak frequency to f ⋆ ≳ 24 nHz (β/H ⋆ ), neglecting the small impact of the change of effective degrees of freedom.Hence, a detection of a CT signal at PTAs would necessarily imply T ⋆ ≳ 300 MeV, or the peak of the signal would lie in the first PTA bins.Future PTA measurements, while rapidly decreasing their sensitivity at frequencies below the inverse of the observation time [107][108][109], will significantly improve their sensitivity at frequencies larger than the presently constraining bins, thanks to the expected improvement in the instrumental noise in detectors like FAST [110] and MeerKat [111].
The possibility of a PT occurring in a DS not far from the SM QCD crossover might be motivated for example in proposals addressing the approximate coincidence of dark matter and baryonic energy density in the context of asymmetric dark matter [112][113][114][115][116][117][118].
b. Constraints on the evolution of the dark sector responsible for the SGWB.The GW signal investigated in our work is largely independent of the nature of the GW source, as long as the latter occurs on sub-horizon scales and is active only up to a definite epoch in the early Universe.While this is obviously a virtue, it also implies that the detection of a CT signal at PTAs would not be able to uniquely identify the microphysical mechanism which acts as a GW source.To this aim, complementary probes would then be essential.
Nonetheless, some general conclusions can be drawn from a detection of a CT, which apply to a broad class of sources.These necessarily involve physics beyond the SM, which can be organized in a non-specified dark sector (DS).Several scenarios of interest can then be described by using effective macroscopic parameters: the temperature T ⋆ at which the GW source is active and the fraction of the (radiation) energy density that it comprises, α ⋆ ≡ [ρ DS /(3H2 M 2 p )] T =T⋆ .The peak frequency of the GW signal at the time of emission is then determined by the inverse characteristic time scale of the source in units of the Hubble rate: this is given in Eq. (A3), where we can define in greater generality a dimensionless parameter δ ⋆ ≡ f ⋆ /H ⋆ (that takes the specific form β/H ⋆ for PTs, while for domain walls in the scaling regime δ ⋆ ≃ 1).Typically, this time scale is smaller than or comparable to one Hubble time, thus δ ⋆ ≳ 1.The peak amplitude of the signal today can be estimated as [59] Ω where g * (T ⋆ ) is normalized to T ⋆ ≃ GeV, and ϵ ⋆ ≲ 1 generically scales as an inverse power of δ ⋆ (depending on the microphysics of the source) and also accounts for additional signal-suppressing effects (such as subluminal velocities).
For simple power-law signals below the peak frequency f ⋆ , combining with Eq. ( 5) one has for α ⋆ ≪ 1.This estimate yields a lower bound on α ⋆ for CT signals, by replacing f ⋆ → f ct , A cgw → A ct .In particular, for the 2σ posteriors shown in Fig. 3 lead to α ⋆ ≳ (0.1 − 0.4)ϵ −1/2 ⋆ , for 15 nHz ≲ f ct ≲ f yr .For larger values of f ct , this bound becomes even stronger.As mentioned in the main text, the required large fraction of energy density to fit the current excess cannot remain secluded from the SM bath below T ⋆ because of CMB and BBN constraints.This leaves only one possibility: the dark sector should deposit a significant amount of its energy density into the SM bath shortly after the emission of GWs.
This conclusion generically remains true when the source is a first order PT in a dark sector, even when the peak frequency lies among the first frequency bins of PTAs.Indeed, for PTs δ ⋆ ≃ (β/H ⋆ ) ≥ 1 is the rate of variation of the bubble nucleation rate and is typically O(100) (or possibly smaller, depending on the exact time-dependence of thermal effective potential), while ϵ ⋆ ∼ δ if the signal is sourced by bubble collisions (sound waves), where v w ≤ 1 is the bubble wall velocity (see e.g.[59]).Therefore, achieving such a large GW signal from a PT at a scale not far from 100 MeV requires a large value of α ⋆ , 2 and a mildly efficient damping of the DS energy density into the SM thermal bath before BBN.Such a coupling between the SM and new particles of O( GeV) mass could lie in a wide range, as it would not suffer from significant constraints from supernova emission for such masses [120], and could be probed in principle by laboratory experiments.
A scenario that can potentially evade the constraints above, while still exhibiting a CT, is that of GWs induced by the collapse of large density fluctuations [70,[121][122][123][124][125].This can occur if the curvature power spectrum is enhanced at scales k ≳ 10 6 Mpc −1 , much smaller than those probed by CMB anisotropies [20,24,29,30,[126][127][128][129][130][131][132][133].However, for the CT to be relevant in the induced GW scenario, the peak of the SGWB should be reached at frequencies higher then those observed with PTA experiments, and with a peak amplitude h 2 Ω gw (f ⋆ ) ≳ 10 −8 (see also [24]).Therefore, the required large scalar perturbation amplitude would generate an unacceptably large relic abundance of PBHs (see e.g.[134,135] for reviews), unless negative non-Gaussianities suppressing PBH formation are present [136] (overproduction of PBHs in Gaussian models has been shown to significantly constrain the interpretation of the excess already in the IPTA-DR2 and NANOgrav 12.5 years datasets in Ref. [24]).Thus, a viable explanation of the PTA signal by second order GW requires the scalar source to be active around the horizon crossing of nHz modes, with a different specific imprint from the QCD crossover era (see e.g.[66,137]).
Finally, let us stress that our discussion does not generically apply to sources that are active on time scales much larger than one Hubble time, such as cosmic strings.These sources are not expected to exhibit a CT, and their GW signal is also not characterized by Eqs.(A3) and (C1).However, the interpretation of the PTA excess in terms of strings generated by the spontaneous breaking of a global symmetry (such as arise in axion models, notice that if the axion is massive then the string network disappears once the corresponding domain walls form and the GW signal then features a CT as usual) faces constraints from cosmology, see [138].
c. Impact of the bound from ∆N eff on the amplitude of Ω gw (f ⋆ ).This consistency condition, that only applies to the scenario where the source for GW emission is a cosmological phase transition, becomes relevant when confronted with the requirement of falling below the bound from ∆N eff .Given the amplitude required to explain the PTA signal, by extrapolating the CT we would infer f ⋆ ≲ 60 nHz.By accounting for the relation between the peak frequency and the duration of the transition β −1 , this would force us to β/H ⋆ ≲ 1.3 • (150 MeV/T ⋆ ).On top of this upper bound, there is a lower bound imposed by the our assumption about the CT lying around the QCD transition, T ⋆ > T QCD .Given that β > H ⋆ , these equations would imply T ⋆ ≃ [150,200] MeV.
This back-of-the-envelope estimate is, however, overly restrictive for various reasons.i) First, the parameter space for our assumptions to be consistent becomes larger if one considers phase transitions with non-relativistic bubble wall velocities v w < 1.Smaller velocities would lead to an additional scaling of the peak frequency that follows f ⋆ ∝ 0.35/(1 + 0.07v w + 0.69v 4 w ) or f ⋆ = 0.536/v w for the bubble or sound waves contribution respectively.ii) Secondly, the approximately cubic scaling dictated by causality arguments is expected to break down before the peak frequency, and the SGWB is much less steep between f ct and f ⋆ .This alleviates the upper bound from ∆N eff .We show this possibility in Fig. 4, where we adopt the following parametrisation for the contribution of bubble 10 -2 10 -1 10 0 10 1 10 2 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 10 1 FIG.4: Spectra of SGWB induced by first order PT.We normalise the frequency in terms of the peak frequency f ⋆ while factorise out the overall spectral amplitude.We show both the expected contribution of bubble collisions and sound waves.In green we indicate the causality tail obtained assuming perfect radiation domination, while the black dashed line indicates the effect of an intervening matter domination (MD) phase that would shift the peak to larger frequencies.
collisions [139] Ω with (a, b, c) = (3, 1, 1.5) [21] and sound waves [140] Ω Here, for simplicity, we are neglecting the effects of the QCD crossover on the CT, as they do not affect this discussion.Due to the breaking of the f 3 scaling close to the peak, Ω gw is lower than the one obtained extrapolating the CT up to peak frequency of a factor O(6).
iii) Lastly, an additional model dependence is still possible.In particular, as also discussed in App.D, one may envision an intermediate matter dominated phase that could effectively break the cubic scaling with an intermediate linear growth Ω gw (f ) ∝ f .As a consequence, depending on the duration of the intermediate matter phase, the bound on the SGWB amplitude in the CT could be relaxed, effectively shifting the onset of the cubic drop-off.

D. Effects of intermediate matter domination
The interpretation of our results presented in the discussion in the main text are based on Eqs.(A3) and (C1), which assume a standard radiation-dominated (RD) Universe until matter-radiation equality.Changes to the expansion and including HD correlations (green).For the bending frequency, the value used in Fig. 5a for IPTADR2 is 25 nHz, corresponding to T i ≃ 800 MeV, while for NG15 we used ≃ 7 nHz, corresponding to T i ≃ 300 MeV.For comparison, we also plot the model-agnostic broken power-law with white noise in A cp at high frequencies, as in the analyses of the IPTA-DR2 collaborations [11] (the curves are obtained by letting the bending frequency, as well as the low-frequency slope and the amplitude free to vary, notice that the difference of our results with respect to those of the collaboration is due to the fact that we are using the maximum posterior values of each parameter, rather than the maximum likelihood signals).
We then compare the models in the IPTA-DR2 dataset with the choices discussed above.The resulting Bayes factor is log 10 B CT+MD, BPL ≃ 1.2, showing that a CT signal affected by matter domination below the QCD crossover fits the full IPTA-DR2 30 bins dataset better than a broken power-law with white noise at high frencies, although the latter might still be considered as a more motivated model.The 2d posteriors resulting from fitting separately the two models to the data are shown in Fig. 5b.We find: log −0.43 .Following the above considerations, it is important to notice that additional suppression factors in both the peak frequency (A3), ∼ r 1/3 , and amplitude in (C1), r 4/3 (if the source is active during MD) arise.Requiring the CT to be in the first bins of the PTA sensitivity band implies that its amplitude is O(10 −7 − 10 −8 ) at f yr .In turn, this requires α ≳ (0.1 − 0.3)r −2/3 .On the other hand, successful BBN imposes T f ≳ 5 MeV, i.e. r ≳ 0.005 for f ⋆ ≃ 10 −7 Hz, with the inequality becoming weaker (stronger) if f ⋆ is larger (smaller).Therefore the GW source could have overcome the radiation background, if the MD phase lasted until temperatures corresponding to the first frequency bins of PTA datasets.

E. Details of the numerical analysis and additional posteriors
Here we report details of our numerical analysis.The analysis to obtain Bayes factors has been performed using using the codes enterprise [100] and enterprise extensions [101], in which we implemented the CT signal reported in Fig. 1, and use PTMCMC [102] to obtain MonteCarlo samples.We derive posterior distributions using GetDist [141].We include noise parameters (white, red and dispersion measure) according to the strategy of the IPTA DR2 [11].For more details, see also the Appendix of [23].For the NG15 dataset, we have used PTArcade [90].
Our prior choices for noise (DM only applies to IPTA DR2) and GW parameters are reported in Table III.In order to produce the results reported in the main text, we have restricted our analyses to the first 13 and 14 frequency bins of the IPTA-DR2 and NG15 datasets respectively, as in [2,11].We obtained more than 10 6 (typically 5 • 10 6 ) samples per each analysis and model in this work, and discarded the first 25% of each chain.With these choices, we reproduce the posteriors for CP parameters presented in [11] with excellent agreement.
The upper prior boundary on the amplitude of cosmological signals is set by the ∆N eff bound reported in Eq. ( 7), setting f ct = f yr for the analyses reported in the main text to set ideas.
For the NG15 dataset, the Bayes factor with respect to SMBHBs model has been computed using the prior on the astrophysical signal provided in PTArcade, which is based on population synthesis assuming circular orbits and

FIG. 1 :
FIG.1: Top: EoS during the QCD crossover.The inset shows w in an expanded range of temperatures.Bottom: Impact of the variation of w(T ), g ⋆ (T ) on the CT of a primordial SGWB (plotted with an arbitrary amplitude).The dashed line shows the f 3 scaling obtained in a pure radiation-dominated universe.
FIG.2:The GW spectrum (upper ) and the induced CP time delays S ab (f )/T (lower ), T being the observation time of each PTA.We compare CT (solid), n T = 3 (dashed) and free power-law (dotted) models, selecting the maximum posterior values obtained with IPTA-DR2[11] (blue) and NG15 [2] (green) datasets.The shaded region shows the CMB+BAO bound on ∆N eff , and the dotted line the reach of CMB-S4 experiments.The lower limits of the posteriors are determined by the priors of[2,11].See App.E for parameter posteriors.

FIG. 3 :
FIG.3: IPTA-DR2[11], NG15 and EPTA-DR2 posteriors for a power-law model, the last two enforcing HD correlations in the analyses[2,4].The vertical lines indicate approximations of the CT.The yellow star shows the best fitting CT amplitude according to NG15 data.We show current and future ∆N eff bounds affecting any power-law signal extending up to f ⋆ = 1, 3, 10 f yr .

TABLE II :
Tabulated values for the Causality Tail computed in this work, accounting for the effect of the SM QCD phase transition.The last column assume negligible f FS .