The second data release from the European Pulsar Timing Array V. Implications for massive black holes, dark matter and the early Universe

The European Pulsar Timing Array (EPTA) and Indian Pulsar Timing Array (InPTA) collaborations have measured a low-frequency common signal in the combination of their second and first data releases respectively, with the correlation properties of a gravitational wave background (GWB). Such signal may have its origin in a number of physical processes including a cosmic population of inspiralling supermassive black hole binaries (SMBHBs); inflation, phase transitions, cosmic strings and tensor mode generation by non-linear evolution of scalar perturbations in the early Universe; oscillations of the Galactic potential in the presence of ultra-light dark matter (ULDM). At the current stage of emerging evidence, it is impossible to discriminate among the di ff erent origins. Therefore, in this paper, we consider each process separately, and investigate the implications of the signal under the hypothesis that it is generated by that specific process. We find that the signal is consistent with a cosmic population of inspiralling SMBHBs, and its relatively high amplitude can be used to place constraints on binary merger timescales and the SMBH-host galaxy scaling relations. If this origin is confirmed, this is the first direct evidence that SMBHBs merge in nature, adding an important observational piece to the puzzle of structure formation and galaxy evolution. As for early Universe processes, the measurement would place tight constraints on the cosmic string tension and on the level of turbulence developed by first-order phase transitions. Other processes would require non-standard scenarios, such as a blue-tilted inflationary spectrum or an excess in the primordial spectrum of scalar perturbations at large wavenumbers. Finally, a ULDM origin of the detected signal is disfavoured, which leads to direct constraints on the abundance of ULDM in our Galaxy.


Introduction
The recent observation of a common signal with excess power in the nanohertz frequency ranges (i.e. a common red signal, as defined in Arzoumanian et al. 2020;Goncharov et al. 2021;Chen et al. 2021) in pulsar timing array (PTA) datasets, with emerging evidence for quadrupolar correlations1 opens a new era in the exploration of the Universe.This important milestone has been achieved thanks to the efforts of the European Pulsar Timing Array (EPTA, Desvignes et al. 2016), the Indian PTA (InPTA Joshi et al. 2022), the North American Nanohertz Observatory for Gravitational Waves (NANOGrav, McLaughlin 2013), the Parkes PTA (PPTA, Manchester et al. 2013), and the Chinese PTA (CPTA, Lee 2016).Although the significance of the signal does not yet reach the 5σ mark, which is generally accepted as the 'golden rule' for a firm detection claim, the evidence reported by the different collaborations ranges between 2σ and 4σ (EPTA and InPTA Collaborations 2023b;The Nanograv Collaboration et al. 2023;Reardon et al. 2023;Xu et al. 2023), strongly suggesting a genuine gravitational wave (GW) origin.Awaiting decisive confirmation within the International PTA (IPTA, Verbiest et al. 2016;Perera et al. 2019) framework, with the additional contribution of the MeerKAT PTA (Miles et al. 2023), we are hearing, for the first time, the faint murmur of the GW Universe at frequencies of 1-to-30 nano-Hz, which is ten orders of magnitude lower than the frequencies currently probed by groundbased interferometers (Abbott et al. 2016).This opens a completely new window onto the Universe, allowing us to look at different phenomena, probe new astrophysical and cosmological sources, and, potentially, new physics as well.
By monitoring an array of millisecond pulsars (MSPs) for decades with a weekly cadence, PTAs are sensitive to GWs in the 10 −9 -10 −7 Hz range (Foster & Backer 1990).At those frequencies, the most anticipated signal to be detected is a stochastic GW background (GWB) produced by the incoherent superposition of waves emitted by adiabatically inspiralling supermassive black hole binaries (SMBHBs, Rajagopal & Romani 1995;Jaffe & Backer 2003).The signal manifests as a stochastic Gaussian process characterised by a power-law Fourier spectrum of delays-advances to pulse arrival times, with characteristic interpulsar correlations of general relativity identified by Hellings & Downs (1983).The statistical properties of the signal are expected to significantly deviate from the typical isotropy, Gaussianity and perhaps even stationarity that is typical of many stochastic signals from the early Universe (e.g.Sesana et al. 2008;Ravi et al. 2012).In fact, due the shape of the SMBHB mass function and their redshift distribution, the overall signal is often dominated by a few sources, and particularly massive, nearby SMBHBs might result in loud enough signals to be individually resolved as continuous GWs (CGWs, Sesana et al. 2009;Babak & Sesana 2012;Kelley et al. 2018) emerging from the GWB.The exact amplitude and spectral shape of the spectrum are intimately related to the cosmological galaxy merger rate and to the dynamical properties of the emitting binaries forming in the aftermath of the merger process (Enoki & Nagashima 2007; Kocsis & Sesana 2011;Sesana 2013a,b;Ravi et al. 2014).Therefore, the demonstration of an SMBHB origin of the signal observed by PTAs provides invaluable insights into the formation, evolution, and dynamics of these objects.Moreover, it brings decisive evidence that SMBHBs merge in nature, thus overcoming the 'final parsec problem' (Milosavljević & Merritt 2003), which is still an open question in our understanding of galaxy and structure formation.
Beyond SMBHBs, a number of processes (potentially) occurring in the early Universe can also produce a stochastic GWB at nanohertz frequencies.Tensor modes can be produced as early as during the first tiny fraction of a second after the Big Bang through quantum fluctuations of the gravitational field stretched by the accelerated expansion of inflation (Grishchuk 1975;Rubakov et al. 1982;Starobinskii 1985;Fabbri & Pollock 1983;Abbott & Wise 1984).In the literature, these GWs are referred to as 'primordial'.In this case, the shape of the power spectrum is defined by the specific model of inflation.Classical tensor mode production invoking the presence of a source term in the GW equation of motion can also take place in the early Universe.There are a plethora of physical processes that can lead to such a source term, and trigger the production of GWs.Topological defects, for example decaying cosmic string loops (Vilenkin & Shellard 2000;Damour & Vilenkin 2001;Jones et al. 2003;Dvali & Vilenkin 2004;Damour & Vilenkin 2005), particle production during inflation (Sorbo 2011; Barnaby et al. 2012;Cook & Sorbo 2013;Anber & Sorbo 2012), (magneto-)hydrodynamic turbulence ((M)HD, Kamionkowski et al. 1994;Kosowsky et al. 2002;Dolgov et al. 2002;Caprini & Durrer 2006;Gogoberidze et al. 2007;Caprini et al. 2009), the collision of bubble walls during a first-order primordial phase transition (Kosowsky et al. 1992;Kosowsky & Turner 1993;Caprini et al. 2008;Huber & Konstandin 2008;Jinno & Takimoto 2017;Cutting et al. 2018), sound waves in the aftermath of a first-order phase transition (Hindmarsh et al. 2014(Hindmarsh et al. , 2015(Hindmarsh et al. , 2017)), as well as scalar perturbations at second order in cosmological perturbation theory (Baumann et al. 2007;Ananda et al. 2007), are among the most commonly considered scenarios.GWs decouple from the rest of the matter and radiation immediately after their generation at essentially any temperature in the Universe, so that in the case of clear observational evidence for these types of signals, we can infer nearly unaltered information on the physical processes occurring during or just after the birth of the Universe (Caprini & Figueroa 2018).
Contrary to the incoherent superposition of GWs from SMB-HBs, the stochastic GWB from sources in the early Universe is usually assumed to be statistically homogeneous and isotropic, unpolarised, and Gaussian (Allen 1996;Maggiore 2000;Caprini & Figueroa 2018).Statistical homogeneity and isotropy are inherited from the same property of the FLRW Universe.The absence of polarisation holds provided no macroscopic source of parity violation is present in the Universe.Gaussianity follows by the central limit theorem in most cases of GWBs generated by processes operating independently in many uncorrelated, subhorizon regions.This also applies to the irreducible GWB generated during inflation in the simplest scenarios, when the tensor metric perturbation can be quantised as a free field, and hence with Gaussian probability distribution for the amplitude.There are, however, notable exceptions, among which for example rare GWB bursts from cosmic strings cusps and kinks (Damour & Vilenkin 2000, 2001), or the GWB from particle production during inflation (Cook & Sorbo 2013;Sorbo 2011;Anber & Sorbo 2012).Therefore, although statistical properties might be useful for discriminating between SMBHBs and several processes in the early Universe, a full assessment of the nature of the GW signal will not be trivial.
Spatially correlated delays of the time of arrivals (TOAs) in an array of MSPs are not a unique imprint of GWs.For example, it is well known (e.g.Tiburzi et al. 2016) that such delays can emerge due to the imperfect fitting of the solar system ephemerides (dipolar correlated noise), or due to a miscalibration of the time standard to which the measured TOAs are referred (monopolar correlated noise).Furthermore, individual Fourier harmonics of a common signal in PTA data may include contributions from the oscillations of the gravitational potential associated with the presence of ultralight dark matter (ULDM, Smarra et al. 2023) 2 , also known as fuzzy dark matter (FDM), in the Galactic halo (Khmelnitsky & Rubakov 2014).The existence of ultralight scalars, generally motivated by string-theoretical frameworks (Green et al. 1988;Svrcek & Witten 2006;Arvanitaki et al. 2010), is also particularly appealing from the astrophysical and cosmological point of view.In fact, several potential issues in the small-scale structure of the Universe, such as the cusp vs core (Flores & Primack 1994;Moore 1994; Karukes, E. V. et al. 2015) or the missing satellite (Klypin et al. 1999;Moore et al. 1999) problems, could be disposed of or, at least, mitigated assuming that dark matter is made of ultralight particles.As predicated by Khmelnitsky & Rubakov (2014), the presence of ULDM induces harmonic delays in the arrival times, with a frequency proportional to the ultralight boson mass.
In this paper, we provide a broad overview of the implications of the signal observed in the second data release of the EPTA+InPTA (EPTA and InPTA Collaboration 2023) for the different physical processes mentioned above.More in-depth analysis of several of these scenarios will be the subject of separate future publications.Unless otherwise stated, we consider each process separately, and we discuss the implications of the detected signal under the hypothesis that it is generated by that specific process.We do not attempt any Bayesian model selection on the signal origin, although a general framework for that is being developed (Moore & Vecchio 2021).The main reason for this choice is that, at this stage of data taking and analysis, the information carried by the signal is not particularly constraining; the evidence of the measurement is still at ≈ 3σ, and the amplitude and spectral shape of the signal are not very well measured.With these premises, the result of any model selection is bound to be severely influenced by the priors employed for each of the models under examination.This exercise becomes more meaningful as data get more informative, which we strive to achieve with the analysis of the third release of the combined IPTA data, which is now being assembled.
The paper is organised as follows.In Sec. 2, we describe the signal observed by EPTA+InPTA and its main features, including its free spectrum and best-fit parameters.We then proceed with detailing the main implications of the detected signals under the assumption that it is generated by a cosmic population of SMBHBs (Sec.3) or by a number of processes occurring in the early Universe (Sec.4).In Sec. 5, we investigate the compatibility of the observed signal with a DM origin and place constraints on ULDM candidates.Finally, in Sec.6, we summarise our main results and discuss future prospects.

The observed signal in the EPTA DR2 dataset
Our investigation is based on the results reported in EPTA and InPTA Collaborations (2023b, hereinafter PaperIII), which analyses the data of 25 MSPs collected by the EPTA using five of the largest radio telescopes in Europe: the Lovell telescope at the Jodrell Bank Observatory, the Nançay decimetric radio telescope, the Westerbork synthesis radio telescope, the Effelsberg 100 m radio telescope, and the Sardinia radio telescope.The dataset also includes the Large European Array for Pulsars (LEAP) data, in which individual telescope observations are coherently phased to form an equivalent dish with a diameter of up to 194 m (Bassa et al. 2016).These data are complemented by low-frequency observations of a subset of 10 MSPs performed by the InPTA using the upgraded Giant Metrewave Radio Telescope (uGMRT) and covering about 3.5 years.
The data of each individual pulsar are combined as described in EPTA and InPTA Collaboration (2023) and the noise properties of each pulsar are then extracted according to the optimised custom noise models presented in EPTA and InPTA Collaborations (2023a).The final result is a dataset of unprecedented sensitivity spanning up to 24.7 years.Four versions of the dataset were analyzed: 1. DR2full.24.7 years of data taken by the EPTA; 2. DR2new.10.3 years of data collected by the EPTA using newgeneration wide-band backends; 3. DR2full+.The same as DR2full, but with the addition of InPTA data; 4. DR2new+.The same as DR2new, but with the addition of InPTA data.
The analysis presented in this paper refers to the DR2new dataset only.We do not consider DR2full and DR2full+ because evidence of quadrupolar correlation (usually referred to as HD correlation, from Hellings & Downs 1983) of the common process is weaker in those datasets, potentially due to the lower quality of early data that were collected with narrowband backends (see discussion in PaperIII).On the other hand, although the analysis of DR2new+ produced results in broad agreement with DR2new, that dataset was assembled relatively recently and has not been analysed as thoroughly.For example, the binned freespectra that we will use in some of the following analyses have only been produced after this work was completed.
Before proceeding with the description of the signal detected in DR2new, here we summarise some notations used in PTA analysis for the benefit of the reader.The perturbation affecting the TOAs, whether produced by GWs or DM, is described in terms of its dimensionless strain h.A broad-band stochastic perturbation is defined by its characteristic dimensionless strain h c ( f ), often modelled as a power law For example, a population of GW-driven circular SMBHBs produces a spectrum with α = −2/3 and amplitude A ≈ 10 −15 , assuming f 0 = 1yr −1 .h c ( f ) is connected to the differential energy content of the signal per logarithmic frequency through the equation: where H 0 is today's Hubble expansion parameter.We note that h c ( f ) and Ω( f ) provide equivalent parametrizations of the spectrum.The former is more popular in the astrophysics domain, whereas the latter is the preferred choice for early Universe and cosmology.
Given h c ( f ), the one-sided power spectral density induced by the GW signal in the timing residuals is given by (Lentati et al. 2015): where γ = 3 − 2α.PTAs search for HD correlated time delays with such a power spectrum in the data, and measure the parameters A and γ.For an observation timespan T , measurements are discretised in frequency bins ∆ f i = f i+1 − f i , where f i = i/T .It is then customary to convert S ( f ) in RMS residual induced in the TOAs in each frequency bin: The main properties of the GWB signal observed in DR2new and examined in this paper are shown in Fig. 1.The length of the dataset is T =10.3 years, and excess common correlated power is detected in several frequency bins up to ≈ 30 nHz (Fig. 1 left panel).Conversely, some bins are unconstrained, which results in a relatively loose determination of the spectral properties of the observed signal.In the literature, h c and S in Eqs.(1) and (3) are usually anchored to the pivotal frequency f 0 = 1yr −1 .The data are, however, most informative at the lowest frequencies, while the common power at 1yr −1 is essentially unconstrained.This naturally leads to a strong degeneracy of the A − γ 2D posterior, as shown for example in Figure 1 of PaperIII.Therefore, unless otherwise stated, we change the reference frequency to f 0 = 10yr −1 , where the data are actually constraining, which results in a weaker dependence of A upon γ, as shown in the right panel of Fig. 1.
In the following three sections, we discuss three possible contributions to the signal, probing completely different epochs and scales of our Universe, and the implications for the associated physical processes.Namely, the cosmic population of SMB-HBs (at redshifts z ≲ 1), the early Universe (z > 1000), and DM (within our Galaxy).

Implications I: supermassive black hole binaries
A cosmic population of SMBHBs is the primary astrophysical candidate to produce a signal in the nanohertz band detectable by PTAs.If we define d 5 N/(dzdm 1 dqdedt r ) as the cosmic merger rate of SMBHBs as a function of redshift, primary black hole mass, mass ratio, and eccentricity, the general form of the generated GWB as a function of observed frequency f can be written as (Sesana 2013a) Here, dt r /dln f K,r is the time spent by the shrinking binary within a given logarithmic orbital frequency bin, which converts the merger rate into the distribution of rest-frame orbital frequencies of the emitting population.This quantity depends on the physical processes driving the evolution of the SMBHBs including, besides GW emission, interaction with the stellar and gaseous environment surrounding them.As such, it is generally a function of the binary parameters m 1 , q, e, and extra parameters describing the environment, such as the stellar density in the nucleus of the galaxy host (for more details, see Sesana 2013a).The second line of Eq. ( 5) is the strain amplitude produced by each individual eccentric SMBHB binary, cast as the sum of harmonics fulfilling the condition f K,r = f (1 + z)/n.In that expression, h( f K,r ) is the strain of the equivalent circular binary emitted at twice the orbital frequency of the system, as given in Eq. 4-7 of Rosado et al. (2015), and g(n, e) is a combination of Bessel functions (see, e.g.Bonetti & Sesana 2020, for details).For a distribution of circular, GW-driven binaries, the only relevant mass parameter is the chirp mass M = (m 1 m 2 ) 3/5 (m 1 + m 2 ) 1/5 , and Eq. ( 5) takes the familiar form (Sesana et al. 2008) This can be recast in terms of the comoving number density of merging binaries d 2 n/(dzdM) (Phinney 2001) which highlights that, in this case, the expected spectrum follows a power law h c ∝ f −2/3 , and the only free parameter is its overall amplitude.The latter is set by the function d 2 n/(dzdM), which contains all the astrophysical knowledge of the cosmic population of merging SMBHBs, and is determined by the cosmological hierarchical assembly of galaxies and their central SMBHs.
Conversely, in the most general case described by Eq. ( 5) there is also information in the spectral shape of the signal, since coupling with the environment as well as eccentricity affect the function dt r /dln f K,r , suppressing signal at the lowest frequencies.Moreover, eccentricity distributes the GW power at several higher harmonics of the orbital frequency, slightly modifying the power-law behaviour at high frequencies.In general, the GWB cannot be cast in term of a simple analytical form, although a broken power-law is a sufficient approximation for most situations (see, e.g.Ravi et al. 2014;Sesana 2015;Kelley et al. 2017;Chen et al. 2017b).
The literature investigating the GWB produced by a population of SMBHBs is vast, dating back to the mid-nineties and early 2000s (Rajagopal & Romani 1995;Jaffe & Backer 2003;Wyithe & Loeb 2003;Sesana et al. 2004), and predictions have been made by employing different models and techniques.Models can be broadly classified into two categories: self-consistent theoretical models for SMBH evolution within their galaxies (Sesana et al. 2008(Sesana et al. , 2009;;Ravi et al. 2012;Kulier et al. 2015;Kelley et al. 2017;Bonetti et al. 2018a;Siwek et al. 2020;Izquierdo-Villalba et al. 2022), and empirical models based on observed properties of galaxy pairs coupled to SMBHhost galaxy relations (Sesana 2013b;Rosado et al. 2015;Ravi et al. 2015;Simon 2023), or on the evolution of the SMBH mass function inferred from observations (McWilliams et al. 2014).Note that we group both semianalytic models (SAMs) and large cosmological simulations in the first class.The main difference between these two classes is that self-consistent models are constructed to reproduce a large array of observations related to galaxies and the SMBH they host, such as the redshift-dependent galaxy mass function, quasar luminosity function, and so on.Conversely, empirical models are, by construction, consistent with the observations upon which they are based, but are usually not tested against independent constraints.As a consequence, they can generally produce a wider distribution of GWB amplitudes, but consistency with other observations is not necessarily guaranteed.
In this Section, we investigate the implications of the signal observed in the DR2new dataset for representative models of the two classes.In Sec.3.1, we perform a semi-quantitative comparison between the measured signal and predictions of an extended version of the Rosado et al. (2015) models (hereinafter RSG15) including binary eccentricity and environmental coupling.In Sec.3.2, we exploit the framework developed in Middleton et al. (2016); Chen et al. (2017aChen et al. ( , 2019) ) to draw inference on SMBHB astrophysics from the data, either by assuming astrophysical priors from independent observations, or by using a completely generic model for the SMBHB mass function with minimal assumptions.In Sec.3.3, we demonstrate how the measured signal can inform galaxy and SMBH formation models by examining its constraining power on two state-of-the-art SAMs, namely L-galaxies (Henriques et al. 2015) and the model developed by Barausse and collaborators (Barausse 2012;Bonetti et al. 2018b;Barausse et al. 2020).We discuss caveats and future directions of improvement in Sec.3.4.

Qualitative analysis of empirical SMBHB population models
To carry out a semi-quantitative comparison between observations and empirical models, we use an extended set of SMBHB populations based on the work of Sesana (2013b) (S13 hereinafter) and RSG15.
Fig. 2: GWB amplitude distributions predicted by the RSG15 models.The thin-dashed yellow line is for the full set of models in RSG15, whereas the thick-dashed orange line is for the subset considered here.The solid blue line is the distribution predicted by the 108 downselected sample used in the analysis.The vertical line marks the median value of A at f 0 = 1yr −1 reported in PaperIII when fixing γ = 13/3 in the search.Note that the lower x-axis scale is for A at f 0 = 1yr −1 , whereas the upper x-axis is for A at f 0 = 10yr −1 (the normalization used in this paper).Since α = −2/3 for circular GW-driven binaries, there is a shift of 0.666 dex between the two.

Description of the models
As described in S13, the models are constructed around the number density of merging SMBHs per unit primary mass, mass ratio, and redshift, d3 n/(dm 1 dqdz) obtained by combining different observations of the galaxy mass function and pair fraction, estimated galaxy pair merger timescales, SMBH-host galaxy relations, and prescription for SMBH accretion during mergers (see Section 2 of S13 for full details).Given the large uncertainties in all of the ingredients, the models predict a broad distribution of expected GW amplitudes, as shown in Fig. 2.
Guided by the relatively large amplitude of the detected signal and by theoretical and observational advancements in the last decade, we select a sub-sample of those models, as we now justify.First, hydrodynamical simulations of merging galaxies at different scales as well as deep X-ray observations of merging galaxies support accretion activation onto the individual SMBHs prior to merger (e.g.Capelo et al. 2017;Koss et al. 2018).Moreover, hydrodynamical simulations of sub-pc scale binaries, have found most of the accretion to occur on the secondary (i.e. less massive) SMBH (Farris et al. 2014).We therefore restrict the analysis to models where SMBHs accrete prior to the merger, either with an equal Eddington ratio or with preferential accretion onto the secondary. 3Second, observations of overmassive black holes in the centre of large ellipticals (McConnell et al. 2011) has led to an upward revision of the SMBH-galaxy relations.Contrary to S13 and RSG15, here we consider only those revised realations, namely the ones reported by Kormendy & Ho (2013); McConnell & Ma (2013); Graham & Scott (2013).Finally, given the large number of models, to save computing power, we perform an ad hoc down-selection of 108 models that preserves the overall distribution of the expected GWB amplitudes, as shown in Fig. 2.
As opposed to RSG15, we go beyond the circular-GW driven binary approximation and consider the self-consistent evolution of SMBHBs within their stellar environment.4This is done by employing the hardening models of Sesana (2010) that selfconsistently evolve the SMBHB semimajor axis and eccentricity under the combined effect of stellar hardening and GW emission, once a given initial eccentricity e 0 at binary formation is given.Those evolutionary tracks allow us to evaluate the term dt/dln f K,r in Eq. ( 5), and thus to reconstruct from the density distribution of merging binaries, d 3 n/(dm 1 dqdz), the numerical distribution of systems emitting at any time in the whole sky as a function of mass, mass ratio, redshift, orbital frequency and eccentricity, d5 N/(dm 1 dqdzd f de).For each model, we consider 10 values of e 0 = 0, 0.1, ..., 0.9 and three different normalizations of the stellar density profile, described as ρ = C × ρ 0 (r/r 0 ) −1.5 , with C = 0.1, 1, 10 (details in Sesana 2010).
In total, we explore 108×10×3= 3240 models, spanning different eccentricities and densities of the stellar environment.For each model, we draw 100 Monte Carlo samplings of the distribution of the emitting binaries, we discretise the frequency domain in bins of ∆ f = 10.3yr−1 , and add the signals of binaries falling in the same bin in quadrature.This leads to the binned characteristic strain spectrum h c ( f ) that we then convert in S ( f ) and RMS residuals using Eqs.( 3) and (4).The full procedure is described in Amaro-Seoane et al. ( 2010); Bonetti & Sesana (2020).In this way, we generate a grand total of 324k Monte Carlo realizations of the predicted GW spectrum in the PTA band.

Comparison with the observed signal
The binned spectrum shown in Fig. 3 contrasts expectations from the 324k models (green) to the measured correlated signal in DR2new (orange).The two sets of violin plots are in good agreement in the few lowest frequency bins, where measurements are the most constraining.Note that the model prediction distributions are highly non-Gaussian and asymmetric, with long tails extending upwards.This is due to the fact that sparse very massive/nearby binaries can sometimes produce exceptionally loud signals, as illustrated by the 100 individual GWBs overplotted to the violins.In fact, this might explain the extra power measured in the 4th and, most strikingly, in the 9th lowest bins compared to the bulk of the model predictions.We caution that the 9th bin is close to the 1yr −1 mark, where PTAs are blind due to fitting for the Earth orbital motion, and leakage from imperfect fitting might affect that measurement.In any case, if this extra power is indeed due to GWs, it can be easily accommodated by theoretical models, as demonstrated by the realization highlighted by the tick grey line.
Our Monte Carlo approach to generate the SMBHB population and its associated GW signal also allows us to investigate the occurrence of CGWs in the data, for which evidence in DR2new is found to be inconclusive (EPTA and InPTA Collaborations 2023c).Since the search performed in that paper was limited to circular binaries, we only carry out this analysis for the 32.4k models with e 0 = 0. 5 A full assessment of the detectability of CGWs requires the evaluation of the detection probability of each individual binary for a given false alarm rate, as detailed in RSG15.For the sake of simplicity, and given the qualitative nature of this analysis, we just compute the signal-to-noise ratio (S/N) of each individual binary according to Eq. (46) of RSG15 (thus also restricting to the Earth term only).When computing the S/N of a source, we model each pulsar noise by using the maximum likelihood values of the single pulsar noise analysis presented in EPTA and InPTA Collaborations (2023a), and add the GWB produced by all of the other binaries to the noise spectral density.We arbitrarily set the detectability threshold at S/N= 3 in the following.
Results are shown in Fig. 4, which compares the power distribution of resolvable CGWs to the binned spectra of the overall predicted GW signal and of the DR2new measurements.In line with RSG15, the probability of detecting a CGW is maximum at the lowest frequency, rapidly decaying to less than 0.01 past the 6th bin.Although this seems to suggests that the feature observed at the 9th frequency bin is unlikely to be due to a CGW, it should be noticed that we considered here a threshold of S /N = 3.The overall GW signal in our data has a total S /N ≈ 3. 5 − 4 (EPTA and InPTA Collaborations 2023b), mostly accumulated at the lowest frequency bins.It might still be the case that the feature at the 9th bin is due to an unresolved CGW with S/N< 3, which would be a more common occurrence in the data.Note that the average S/N of CGWs slightly increases at higher frequencies, which is primarily due to the frequency dependence of the CGW characteristic strain, h c ∝ f 7/6 .If p i is the probability of having a CGW of S/N > 3 in the i-th bin, we can compute the probability of detecting at least one CGW with S/N > 3 in DR2new according to these models as p = 1 − (1 − p i ), which gives p = 0.49.It is therefore quite possible that the observed signal is dominated by few bright sources, which might be resolvable in the future with more sensitive datasets.Thus far, searches for CGWs on the current dataset yielded inconclusive evidence (EPTA and InPTA Collaborations 2023c).This probability is obviously S/N threshold dependent.For example, by increasing this threshold to S/N> 5, we get p = 0.13.This is comparable to the 6% chance found by Bécsy et al. (2022).and the slightly larger probability in our models is likely due to the louder overall amplitudes of the signals considered here.We stress, however, that these findings apply to models where binaries remain essentially circular.The number of resolvable CGWs, tends, in fact, to slightly decrease when the eccentricity increases (Truant et al. in preparation).
Finally, we once again propose the comparison first shown by Middleton et al. (2021), who contrasted the measured 2D A − γ posterior to model expectations.For the latter, we just fit the 9 lowest frequency bins of the GWB spectrum of each Monte Carlo realization of the Universe with a straight line in the logA − log f plane.As described in Sec. 2, we normalise the measurement to f 0 = 10yr −1 , where our data are informative, which alleviates the A − γ degeneracy in the posterior.Results are shown in Fig. 5.Although the measured spectrum tends to be shallower than the theoretical one (see also Sec. 3.4), the contours overlap at 2σ and the marginalised amplitudes are broadly consistent.

Inference on the SMBHB population.
After checking the general compatibility of the observed signal with expectations from empirical SMBHB population models, we take a more quantitative approach and exploit Bayesian inference to constrain the properties of SMBHBs from the data.We repeat the analysis carried out by Middleton et al. (2021) on the common signal detected in the NANOGrav 12.5 year dataset (Arzoumanian et al. 2020), exploiting the framework laid out in Middleton et al. (2016); Chen et al. (2017aChen et al. ( , 2019)).The SMBHB population of a given model M is described by a set of parameters {θ}, we then use Bayesian inference to find the posterior distribution p(θ|d, M) for the parameters θ of model M given the observation data d: Here, p(θ|M) is the prior distribution on the model parameters, p(d|θ, M) is the likelihood of model M with parameters θ of producing the data, and p(d|M) is the evidence.For each set of θ, the likelihood is computed by comparing the signal amplitude produced by M at frequencies f i = i/(10.3yr),i = 1, ..., 9, to the posterior distribution of the correlated signal measured in DR2new at the same frequencies.We treat each bin as independent, therefore the likelihood takes the form where p(A, f i ) is the probability density of the amplitude A of the correlated signal measured in the i-th bin, and it is evaluated at the value A M predicted by the model.We estimate the likelihood in each bin using a KDE estimate of the DR2new posteriors, similar to the method used in Moore & Vecchio (2021).
Note that the models we use are deterministic in the sense that they have a 1:1 correspondence between θ and the predicted h c ( f ).In reality, given θ, the ensemble of SMBHBs generating the signal is statistically drawn from the deterministic distribution function, which results in a significant scatter of h c ( f ), as demonstrated by the individual spectra shown in Fig. 3.We caution that this variance is not captured by these models, and its inclusion in the Bayesian inference pipeline is the subject of ongoing work.
As in Middleton et al. (2021), we use two models to describe the SMBHB population, which are described in turn below.

Agnostic SMBHB population model
The agnostic model, developed in Middleton et al. (2016), makes minimal assumptions about the underlying population of SMB-HBs.Binaries are assumed to be circular, GW-driven and the characteristic strain is computed according to Eq. ( 7) where the source distribution is described by a Schechter function (Schechter 1976) in z and M, where t R is the time in the source frame and we assume cosmological parameters from Planck18 (Planck Collaboration et al. 2020a).The five model parameters are θ = {ṅ 0 , α M , M ⋆ , β z , z 0 }, where ṅ0 is the merger rate per unit rest-frame time, co-moving volume, and logarithmic M interval, and the parameter pairs {α M , M ⋆ } and {β z , z 0 } control the shape of the M and z distributions, respectively.The integration limits in M and z are 10 6 ≤ M/M ⊙ ≤ 10 11 and 0 ≤ z ≤ 5, respectively.The prior ranges of the five model parameters are identical to those used in (Middleton et al. 2021).

Astrophysically-informed SMBHB population model
The astrophysically-informed model was developed in Chen et al. (2019).This model captures the interaction between the SMBHBs and their environment and allows for eccentric orbits, both of which lead to a characteristic amplitude that does not follow a simple single power law, as in Eq. ( 5).The model has 18 parameters, 16 of which describe astrophysical observables linking the number of SMBHB mergers to the number of galaxy mergers.The galaxy stellar mass function is modelled as a redshift-dependent Schechter function defined by five parameters: {Φ 0 , Φ I , M 0 , α 0 , α I }.Both the galaxy pair fraction and merger timescales have power law dependencies on the primary galaxy stellar mass M, mass ratio q and redshift function (1 + z) and each of them is defined by a set of four parameters: { f 0 , α f , β f , γ f } for the pair fraction, and {τ 0 , α τ , β τ , γ τ } for the merger timescale: Galaxy pairs are then populated with SMBHs of mass m following a standard black hole to stellar bulge mass relation of the form where N{x, y} is a log normal distribution with mean value x and standard deviation y, which adds three further parameters, {M * , α * , ϵ}, to the model.The final two parameters describe the eccentricity at SMBHB pairing {e 0 } and the density of the stellar environment {ζ 0 }.For the 18 parameters listed above, in the analysis presented here, we use the extended prior intervals listed in Table I of Chen et al. (2019).

Results of the inference
The main results of the inference are shown in Figs. 6, 7 and 8. Fig. 6 shows the marginalised posterior distribution for the normalization of the merger rate density ṅ0 from the  the posterior is very similar to the prior (see Appendix A for full posterior distributions for both models).Fig. 7 displays the posterior on the SMBHB chirp mass function for the two models integrated over the redshift range 0 < z < 1.5.Although the agnostic model results in a loosely constrained mass function, the measured PTA signal alone places interesting lower limits on the SMBHB binary merger rate in the Universe.For example, we can say at 95% confidence that for each comoving Gpc 3 , there have been at least 10 4 SMBHB mergers with M ≈ 10 7 M ⊙ since cosmic noon.When astrophysical priors weights in, the mass function of merging SMBHBs is well constrained by the PTA signal and, as expected from Fig. 6, it lies in the upper range of the astrophysically informed prior.
Within this model, the DR2new measured signal implies there have been about 10 6 SMBHB mergers for each comoving Gpc 3 , with M ≈ 10 9 M ⊙ since z = 1.5.This points towards a very active merger history of massive galaxies and a very efficient dynamical evolution of the SMBHBs forming in the merger process.
For the astrophysically informed model, DR2new also provides interesting information on several model parameters.This is because the astrophysical prior considerably narrows down the range of signal amplitudes allowed by the model, and the measured signal pushes towards the upper bound of this range.This results in informative constraints on several key parameters, related in particular to the SMBHB merger timescale and the SMBH-bulge mass relation.As shown in Fig. 8, the SMBH merger timescale τ 0 following galaxy pairing must be shorter than ≈ 1 Gyr (90% confidence), with the data mildly favouring shorter merger times for massive galaxies at low redshifts (i.e.α τ , β τ < 0).Moreover, the data favour a high normalization of the SMBH-bulge mass relation log 10 M * ≈ 8.4 +0.18  −0.32 , compared to a much wider prior range extending all the way down to log 10 M * = 7.8.This is in line with the qualitative analysis of Sec.3.1, which showed that the signal is consistent with recent, upward-revised, SMBH-galaxy relations.There is also a slight preference for a high normalization of the pair fraction f 0 with a positive z dependence, β f > 1.Despite the low γ value favoured by the data, indicative of a flatter spectrum compared to the canonical value predicted by circular GW-driven binaries, SMBHB dynamics is largely unconstrained, perhaps with a marginal preference for eccentric binaries evolving in dense environments (e 0 and ζ 0 posteriors slightly rising towards the right bound of the prior).
Altogether, these results point towards efficient orbital decay of SMBHBs in the aftermath of galaxy mergers, providing direct evidence that the 'final parsec problem' is solved in nature and that compact sub-parsec SMBHBs must be common in the centre of massive galaxies.

Implications for SAMs
We now explore the implication of this signal for the joint modelling of the galaxy and SMBH formation and evolution by taking a close look at two state-of-the-art SAMs: the model constructed by Barausse and collaborators (Barausse 2012;Sesana et al. 2014;Antonini et al. 2015;Klein et al. 2016;Bonetti et al. 2018b; Barausse et al. 2020) and L-Galaxies (Henriques et al. 2015;Izquierdo-Villalba et al. 2022).In this preliminary study, we do not model the dynamical evolution of the binaries and we assume them to be circular and GW driven, thus resulting in a characteristic strain spectrum with α = −2/3.

SAMs and SMBHB delays
In Fig. 9, we show this comparison for the model of Barausse (2012) in its original version (B12) and its subsequent evolutions, which were used to produce the results of Klein et al. (2016) Barausse et al. (2020) (B+20).Besides the specific SAM implementation and (astro)physics, these models mainly differ for the physical prescriptions describing the delays between galaxy and MBH mergers, with (i) models "LS-nod (B12)", "HS-nod (B12)", "Q3nod (K+16)", "LS-nod-noSN (B+20)", "HS-nod-noSN (B+20)" and "HS-nod-SN-high-accr (B+20)" assuming no such delays (except for the delays between the mergers of the halos and those of the corresponding baryonic components)6 ; (ii) models "popIII-d (K+16)", "Q3-d (K+16)", "LS-d (B+18)", "HSd (B+18)" additionally introducing the effect of stellar hardening, triple MBH interactions and gas-driven migration; and (iii) models "LS-noSN-d (B+20)", "LS-SN-d (B+20)", "HS-noSNd (B+20)" and "HS-SN-d (B+20)" accounting for even longer delays (including large contributions from SMBHB separations of hundreds of pc).Note that the labels "SN" (and "noSN") refer respectively to a putative effect of SN feedback on black Fig. 9: Predictions for the GWB characteristic strain amplitude at f = 1/10yr from a range of SAMs published in the literature, assuming quasicircular orbits and no environmental interactions (i.e.γ = 13/3), but different physical prescriptions for the delays (increasing from left to right) between galaxy mergers and black hole mergers.The "no delays", "medium delays" and "long delays" models correspond respectively to the classes of models (i), (ii) and (iii) defined in the text.The ranges account for the finite resolution of the models.The shaded area is the DR2new 95% confidence bound.More details about the models are provided in the text.
hole accretion (and the absence thereof), while "LS"/"popIII" and "HS"/"Q3" denote respectively light and heavy high-redshift seeds for the black hole population.The predictions are computed by summing the gravitational energy spectra of all the SMBHBs in each model's theoretical population, assuming quasi-circular orbits.As a result, the spectrum has a slope of γ = 13/3 and has no cosmic variance (i.e.we do not account at this stage for the scatter from one realization of the SMBHB population to another).The range shown for each model represents the uncertainty due to the correction for the finite resolution of the SAM's merger tree.In more detail, the lower end of the range represents a model's prediction at the finite resolution, while the upper end is the extrapolationperformed as described in Figure 4 of Klein et al. (2016) -to infinite resolution.The lower arrow (pointing up) should therefore be interpreted as a lower limit, while the upper arrow (pointing down) should be understood as an upper bound (due to the uncertainty of the extrapolation procedure).The extrapolation has not been performed for the model HS-nod-SN-high-accr (B+20), for which we report only the (more robust) prediction at finite resolution.The latter already agrees with the measured amplitude, as a result of a higher accretion rate (by a factor ∼ 4) for SMBHs.
For two of the models in better qualitative agreement with the data (i.e."HS-nod-SN-high-accr (B+20)" and "HS-nod-noSN (B+20)"), we compare the predicted signal with the measured DR2new free spectrum in Fig. 10.Unlike in the case of Fig. 9, these predictions were obtained for multiple specific realizations of the SMBHB population, following Sesana et al. els "HS-nod-SN-high-accr (B+20)" and "HS-nod-noSN (B+20)".The distribution of the predictions represents the scatter among different realizations of the SMBHB population ("cosmic variance").Also shown are power-law fits to the predictions.
(2008) 7 .The probability distribution function plotted in each bin represents the scatter of these multiple realizations, and should therefore be interpreted as a "cosmic variance".Similarly, in Fig. 11 we show the theoretical forecasts for A( f = 1/10yr) from a subset of the models presented above (namely those in qualitative agreement with the data in Fig. 9).These forecasts are obtained by fitting the model predictions (from multiple realizations of the SMBHB population) in the first 9 frequency bins with a power law with γ = 13/3.The error bars represent the 95% confidence regions (accounting for the scatter due to cosmic variance), while the shaded area indicates the 95% confidence region of the posterior for A( f = 1/10yr) (assuming again γ = 13/3).Overall, this qualitative comparison, while somewhat dependent on the details of the SAM implementation, suggests that (i) large delays arising from the dynamics of black hole pairs at large ∼ 100 pc separations are disfavoured, (ii) SMBHB mergers proceed efficiently after their host galaxies have coalesced.Moreover, these results seem to suggest that (iii) accretion onto SMBHs proceeds efficiently, possibly resulting in a larger local SMBH mass function at high masses.

L-Galaxies
Next, we explore the implications that the EPTA results have for L-Galaxies (Henriques et al. 2015;Izquierdo-Villalba et al. 2022), a sophisticated SAM constructed on the dark matter merger trees extracted from the Millennium simulation suite (Springel et al. 2005;Boylan-Kolchin et al. 2009).On top of the galaxy physics, L-Galaxies features a comprehensive modelling for the assembly of SMBHs, including gas accretion triggered by galactic mergers and disc instabilities and dynamical evolution of SMBHBs within the galaxy merger process.The latter accounts for dynamical friction (DF), stellar and gas hardening and, eventually, GW emission.All of these processes are governed by a set of parameters that are tuned to reproduce a vast model (delays increasing left to right) fitting the spectrum in the first 9 frequency bins with γ = 13/3 for multiple realizations of the SMBHB population.The error bars represent the 95% confidence interval for the predictions, and account for the scatter due to cosmic variance.For each model (except for the boosted accretion model HS-nod-SN-high-accr (B+20)), the higher prediction is the extrapolation to infinite SAM resolution, while the lower one is the finite-resolution prediction.The shaded area is the 95% confidence interval for the measurement of A( f = 1/10yr), fixing γ = 13/3.For HS-nod-SN-high-accr (B+20) we only show the result uncorrected for resolution.2022) found that the standard L-Galaxies tuning results in a GWB with log 10 A = −14.9 at f 0 = 1yr −1 , lower than that measured in DR2new.Here we perform a systematic investigation of how the stochastic GWB at nanohertz frequencies depends upon the parameters governing the physics of SMBHs and SMBHBs in the SAM.To this aim, we run L-Galaxies in the following configurations:

Further considerations on the measured spectrum:
eccentricity and statistical biases.
The analyses presented so far give strong indications that the signal is compatible with a cosmic population of SMBHBs swiftly coalescing in the aftermath of galaxy mergers.The relatively flat slope of the measured spectrum might be indicative of strong environmental coupling and high eccentricities, although inference from the data is inconclusive in this respect (see Fig. 8).
The eccentricity of SMBHBs is of particular relevance as it might carry important information on the dynamical processes driving binary evolution at sub-pc scales (see e.g.Roedig & Sesana 2012).While gas-driven dynamics is expected to saturate the binary eccentricity at a value e ≈ 0.4 − 0.6 (Roedig et al. 2011;D'Orazio & Duffell 2021), stellar hardening is known to statistically increase eccentricity without any saturation point (Quinlan 1996), potentially leading to extremely eccentric systems (Sesana 2010).A large binary eccentricity has two important implications for the interpretation of the current data: it flattens the low-frequency spectrum and it speeds up the SMBHB merger process, as inferred by the small τ 0 derived in Sec.3.2.Fig. 13: Orbital parameters (distance between the SMBHs, semi-major axis and eccentricity) of a SMBHB formed in a representative N-body simulation of a galactic merger with parameters drawn from progenitors of likely PTA sources in the IllustrisTNG100-1 cosmological simulation.Mergers are selected from the merger trees of the 100 most massive galaxies at z = 0, based on galaxy mass ratio (major mergers with mass ratio 1 : 4 or higher) and redshift (z ≤ 2).The dashed lines indicate the critical separation a f and the corresponding eccentricity e f at the time in the evolution marking approximately the end of the SMBH inspiral due to DF and the beginning of the hardening phase.Dots represent a and e computed from the apoastron-periastron separation of the two SMBHs before pairing in a bound binary.
Low-redshift massive galaxies are generally gas-poor, and stellar-driven hardening represents the main mechanism driving the evolution of the binaries comprising the bulk of the PTA GW signal.Modelling the whole dynamical evolution, from the first galaxy encounter to black hole pairing and final coalescence, is theoretically and numerically challenging and has been the subject of many studies (e.g.Preto et al. 2011;Khan et al. 2012Khan et al. , 2016;;Nasim et al. 2021;Gualandris et al. 2022).In particular, the binary eccentricity is very sensitive to the initial orbits of the merging galaxies (Gualandris et al. 2022).In ongoing work (Fastidio et al. in prep), we are connecting the sub-pc dynamics of SMBHBs to the large-scale parameters of the galactic encounters extracted from the IllustrisTNG100 simulation (Pillepich et al. 2018).Preliminary results show that mergers of massive galaxies occur preferentially on nearly radial orbits, potentially resulting in highly eccentric binaries.Fig. 13 shows the orbital parameters of a SMBHB formed in a high-accuracy N-body simulation of a representative galactic merger with properties taken from a merger tree in IllustrisTNG100.Merger trees are selected to represent likely PTA sources at low redshifts.The galactic merger is followed from early times through the inspi- ral, pairing and hardening of the SMBHs via a hybrid numerical scheme able to model the evolution self-consistently from kpc to mpc scales (Dehnen 2014).Despite the intrinsic stochasticity of the processes driving binary formation and hardening (Nasim et al. 2020), binaries tend to form with a large eccentricity, which then further grows due to encounters with background stars, as also found by scattering experiments (Sesana 2010).Although more work is needed to determine the distribution of expected binary eccentricities and current EPTA data are not yet strongly informative, this pilot study shows the potential of using these measurements in the near future to constrain the physics of galaxy and SMBH mergers.
When drawing conclusions on the physical properties of the sources of the GW signal, it is useful to bear in mind not only that the constraints on the spectrum are relatively loose (see Fig. 1), but also that the measured parameters can be subject to statistical and systematic biases.To address this, we are conducting an extensive campaign of simulations by injecting and recovering different types of signals in synthetic PTAs mimicking the properties of the EPTA DR2new dataset (Valtolina et al. 2024).We generated individual noises for 25 pulsars using the maximum likelihood values of the measured white noise and drew the red noise and dispersion measure parameters from the posterior distribution of the customised noise models (EPTA and InPTA Collaborations 2023a).We simulated TOAs from multifrequency observations and added a GWB spectrum from an astrophysical population of circular SMBHBs producing a nominal GWB with A 1yr −1 = 2.5 × 10 −15 , consistent with the DR2new measurement at γ = 13/3.We performed 100 experiments by changing the specific noise realization and the sampling of the injected SMBHB population.The analysis was carried using the ENTERPRISE software package (Ellis et al. 2020).
Two illustrative examples of injection-recovery mismatch are shown in Fig. 14.The top panel shows one of the GWB recoveries.Although the injection did not present particular features (e.g.loud CGWs), for this specific noise realization, the recovered signal has a shallow slope with a median value of γ = 3.10.Similar cases have been observed when injecting a pure γ = 13/3 power law with the createGWB function of libstempo (Vallisneri 2020).This shows that even with an ideal setup when simultaneously fitting multiple parameters (102 in this case) in a complex problem, the stochastic properties of the noise can easily bias the recovered signal, especially if the S/N is low (S/N≈ 3.5 for DR2new).Multiple injections with the realistic GWB model and createGWB show systematic biases of recovery of the realistic GWB signal, when modelled with an ideal power law.This is explored in detail in a follow-up work.In the bottom panel of Fig. 14, we show how the presence of some extra highfrequency noise unaccounted for in the MSP noise model can also influence the recovery of the parameters.The setup is the same as in the left panel, but we simulate high-frequency noise mismatch by setting different values of EFAC= 0.8, 1, 1.2 in the recovery.Although the posterior of the signal amplitude is hardly affected, γ can shift significantly depending on whether the highfrequency noise is slightly over-or under-estimated.While these are only two examples, they highlight the complexity of PTA measurements and invite us to be cautious when drawing conclusions that might strongly be influenced by potential biases in the recovered signal.

Implications II: physics of the early Universe
Although a GWB generated by an ensemble of the putative SMBHBs is the most plausible source of the observed commonred noise process in pulsar data, more exotic explanations are possible, such as signals originated in the early Universe.The various possible types of cosmological backgrounds of GWs associated with early Universe physics are reviewed in Caprini & Figueroa (2018) and are found to be stochastic.Similarly to the traditional case invoking SMBHBs, the angular spatial correlation for those scenarios follows the HD curve.However, the spectral shapes of the predicted GW spectra are generally different, which can help to disentangle between different types of backgrounds.In this work, we focus on four possible scenarios: 1. an inflationary GWB from the amplification of quantum fluctuations of the gravitational field, 2. a GWB from a network of cosmic string loops, 3. a GWB from vortical (M)HD turbulence at the QCD energy scale, 4. a scalar-induced GWB arising from inflationary scalar perturbations at the 2nd order in perturbation theory.
Given the low significance of the detected signal and the limited number of probed frequency bins due to the short timespan of the data, one cannot currently perform a reliable model selection.Therefore, throughout the section, we consider these scenarios separately and assume that each of them can fully explain the detected signal independently.Analysis invoking more complex models with simultaneous fits for multiple scenarios as well as opportunities to disentangle between those (e.g.Goncharov et al. 2022;Kaiser et al. 2022) will be considered in a number of future works.

Implications on a stochastic background of primordial (inflationary) gravitational waves
Here we address the GWB possibly generated during inflation (Grishchuk 1975;Rubakov et al. 1982;Starobinskii 1985;Fabbri & Pollock 1983;Abbott & Wise 1984).In the standard inflationary scenario, tensor quantum vacuum fluctuations of the metric are amplified by the accelerated expansion, leading to a GWB as they subsequently re-enter the horizon during the radiation (or matter) era.The cosmic microwave background (CMB) and Big Bang Nucleosynthesis (BBN) provide precise measurements of the radiation energy density, from which one can derive weak upper bounds on the amplitude of such a GWB (see e.g.Caprini & Figueroa 2018, and references therein).Furthermore, tensor metric perturbations lead to CMB temperature anisotropies and polarisation at large angular scales (Sachs & Wolfe 1967;Starobinskii 1985;Kosowsky 1996;Allen & Koranda 1994).Since the anisotropies and polarisation detected so far are due to scalar perturbations, it is possible to constrain the energy density of a GW background by placing an upper limit on the tensor-to-scalar ratio r at CMB scales: recent upper bounds are given e.Within the slow-roll consistency relation, the GWB spectral slope is therefore slightly red-tilted, causing this signal to be out of reach of most current and planned GW detectors such as PTAs, LISA, aLIGO, aVirgo or the Einstein Telescope.On the other hand, it is fair to consider that the spectral slope could vary over the large frequency span ranging from CMB scales to those probed by GW detectors (Lasky et al. 2016).Inflationary scenarios breaking the consistency relation at small scales might therefore produce a detectable GWB, if they lead to a blue-tilted spectrum.In this case, PTAs, LISA and ground-based devices can place upper limits on n T (see e.g.Abbott et al. 2017).It is interesting to investigate which kinds of processes could give rise to a blue-tilted GW background while still obeying CMB constraints at large scales.One possibility is the presence, after inflation, of a stiff component in the Universe, with an equation of state w > 1/3 (Boyle & Buonanno 2008;Giovannini 1998;Boyle & Steinhardt 2008).The enhancement of the tensor spectra can also be produced during inflation thanks to processes such as, for example, (i) particle production during inflation (see e.g.Sorbo (2011); Barnaby et al. (2012); Cook & Sorbo (2013); Anber & Sorbo (2012); Bartolo et al. ( 2016)) (ii) enhancement of tensor perturbations for example by spectator fields, or spacedependent inflation (see e.g.Bartolo et al. (2007); Biagetti et al. (2013Biagetti et al. ( , 2015)); Fujita et al. ( 2015)) (iii) modified gravity theories such as f (R) or Horndeski gravity (Horndeski (1974); Sotiriou & Faraoni (2010)) and iv) enhanced scalar perturbations at small scales and/or primordial black holes, which are treated in Sec.4.4.

Analysis
Similarly to what was done in Lasky et al. ( 2016), we constrain the key parameters defining the GWB, r and n T , while being agnostic on the underlying process generating the bluetilted spectrum.If we assume that the common quadrupolar red noise signal present in EPTA data is of inflationary origin, these two parameters can be estimated using the DR2new dataset.Note that the spectral index n T is expected to vary with the frequency scale considered (see e.g.Giarè & Melchiorri 2021;Giarè et al. 2023;Auclair & Ringeval 2022).However, for simplicity, here we consider a constant n T all the way from CMB scales to those corresponding to the (narrow) EPTA frequency band.It is then possible to approximate the fractional characteristic GW energy density using (Lasky et al. 2016;Caprini & Figueroa 2018) where the second line is valid in the PTA frequency band, and has been obtained by setting h 2 Ω rad = 2.47 • 10 −5 with h = 0.67, the amplitude of the scalar spectrum P * R = 2•10 −9 , and f * ≈ 7.7× 10 −17 Hz related to the CMB pivot scale k * = 0.05/Mpc (Planck Collaboration et al. 2014).f eq denotes the frequency entering the horizon at matter-radiation equality.
We then use the nine lowest frequency posteriors of the RMS free spectrum shown in Fig. 1 (see Moore & Vecchio 2021;Lamb et al. 2023;Leclere et al. 2023, for details on the method) to fit the inflationary spectrum of Eq. 14 and obtain posteriors on log 10 r and n T .Results are reported in Fig. 15.Note that, since γ = 5 − n T , the correlation between the amplitude and spectral index of the signal is compatible with Fig. 5.The 90% credible (symmetric) intervals are log 10 r = −12.18+8.81 −7.00 and n T = 2.29 +0.87  −1.11 .The obtained value of n T corresponds to a PSD spectral index of γ ≃ 2.7, as in Fig. 1.The excessively small value of r is a consequence of the simplistic parameterisation of Eq. 14, which assumes a constant n T at all scales.The fractional energy density spectrum obtained from the maximum a posteriori parameter values is plotted in Fig. 17.
We have so far considered a primordial background to be the only source of GWs in our data.We now recall that the most plausible and loud source of a GW background at these frequencies remains that of a SMBHB background.It is therefore likely that any signature for a cosmological background needs to be considered in parallel with a SMBHB background, or in this case more accurately termed 'foreground'.Kaiser et al. (2022) have explored the likelihood of detecting a cosmological background in the presence of a SMBHB foreground using simulations, and found that the shallower the slope of the cosmological background (for example γ = 4 as opposed to γ = 5), the harder it is to detect (and the longer it takes, possibly more than 20 years).According to these simulations, this does not bode well for an even shallower slope, similar to the one detected in DR2new with a possible γ < 3.
Here we explore a superposition of these two backgrounds in the DR2new dataset.Considering a two-component GWB for the common red noise model, we place constraints on log 10 r for given values of n T spanning the range [−1, 3].In this case, our null hypothesis is a GWB from a population of GW-driven circular SMBHBs parameterised only by the PSD amplitude log 10 A of Eq. (3) (we fix γ = 13/3).We run several analyses with a fixed n T for the inflationary background, sampling over (log 10 r, log 10 A).For each of the n T values, we obtain a distribution for log 10 r and take the 95% quantile as an upper bound.As found in Lasky et al. (2016), n T and the log 10 r upper bounds are related with good precision by a linear relation: Our analysis gives a = −0.16 and b = 0.70, which is comparable to the forecast values given in Lasky et al. ( 2016) (note that they normalise r to 0.11).

Discussion
From the analysis of the DR2new dataset above, we have obtained credible intervals for the tensor-to-scalar ratio r and the spectral index n T .This was performed assuming that reheating is instantaneous, and that inflation is followed directly by the radiation-dominated era, for which the equation of state parameter of the Universe is w = 1/3.Under this assumption, one finds that the best-fit value for the tensor spectral index is n T = 5 − γ ≃ 2.3, which is directly linked to the best-fit PSD spectral index γ ≃ 2.7.This high value of n T is not consistent with slow roll inflation.However, if inflation is followed by a stage in which w 1/3, the relation between the PSD spectral index γ and the primordial tensor spectral index n T changes to (Arzoumanian et al. 2016;Caprini & Figueroa 2018) again with γ ≃ 2.7.If a stiff fluid component (w > 1/3) were to dominate the Universe for a finite amount of time after inflation, the last term in Eq. 16 would be bounded between 0 and −2.
Hence, n T ≳ 0.3, meaning that even allowing for the presence of a stiff component after inflation, it does not seem possible to explain the common red noise in the context of slow roll inflation (n T ≃ 0) for the best fit value γ ≃ 2.7.However, by broadening the range of possible values to γ ≥ 3, n T ≃ 0 does become compatible with the common red noise.

Implications on a background of cosmic strings
Cosmic strings are line-like topological defects that may form after a symmetry-breaking phase transition in the early Universe (Kibble 1976;Hindmarsh & Kibble 1995); they are generic predictions of most Grand Unification Theories scenarios (Jeannerot et al. 2003).These 1D objects are characterised by the string tension Gµ (or equivalently their energy per unit length) which is related to the energy scale of the phase transition.
Overall, cosmic strings combine relativistic velocities and large energy densities, making them natural sources of GWs.These GWs may take the form of bursts from cusps, kinks and kink-kink collisions on the loops (Damour & Vilenkin 2001), and have been searched for in the LIGO/Virgo/KAGRA (LVK) detectors (Abbott et al. 2018(Abbott et al. , 2021)).Cusps are points on the string which, in the Nambu-Goto approximation, propagate at the speed of light, and the string doubles back on itself.On the other hand, kinks are discontinuities in the tangent vector of a string and are formed at each intercommutation of strings.The future space-based detector LISA will be sensitive to cosmic string bursts for tensions as low as Gµ ∼ 10 −11 (Auclair et al. 2023a).Most notably, in the event that LISA detects GW bursts from cosmic strings, they will likely repeat multiple times a year due to the periodic nature of the cosmic string loops (Auclair et al. 2023b).
The incoherent sum of these GW bursts also produces a stochastic GWB (SGWB), which has been looked for by LVK (Abbott et al. 2018(Abbott et al. , 2021)).LISA is expected to be able to detect a SGWB from cosmic strings with tensions as low as Gµ ∼ 10 −17 (Auclair et al. 2020).The strings background has already been looked for in EPTA (Sanidas et al. 2012;Leclere et al. 2023), in NANOGrav (Blasi et al. 2021;Ellis & Lewicki 2021) and PPTA (Bian et al. 2022;Chen et al. 2022).These earlier analyses consistently found a good agreement between the common uncorrelated red signal (CURN) present in the data and a SGWB from cosmic strings with tensions Gµ ∼ 10 −10 − 10 −11 .

Description of the models
Since the SGWB is sourced by sub-horizon cosmic string loops, the central quantity in our analysis is the loop density distribution, n(ℓ, t), where ℓ is the invariant length of the loop (defined by its total energy divided by µ) and t is cosmic time.We focus on the most up-to-date loop distributions, calibrated on large scales with Nambu-Goto simulations (Lorenz et al. 2010;Blanco-Pillado et al. 2014).Loops are formed from the intercommutation of super-Horizon infinite strings, with a maximal size ℓ ≈ 0.1t at time t.(Note that the simulations assume an intercommutation probability of 1, as is the case for field theory strings.)The two models -which we refer to as BOS (Blanco-Pillado et al. 2015) and LRS (Lorenz et al. 2010), after the names of their authors -differ in their treatment of small loops, particularly on scales at which gravitational radiation (not included in numerical simulations) can have important effects.Indeed, compared to the BOS model, the LRS model has an additional population of small loops, which leaves an imprint in the SGWB spectrum (Auclair 2020).The explicit expressions for the two loop distributions are given in Abbott et al. (2018), and both have been considered by the LVK (Abbott et al. 2021) and LISA (Auclair et al. 2020) collaborations.The third observing run of LVK placed constraints on Gµ, based on the non-detection of a SGWB.These are Gµ ≲ 10 −8 for BOS and Gµ ≲ 10 −15 for LRS.At the frequency of ground-based detectors, the SGWB signal is produced by loops formed during the radiation era.At low PTA frequencies, the SGWB signal is dominated by larger loops, namely those formed at recent times, in transition from the radiation to matter era and also in the matter era (Ringeval & Suyama 2017;Auclair 2020;Leclere et al. 2023).
Other than Gµ, a further parameter appears in n(ℓ, t), namely Γ, which describes the rate at which loops lose energy through gravitational radiation: l = −ΓGµ.If Γ is large, fewer loops are present in the distribution since loops decay more rapidly.On the other hand, GWs are emitted more intensively: the final effect is a combination of these two, which impacts the SGWB of the BOS and LRS models in subtly different ways (Auclair 2020).The value of Γ depends on the properties of loops, and in particular on how many cusps, kinks, and kink-kink they contain as (Damour & Vilenkin 2001;Siemens et al. 2007;Ringeval & Suyama 2017) where N c,k is the number of cusps/kink events per oscillation period of the loop, N 2 k /4 the number of kink-kink collisions, and where g 2 = √ 3/4, g c 1 ≈ 0.85, g k 1 ≈ 0.29, g kk 1 ≈ 0.1.In this paper, we consider 2 cases.For the first model (i) N c = 2 and N k = 0, for which Γ = 57 (Leclere et al. 2023), a value close to that observed in numerical simulations of individual loops (Vachaspati & Vilenkin 1985;Blanco-Pillado & Olum 2017).Therefore, the only free parameter for this model is Gµ.As for the second model (ii), we fix N c = 1 and allow N k to vary between 1 and 200 (with a uniform prior) so as to account for theoretical uncertainties on the initial number of kinks at loop creation, and on the efficiency of the gravitational backreaction that should smooth out the loops (see the LVK analysis, Abbott et al. 2021, for a similar approach).Therefore, this is a two-parameter model (Gµ, N k ).
The fractional energy density of the SGWB per logarithmic interval of frequency is given by where H 0 is the Hubble constant, H(z) is the Hubble parameter for which we assume standard ΛCDM cosmology with the Planck-2018 fiducial parameters (Aghanim et al. 2020), and t 0 is the cosmic time today.The sum is over the cusp, kink and kink-kink contributions (labelled with the index b) for which q b = 4/3, 5/3 and 2 respectively.

Analysis results
Our analysis follows the one presented in Leclere et al. (2023), which was based on six pulsars with the best timing precision from the EPTA early Data Release 2 (Chen et al. 2021).We now analyse the DR2new dataset with 25 pulsars, using the (computationally heavy) hierarchical likelihood data analysis method described in PaperIII.We sample the SGWB parameters (Gµ) or (Gµ, N k ) as well as the individual pulsar noise parameters.
Solid lines in Fig. 16 show the posterior distribution on log 10 (Gµ) in case (i), for which N k = 0 and Γ = 57.The string tension 90% credible (symmetric) interval is log 10 Gµ = −10.07+0.47  −0.36 (resp.−10.63 +0.24 −0.22 ) for the BOS (resp.LRS) model.The corresponding SGWB spectra are shown in Fig. 17 where, for each model, we set Gµ to their median values.We also consider the two-component SGWB model made of the sum of a background from cosmic strings and one from a population of GW-driven circular SMBHBs characterised by the PSD of Eq. ( 3).The posteriors of the two background parameters (log 10 A, log 10 Gµ) are highly correlated, since both provide a possible explanation for the detected signal.As a result, the posterior on log 10 Gµ no longer has compact support, but a tail to lower values (see the dashed lines in Fig. 17).We therefore extract the 95-quantile of the string tension posterior to obtain an upper bound of log 10 Gµ < −9.77 (resp.−10.44) for the BOS (resp.LRS) models.
The DR2new dataset exhibits a shallower slope for the PSD of the common red signal than DR2full.While cosmic strings were a good fit to the common red signal of 6 pulsars of DR2full (Leclere et al. 2023), this is no longer true for DR2new.This is because the predicted SGWB PSD is generally steeper than the measured correlated red signal in the data, as can be seen in Fig. 17.
For case (ii), we obtain very similar results to those discussed in Leclere et al. (2023).Namely, we obtain quasi noninformative posteriors for N k , showing that the data can be equally explained by a population of kinky loops with N k ≳ 120.In other words, we cannot extract any upper bound on the number of kinks, since this quantity is degenerate with Gµ.

Implications on background from turbulence around the QCD energy scale
Turbulence can arise in the early Universe in the aftermath of a first-order phase transition (Witten 1984;Kamionkowski et al. 1994), or can be driven by pre-existing primordial magnetic fields (Quashnock et al. 1989;Brandenburg et al. 1996).If the (magneto-)hydrodynamic turbulence were present around the QCD epoch, when the Universe had a temperature of T * ∼ 100 MeV, it would generate a GWB in the PTA band.The characteristic scale of the turbulence, determining the characteristic GW frequency, is in fact related to the (comoving) Hubble radius at that epoch λ * ≃ O(H −1 * ), where and g * denotes the number of relativistic degrees of freedom.If a large lepton asymmetry and/or primordial magnetic fields were present in the early Universe, the QCD phase transition might have been of first order (Schwarz & Stuke 2009;Wygas et al. 2018;Middeldorf-Wygas et al. 2020, 2022;Vovchenko et al. 2021;Cao 2023).In this case, one would expect additional sources of GWs, from the collision of broken phase bubbles and the subsequent development of sound waves in the primordial fluid (Kosowsky et al. 1992;Kosowsky & Turner 1993;Caprini et al. 2008;Huber & Konstandin 2008;Jinno & Takimoto 2017;Cutting et al. 2018;Hindmarsh et al. 2014Hindmarsh et al. , 2015Hindmarsh et al. , 2017)).This was analysed for PTAs, e.g. in Moore & Vecchio (2021); Arzoumanian et al. ( 2021); Xue et al. (2021).In what follows, we focus on the GWB generated by decaying (M)HD turbulence.

Description of the model
The presence of bulk velocity and magnetic fields produce anisotropic stresses, which in turn act as a source of GWs (Kamionkowski et al. 1994;Kosowsky et al. 2002;Dolgov et al. 2002;Caprini & Durrer 2006;Gogoberidze et al. 2007;Caprini et al. 2009).This has been recently studied via numerical simulations in Roper Pol et al. (2020a,b); Brandenburg et al. (2021); Roper Pol et al. (2022a).In particular, Roper Pol et al. (2022a) show that the envelope of the GWB produced by decaying MHD turbulence can be estimated analytically, assuming that the anisotropic stresses from the velocity and magnetic fields vary more slowly than the dynamical production of GWs.This was also validated by numerical simulations of purely kinetic turbulence in Auclair et al. (2022).This assumption leads to the following GWB signal: where Ω * is the ratio of the (M)HD turbulent energy density to the radiation one, and λ * H * is the ratio of the characteristic length scale of the turbulence, λ * , to the comoving Hubble horizon H −1 * at the QCD epoch.The parameter A ≃ 1.75 × 10 −3 is the efficiency of GW production,8 estimated in Roper Pol et al. (2022a).The function F GW,0 is the fractional radiation energy density at the epoch of GW generation to its value at the present time.It depends on the temperature scale T * via the number of degrees of freedom g * , The spectral shape of the GWB signal, S turb ( f ), is where B ≃ 50 λ * H * −2 is a normalising factor, and δt fin denotes the effective duration of the turbulence.The latter can be estimated, from the numerical simulations performed in Roper Pol et al. (2022a), to be δt fin ≃ 2λ * / √ 1.5 Ω * .The function p Π (λ * f ) in Eq. ( 23) denotes the spectrum of the anisotropic stresses.For solenoidal fields (e.g. a primordial magnetic field or vortical bulk fluid motion) characterised by a typical correlation scale of the order of the turbulence scale λ * , it is constant for f < 1/λ * .Furthermore, it decays as f −11/3 for f ≳ 1/λ * , if the turbulence is of the Kolmogorov type, as we assume here.Hence, the resulting spectral shape of the GWB in Eq. ( 21) presents three power laws: f 3 at frequencies below the inverse effective duration of the turbulence f < 1/δt fin , f at intermediate frequencies 1/δt fin < f < 1/λ * , and f −8/3 at large frequencies f > 1/λ * .
The GWB produced from vortical (M)HD turbulence is therefore determined by three parameters: the temperature scale T * , the turbulence strength Ω * , and the turbulence characteristic length scale λ * H * .By causality, λ * H * is bound to be smaller than one.In general, also Ω * ≲ 1, otherwise turbulence would change the dynamics of the Universe.However, note that the template described above has been validated in principle only for non-relativistic plasma motions, for which Ω * ≲ O(0.1).

Analysis results
As in subsection 4.1, here we use the fast free spectrum analysis method on DR2new data to constrain the model, considering the nine first Fourier bins of the RMS spectrum of For values of Ω * below 0.1, the model can only explain the level of correlated noise at the lowest frequency bin if the amplitude of the spectrum is sufficiently high.This can be achieved only if λ * H * is close to 1 and the peak frequency lies within the PTA frequency range, implying T * ∼ 60 MeV.However, at frequencies around the peak, the signal corresponds to a power spectral density for the residuals steeper than γ ∼ 4, which cannot fit the data well.For this reason, values of Ω * ≲ 0.1 are disfavoured.
For larger values of Ω * , the f 3 part of the spectrum at frequencies below δt −1 fin ∝ √ Ω * /(λ * H * ) × H * (T * ) can enter the PTA band with a sufficiently high amplitude.Furthermore, the distance between the break at δt −1 fin and the spectral peak at 1/λ * becomes minimal in the limit Ω * ∼ 1.Both of these characteristics lead to a better fit to the data.This is recovered in the posteriors of Figure 18, together with the degeneracy between λ * H * and Ω * from the signal amplitude (see Equation 21), and the degeneracy between λ * H * and T * from the break at 1/δt fin (note that the dependence of the latter on √ Ω * is subdominant).The model therefore provides a good fit to the data in the limit of large Ω * , close to the upper bound of the prior.The extension of the dataset to longer observation time will be crucial for further constraining this model at low frequencies.

Implications on the 2nd-order GWB produced by primordial curvature perturbations
It is well-known that scalar, vector and tensor modes of the perturbed metric do not mix at linear order of the Einstein equations (Lifshitz 1946;Baumann 2022).However, scalar curvature perturbations will source propagating tensorial modes (GWs) at the Fig. 18: 2D posteriors for the parameters of the background from turbulence around the QCD energy scale obtained using a free spectrum fit on DR2new data.The 68% and 95% credible regions are displayed.
2nd order in perturbation theory (Tomita 1967;Matarrese et al. 1993Matarrese et al. , 1998;;Noh & Hwang 2004;Carbone & Matarrese 2005;Ananda et al. 2007;Baumann et al. 2007).Such scalar curvature perturbations and associated primordial density fluctuations inevitably exist in the Universe and can be directly constrained by observations of the CMB.The latest Planck data (Planck Collaboration et al. 2020b) suggests that the power spectrum of the curvature perturbations is nearly scale-invariant with the amplitude A ζ ∼ 2 × 10 −9 , which implies a marginal energy-density of the generated GWB.Specifically, when projected to the PTA sensitivity band, the fractional contribution of the energy density in the associated GWs becomes Ω gw ∼ 10 −24 , which is practically non-detectable by current experiments.On the other hand, some models of inflation (see, for example, Di & Gong 2018;Byrnes et al. 2019;Braglia et al. 2020;Yi & Fei 2023, and references therein) make it possible to produce a sharp increase in the power spectrum of primordial curvature perturbations over many orders of magnitude at small scales.While the CMB is only capable of directly sampling large cosmological scales with k ∼ 10 −3 − 10 −1 Mpc −1 , small scales stay largely uncovered.PTAs provide a unique opportunity to complement the CMB measurements by indirectly probing the scalar curvature perturbations in a scale range k ∼ 10 6 − 10 8 Mpc −1 through the second-order generated GWB, and to place bounds on the steepest possible growth of the power spectrum as well as corresponding models of inflation (Saito & Yokoyama 2009;Bugaev & Klimai 2011;Chen et al. 2020;Dandoy et al. 2023;Zhao & Wang 2023).
In this work, we consider two models of the primordial curvature power spectrum: i) monochromatic and ii) powerlaw.For the former, the primordial spectrum is modelled as: where A ζ is a dimensionless amplitude and k * is a wavenumber at which the monochromatic power spectrum has a Dirac-delta Article number, page 18 of 30 EPTA+InPTA: GWB Interpretation The posterior distribution is overlaid with the current constraints on the primordial power spectrum using Planck data (CMB).The grey colour depicts the 2-σ-confidence intervals.The purple shaded area represents the bounds from spectral distortions (Chluba et al. 2012).For comparison in green we place the prediction of the primordial spectrum of scalar perturbations in the two-field model of inflation described in Braglia et al. (2020) for a range of the model parameters.All three models result in PBH mass functions peaked at ∼ 35 M ⊙ with the brightest line corresponding to the dark matter fraction of PBHs of ∼ 0.01.25).Left panel: 1σ and 2σ contours of the posterior distributions on the amplitude A ζ and the slope of the power spectrum n s obtained by the analysis of DR2New.Right panel: 1σ and 2σ contours of the power spectra inferred from the DR2New analysis by picking 1000 random samples from the posteriors overlaid with the current constraints on the primordial power spectrum using the latest Planck data.The grey colour depicts the 2σ-confidence intervals.The green lines and purple shaded areas are the same as in Fig. 19.peak.For the second case: where n s characterises the slope and k 10yr is a normalizing scale k 10yr = 2π/(10yr × c), so that A 10yr ζ corresponds to dimensionless amplitude at ten years.
In the first scenario, a semi-analytical solution for the induced spectrum of GWB exists and is given by (Kohri & Terada 2018;Espinosa et al. 2018): where k = k/k * and Θ is the Heaviside theta function.In spite of being nonphysical, the δ-function peak approximately describes the maximum of the produced GWB in the inflationary model with the steepest possible k 4 growth of a spectral peak in the single-field inflation at small scales (see Figure 7 in Byrnes et al. 2019).
In the second case of a more general (and more realistic) power-law spectrum, the result can only be obtained numerically (Kohri & Terada 2018): where Q(n s ) is the scaling factor which can be evaluated in a range of n s using interpolation points from Table 1 of Kohri & Terada (2018).After its production, the GWB is damped due to quantum interactions with the particles of the primordial plasma at the radiation-dominated epoch, and redshifted inversely proportionally to the scale factor (as it also occurs to radiation) starting from the epoch of matter-radiation equality (Saikawa & Shirai 2018).The present value of the fractional energy density is then: where T is the temperature of the Universe at the moment when structures of a typical size 1/k re-enter the horizon 9 , T eq is the temperature of the Universe at the epoch of matter-radiation equality, g * and g * s are relativistic degrees of freedom and degrees of freedom in entropy, respectively.The final expression for the auto-power spectral density of the timing residuals is: where H 0 is the Hubble constant at the present epoch.
The outlined formalism was applied to the DR2new version of the latest EPTA dataset.The number of frequency components which was used for the Fourier representation of the signal was fixed to 9. We have chosen broad uninformative priors for the parameters: uniform in [−6, 3] for log 10 A ζ and log 10 A 10yr ζ , uniform in [4,12] for log 10 (k * /Mpc −1 ), and uniform in [0.4,2.6] for n s .Boundaries for the latter are constrained by the limitations of the numerical approximation of the power law model.For this analysis, we assumed that the common red noise process detected in the latest EPTA dataset can be fully explained by the 2nd-order scalar-induced GWs.
Results for monochromatic and power law models are shown in Figures 19 and 20, respectively.The 2D posterior distribution of the model parameters of the monochromatic model is depicted with orange contours on the right panel of Fig. 19.One may notice that the regions of the highest probability are strongly elongated due to a strong positive correlation between A ζ and k * ; these parameters are essentially degenerate.Therefore, DR2new can only provide lower limits on the characteristic scale and amplitude of the monochromatic model: log 10 (k * ,0.05 /Mpc −1 ) = 7.6 and log 10 A 0.05 ζ = −1.7 meaning that a whole range of models predicting Dirac-delta power spectrum can equally good describe the signal of DR2new.This behaviour is explained in the left panel of Fig. 19, for which we have simulated GWB signal generated by the monochromatic primordial scalar perturbations, Eq. ( 26), and attempted to recover a more general power law model of the form S ( f ) ∼ f −γ used in PaperIII to model an arbitrary common red noise process.After a rapid decrease, the recovered slope stabilises at γ ∼ 3.5 in the limit of large k * , and becomes consistent with both values obtained with DR2new and the theoretically predicted 13/3 for the background from circular, GW-driven SMBHBs.This degeneracy can raise important issues when one tries to disentangle one signal from another.For the power-law case, the 2D posteriors are shown in the left panel of Fig. 20 with the following means and 1-σ uncertainties: log 10 A ζ = −2.94+0.42  −0.46 and n s = 2.11 +0.25 −0.32 .On the right panels of Figs.19 and 20, we also overplot the inferred primordial power spectrum with the one obtained from the Planck data.The orange areas on the right panel of Fig. 20 are 1σ and 2σ contours of the power-law model described by Eq. 25 reconstructed using 1000 random draws from the posteriors.To explain the observed signatures of the DR2new in terms of the second-order GWB from the primordial scalar perturbations, an excess in the primordial spectrum at low scales should be invoked without violating the CMB inflationary parameters.Such excess has been proposed in many papers, e.g. in the aforementioned works by Ivanov et al. (1994); Germani & Prokopec (2017); Di & Gong (2018); Cai et al. (2018); Biagetti et al. (2018); Byrnes et al. (2019); Motohashi et al. (2020); Braglia et al. (2020); Yi & Fei (2023).Notably, this power excess would lead to a copious production of primordial black holes (PBHs) at the radiation-dominated stage, which is sometimes taken by the authors as the motivation to introduce them as cold dark matter candidates (Carr et al. 2016, and reference therein).The PBH formation from cosmological perturbations has been extensively explored Sasaki et al. (2018).On the radiation-dominated stage, the PBH mass is related to the mass inside the horizon at the time of the perturbation entering, , where m Pl is the Planckian mass, T is the temperature.In terms of the mode comoving wavenumber at the moment of the horizon crossing, k = aH, the part of the horizon mass collapsing into a PBH reads For example, in the two-field inflationary model by Braglia et al. (2020), a peak around k ∼ 2 × 10 6 Mpc −1 could explain the power A ζ ∼ 10 −2 and simultaneously lead to the production of PBHs with masses peaked at ∼ 35M ⊙ (see also the analysis of the NANOGrav results in Vaskonen & Veermäe 2021).Interestingly, such a peak seems to be observed in the chirp mass distribution of LVK merging binary black holes (Abbott et al. 2023).A more general list of primordial power spectra, as well as a careful retranslation of them to the PBH abundance and their mass function, will be considered in a follow-up paper (Porayko et al., in prep.).

Implications III: dark matter
Unlike spatially and temporally-correlated stochastic processes discussed in Secs. 3 and 4, in this Section, we explore a possible deterministic contribution to the EPTA signal from ultralight scalar-field dark matter (ULDM).For comparison, the morphology of a putative ULDM signal, predicted by Khmelnitsky & Rubakov (2014), is similar to a CGW from a SMBHB, that is, it is prominent only in one frequency bin.Given that the signal observed with DR2new is mostly apparent in the first two fundamental T −1 frequency bins, it is of interest to consider possible contributions from physical processes with narrowband spectra.Therefore, the analysis presented here complements the CGW interpretation of the signal by EPTA and InPTA Collaborations (2023c).Additionally, the ULDM search with DR2full is performed in Smarra et al. (2023).Dark matter currently constitutes approximately 26% of the energy density of the Universe, as confirmed by, e.g., galactic rotation curves (Rubin & Ford 1970;Rubin et al. 1980;de Salas et al. 2019), baryonic acoustic oscillations and cosmic microwave background measurements (Bennett et al. 2013;Aghanim et al. 2020) as well as galaxy surveys (Escudero et al. 2015).The standard cold dark matter (CDM) picture, whose leading candidates are the Weakly Interacting Massive Particles (WIMPs) (Arcadi et al. 2018) and the QCD axion (Luzio et al. 2020), successfully grasps the large-scale structure of the Universe.However, it presents some well-known issues, when it comes to explaining observations at scales smaller than O(kpc).Among these, the cusp-core problem (Flores & Primack 1994;Moore 1994; Karukes, E. V. et al. 2015) concerns the inconsistency between the observation of a flat density profile in the centre of galaxies and the power-law-like behavior predicated by CDM, while the mismatch between the simulated and observed number of dwarf galaxies in the proximity of our Milky Way is often referred to as the missing satellite problem (Klypin et al. 1999;Moore 1994).
While a thorough understanding of baryonic physics feedback mechanisms (Navarro et al. 1996;Governato et al. 2012;Brooks et al. 2013;Chan et al. 2015;Oñorbe et al. 2015;Read et al. 2016;Morganti 2017) might help to alleviate some of these issues, they can be more easily disposed of assuming that DM is an ultralight (m ϕ ∼ 10 −22 eV) scalar/pseudoscalar or axionlike field, whose astrophysically large (O(kpc)) de Broglie wavelength suppresses power on small scales.Moreover, ultralight scalars are also generally present in string theory compactifications (Green et al. 1988;Svrcek & Witten 2006;Arvanitaki et al. 2010), which makes them interesting candidates for new physics as well.CMB-based arguments are used to constrain m ϕ ≳ 10 −24 eV (Hlozek et al. 2015), while Lyman-α bounds push the limit up to m ϕ ≳ 10 −21 eV, provided that the ultralight particles account for more than 30% of the full dark matter budget (Iršič et al. 2017;Armengaud et al. 2017;Kobayashi et al. 2017;Nori et al. 2018;Rogers & Peiris 2021).A lower limit of m ϕ ∼ 10 −19 eV is claimed by studies of ultra-faint dwarf (UFD) galaxies (Hayashi et al. 2021;Dalal & Kravtsov 2022), but a wide consensus is yet to be reached.In a seminal work, Khmelnitsky & Rubakov (2014) showed that the travel time of pulsar radio beams is affected by the gravitational potential induced by ULDM particles, making thus PTAs excellent facilities to investigate the existence of ULDM particles.Moreover, they represent complementary probes which do not suffer from the small-scale structure modelling uncertainties that affect non-CMB bounds (Schive et al. 2014;Zhang et al. 2019), as the aforementioned Lyman-α or UFD limits.In the following, we robustly assume that ULDM interacts only gravitationally, therefore giving rise to periodic oscillations in the TOAs of radio pulses as described in Khmelnitsky & Rubakov (2014).
Being non-relativistic and with a very large characteristic occupation number, the ULDM scalar field can be described as a classical wave whose oscillation frequency is twice its mass m ϕ (Khmelnitsky & Rubakov 2014).The periodic displacement that ULDM induces on the TOAs of signal from a pulsar P can be written as follows (Porayko et al. 2018): where where γ P (γ E ) is a pulsar (Earth) dependent phase and φ2 P ( φ2 E ) takes into account the stochastic nature of the axion-like field near the pulsar (Earth).The parameters and their priors are summarised in Table 1.Considering a typical value of v ϕ ∼ 10 −3 for the ULDM velocity, the region in which the scalar field oscillates coherently, i.e. with the same value of φ, is spanned by the coherence length: In particular, notice that φ2 E and φ2 P should be taken as: different parameters when the average pulsar-Earth and pulsar-pulsar distance is larger than the coherence length.the same parameter when the average pulsar-Earth and pulsar-pulsar distance is smaller than the coherence length.
Following the procedure in Smarra et al. (2023), we analyze three separate cases, which we refer to as the uncorrelated, the pulsar correlated and the correlated limit.As the average inter-pulsar and Earth-pulsar separation is ∼ kpc, the correlated and uncorrelated scenarios stand out as exact limits at the low mass and high mass end of the PTAs band, respectively.Instead, the pulsar correlated limit holds when the coherence length of ULDM is smaller than the Galacto-centric radius probed by rotation curves (inner ∼ 20 kpc), but larger than the average interpulsar and pulsar-Earth distance.More specifically, the correlated regime holds for masses lower than m ϕ ∼ 2 × 10 −24 eV; the pulsar correlated regime for 2 × 10 −24 eV ≲ m ϕ ≲ 5 × 10 −23 eV and the uncorrelated limit for m ϕ ≳ 5 × 10 −23 eV.We defer a more detailed study to future analysis.Based on the above, we fit the model from Eq. ( 31) to the DR2new as in Smarra et al. (2023).As an example, marginalised posterior distributions for Ψ c and m ϕ are plotted in Fig. 21.All the three limits peak at a ULDM mass m ϕ ∼ 10 −23 eV with amplitude Ψ c ∼ 10 −15 , which translates into a density ρ ϕ ≃ 0.6 GeV/cm 3 of the scalar field.While this value is higher than the fiducial local DM density ρ DM ≈ 0.4GeV/cm 3 , it is well within the measurement uncertainties (Bovy & Tremaine 2012;Read 2014;Sivertsson et al. 2018;de Salas 2020).
Additionally, we search for a potential ULDM signature alongside the SMBHB gravitational-wave background.Thus, we introduce a new model that contains ULDM contributions alongside a common red signal to account for gravitational wave contributions.We find no ULDM signal under this hypothesis, in   2018) for further details.The blue, orange and brown lines represent the 95% Bayesian upper limit on Ψ c obtained from the EPTA DR2new dataset with the correlated, uncorrelated and pulsar correlated analysis, respectively.The purple line shows the expected ULDM abundance computed from Eq. ( 32).
we put 95% upper limit on the amplitude Ψ c and the density ρ ϕ of the scalar field.2023) pushes these limits down, thus implying that such ULDM candidates cannot constitute the entire DM.We highlight that our bounds extend below the DR2new sampling frequency f = 1/T ≈ 3 • 10 −9 Hz, with T ≈ 10yr.In fact, while it might naively be thought that the ULDM field needs to complete an oscillation in the timescale τ PTA of the PTA experiment to produce a detectable effect, we point out that an ULDM wave with frequency m ϕ < 1/τ PTA can still be approximated by an expansion in powers of m ϕ t (Kaplan et al. 2022), though the sensitivity of PTAs will be reduced due to degeneracies with the timing model (Ramani et al. 2020).Relying on the robust CMB bounds mentioned before, we fix the lower end of our search at m ϕ = 10 −24 eV.Importantly, we remove the φE = φP parameter from the search in the correlated limit, as it is fully degenerate with Ψ c .In other words, building upon the observation that our Galaxy rotation curve measurements are performed within an ULDM coherence length in the correlated limit, we redefine a new variable, Ψ 0 c = Ψ c ϕ 2 , which represents the instance of DM in the Milky Way.Finally, as shown by Fig. 22, we report a bump in upper limits at m ϕ ∼ 10 −23 eV, which is at around the maximum-a-posteriori boson mass in Fig. 21.This mass further corresponds to the frequency of the CGW analyzed in EPTA and InPTA Collaborations (2023c).Looking at the posteriors in this mass range, we also find an additional contribution, on top of the CURN process.A similar bump in upper limits is also present in the DR2Full dataset, as discussed in Smarra et al. (2023).However, the Bayes factor of lnB ∼ 0.3 we find in favour of the presence of ULDM signal is still inconclusive.We recommend following up on this bump in future work.

Discussion and outlook
In this paper, we have explored the implications of the common, correlated low-frequency signal observed in the latest data release of the EPTA+InPTA collaboration.Four different datasets were assembled, and the signal was more significant in those including only broadband, high-quality data taken with telescope backends of the new generation.Therefore, we took, as benchmark for our analysis, the signal measured in the DR2new dataset, for which the HD correlation is detected at high significance with a Bayes factor of ≈ 60 or, equivalently, at a p−value of ≈ 0.001, indicative of a ≳ 3σ confidence.The signal can be modelled by a single power law spectrum S h ( f ) as in Eq. ( 3), with best-fit parameters A f =1yr −1 = −13.94+0.23  −0.48 and γ = 2.71 +1.18 −0.71 (PaperIII).We considered several physical processes separately, investigating the implications of the detected signal under the hypothesis that it is generated by that specific process.Our main findings can be summarised as follows.
SMBHBs.The signal is consistent with a cosmic population of merging SMBHBs.Phenomenological models based on galaxy pairs observations can account for the power spectral distribution of the correlated signal.Those models can also be used to predict the chance to detect CGWs in the current data, for which the search has given inconclusive results so far (EPTA and InPTA Collaborations 2023c).According to those models, there is roughly a 50% chance to detect a CGW in DR2new at S/N> 3. The relatively high amplitude of the signal can be used to place constraints on key properties of the cosmic SMBHB population.By exploiting the inference framework developed in Middleton et al. (2016); Chen et al. (2019), we can infer that SMBHBs merge in less than 1Gyr following galaxy mergers, and that the SMBH-stellar bulge relation has a normalization log 10 (M * /M ⊙ ) ≈ 8.4 at a reference stellar mass bulge of 10 11 M ⊙ , in line with recent compiled results (e.g.Kormendy & Ho 2013).Finally, we investigated how the detected signal compares with predictions from state-of-the-art galaxy formation SAMs (Barausse 2012;Izquierdo-Villalba et al. 2022).We found that the high amplitude of the signal imposes nontrivial constraints on galaxy and SMBH evolution models.In particular, boosted SMBH accretion and short SMBHB evolution timescales are needed to match the observed signal, which might lead to difficulties when trying to reproduce independent observables such as the quasar luminosity function.If the SMBHB origin of the observed signal is confirmed, it would be a major breakthrough for observational astrophysics and for our understanding of galaxy formation.This would be the first direct evidence that SMBHBs merge in nature, adding an important observational piece to the puzzle of structure formation and galaxy evolution.
Early Universe.The measured signal would have different implications for each of the physical processes considered in this study.An inflationary origin requires non-standard inflationary scenarios breaking the slow-roll consistency relation, leading to a blue-tilted spectrum.In particular, the measured γ ≈ 2.7 implies n T ≈ 2.3 for a radiation-dominated universe with w = 1/3.A cosmic string origin would allow narrowing down the string tension to values of −11 ≲ log 10 Gµ ≲ −9.5, depending on the specific distribution of loops in the string network.Conversely, the number of kinks cannot be constrained.A GWB induced by (M)HD turbulence at the QCD energy scale can also potentially explain the common red noise, but requires either high turbulent energy densities, of the same order of the radiation energy density, or a characteristic turbulent scale close to the horizon at the QCD epoch.Finally, the measured signal can be produced by the evolution of scalar perturbations at second order only if an excess of their primordial spectrum is present at large wavenumbers, compared to the level derived from CMB observations at small wave numbers.Notably, such an excess would lead to the production of PBHs which can non-negligibly contribute to the CDM density.
ULDM.Finally, we searched for ULDM signatures in DR2new.The search returned a prominent peak in the posterior distribution of the ULDM particle mass around log 10 (m ϕ /eV) ≈ −23.This corresponds to a field oscillation frequency of log 10 ( f /Hz) ≈ −8.3, which is consistent with the CGW candidate examined in EPTA and InPTA Collaborations (2023c).
When a joint ULDM+GWB search is performed, however, the peak in the ULDM mass posterior distribution vanishes, as the data strongly prefer the presence of an HD correlation in the common power.We therefore conclude that an ULDM origin of the detected signal is disfavoured, placing a direct constraint on the abundance of ULDM in our Galaxy.The non-detection in DR2new implies that only about 80% of the DM density in the solar neighbourhood can be attributed to ULDM with −24 < log 10 (m ϕ /eV) < −23.7.More stringent constraints are obtained from DR2full and are presented in a separate paper (Smarra et al. 2023).
It is interesting to remark that the best fit to the measured power-law slope is γ = 2.71, well below the predicted γ = 13/3 expected from a cosmic population of SMBHBs, often indicated as the primary candidate for generating nanohertz GWs.Before leaving the floor to speculations, we caution that this 'inconsistency' might have multiple roots, as mentioned in Sec.3.4.First, γ = 13/3 assumes a population of circular GW-driven binaries.Eccentricity and environmental coupling can easily flatten the low-frequency spectrum.Even for circular binaries, a simple power-law of γ = 13/3 is an ideal limit; small number statistics due to the sparseness of massive, nearby systems results in noisy spectra with a large variance, which can produce different spectral indices when fitted by a power-law.Second, the measured value can be subject to statistical and systematic biases.We have in fact shown that even when a power law with γ = 13/3 is injected into the data, the measured values can be biased due to the statistical realization of the noise.Moreover, the mis-modeling of high-frequency noise can systematically bias the recovered γ, leading to flatter spectra if unaccounted high-frequency noise is present in the data.Therefore, caution should be taken when drawing conclusions from this measurement.
Conversely, a shallow spectral slope might indicate that the GW signal has a different origin, perhaps from some key physical process occurring in the early Universe, as investigated in Section 4. Another possibility, proposed by Lieu et al. (2022) is that GWs can lose energy while propagating through the intergalactic medium via the acceleration of charged particles.Owing to the much higher thermal velocity of electrons over protons and the presence of an intergalactic medium magnetic field, a significant fraction of the GWs can be converted to EM radiation.The ensuing electromagnetic intensity is inversely proportional to frequency.Partial GW loss via this mechanism could explain a flatter-than-expected GWB spectrum.The efficiency of such a mechanism depends on the strength of the intergalactic medium magnetic field.
It is also possible that the observed signal is coming from multiple overlapping GWBs.In our analyses, we have assumed a single origin for all the observed signal power in the EPTA+InPTA data, however, an overlap of GWBs of different origins can cause the spectral shape of the recovered signal to deviate from the expected value of any single GWB.Searching, disentangling and identifying the underlying physical processes will be part of the spectral characterisation moving forward (Moore & Vecchio 2021;Kaiser et al. 2022).
As time goes by and the PTA experiment improves, the lowfrequency GWB will leave an increasingly distinctive signature in the data.Interpreting all of the details of this signature will be necessary to understand the nature of the signal and exploit the full potential of this new window into the Universe.As mentioned in the introduction, a population of SMBHBs is expected to generate a highly non-Gaussian, partially anisotropic and perhaps non-stationary signal.As we increase the timespan of our data, improve our instrumentation and combine more pulsars, the detailed properties of the signal will eventually reveal themselves in the data (e.g.Cornish & Sesana 2013;Taylor & Gair 2013;Taylor et al. 2020;Pol et al. 2022).While many early Universe signals are expected to be isotropic and Gaussian, noticeable exceptions exist, such as GWBs from bursts of cosmic strings.Cross-correlating the power distribution of the nanohertz GW sky to the distribution of massive galaxies and large-scale structures in the low-z universe can eventually provide the key to determining the true nature -astrophysical vs. early Universeof this signal (Rosado & Sesana 2014;Mingarelli et al. 2017).In this respect, combining all of the available high-quality datasets within the IPTA framework is the next step towards the fulfilment of the promises of nanohertz GW science.

Fig. 1 :
Fig. 1: Properties of the common correlated signal detected in DR2new.Left panel: free spectrum of the RMS induced by the excess correlated signal in each frequency resolution bin (with width defined by the inverse of the data span, ∆ f = T −1 ).The straight line is the best power-law fit to the data.Right panel: joint posterior distribution in the A − γ plane.Note that we normalise A to a pivotal frequency f 0 = 10yr −1 .

Fig. 3 :
Fig. 3: Free spectrum violin plot comparing measured (orange) and expected (green) signals.Overlaid to the violins are the 100 Monte Carlo realizations of one specific model; among those, the thick one represents an example of a SMBHB signal consistent with the excess power measured in the data at all frequencies.

Fig. 4 :
Fig. 4: Expected properties of CGWs as a function of frequency.Top panel: free spectrum violin plot comparing the measured signal (orange) to the power distribution of CGWs (green).Empty violins show the full GWB produced by the models for comparison.Bottom panel: the probability of detecting a CGW with S/N> 3 as a function of frequency (green circles, left y−axis scale).The average S/N of CGWs is also shown as red crosses (right y−axis scale).

Fig. 5 :
Fig. 5: A − γ distribution of the measured signal (orange) compared to model predictions (green).1σ and 2σ contours are displayed.Shown are also the marginalised A (left) and γ (right) distributions, with their 1σ credible intervals highlighted as shaded areas.

Fig. 6 :Fig. 7 :
Fig.6: Marginalised posterior distributions for ṅ0 using two SMBHB population models.The orange and green open histograms show marginalised posteriors for the agnostic and astrophysically-informed models, respectively.The filled-green histogram shows the prior for the astrophysically-informed model (the prior for the agnostic model is uniform in the range[−20, 3]).The vertical dotted lines show the 5th and 95th percentiles of the posteriors.

Fig. 8 :
Fig. 8: Posterior distribution of selected parameters for the astrophysically-informed model using nine frequency bins of the free spectrum for the inference.Parameters are described in Sec.3.2.2.
Model (delays increasing left to right) Fig. 10: Binned spectrum of the predicted GWB amplitude for mod-

Fig. 12 :
Fig.12: Predictions for the GWB characteristic strain amplitude at f = 10/yr −1 from a range of L-Galaxies semi-analytical model versions, assuming that h c ( f ) ∝ f −2/3 .The error bars are computed taking into account the cosmic variance.To this end, we divided the Millennium box into 125 sub-boxes and we compute the GWB in each one.The standard deviation provided by the 125 GWBs corresponds to the extension of our error bars.

-
std: standard configuration (Izquierdo-Villalba et al. 2022); -t_DF_x0.1:SMBH dynamical friction (DF) time reduced by a factor of ten; -NO_GAS_HARD: only stellar hardening; -growthDF_x10: accretion boosted by ten in the DF phase; -boostBH1: gas accretion doubled after galaxy mergers and disc instabilities; -boostBH2 gas accretion doubled after galaxy mergers and tripled after disc instabilities.-boostBH1_growthDF_x10: adding accretion boost in the DF phase to model boostBH1; -boostBH2_growthDF_x10: adding accretion boost in the DF phase to model boostBH2; Results are shown in Fig. 12. Changes to the dynamics of SMBHs appear to have a minor effect on the amplitude of the GWB.While shortening the DF time (t_DF_x0.1)allows more SMBHBs to merge within the Hubble time, the most massive ones, which are responsible for the bulk of the GW signal, already have short DF timescales, and the overall GW signal is only mildly increased.Turning off gas hardening results in longer-lived SMBHBs that tend to merge at lower redshifts, resulting in louder GW signals.The effect is, however, negligible.Conversely, the tuning of gas accretion can significantly change the masses of the SMBHBs, resulting in a larger impact on the GWB.Model growthDF_x10 leaves the general population of SMBHs untouched, only boosting the growth of those in the dynamical friction phase.This alone increases the level of the GWB by a factor ≈ 1.5 compared to model std.Finally, the models boostBH1 and boostBH2 increase the gas accretion onto the whole population of SMBHs, making the GWB a factor of 2.5 louder with respect to the baseline model.Boosting accretion onto the whole population and in the hardening phase further amplifies the expected GWB, reaching the upper bound of the measured value (models boostBH1_growthDF_x10 and boostBH2_growthDF_x10 in the figure).Although these models can accommodate the GWB signal measured in PTAs, the boosted accretion and subsequently larger SMBH masses can make it harder to reproduce the observed SMBH mass and quasar luminosity functions (Izquierdo-Villalba et al. 2022).Additionally, more work is required to find a model that can reproduce all observational constraints in the light of the PTA GW signal (Izquierdo-Villalba et al. in preparation).

Fig. 14 :
Fig. 14: Posterior distributions of the recovered GWB from injections on synthetic data mimicking DR2new.Top panel: statistical offset in an ideal dataset.The square marks the injected value and the blue contours are 1σ and 2σ of the recovered posterior.Bottom panel: effect of highfrequency noise mismodeling on the recovery.The orange, blue and green contours are respectively obtained when EFAC= 0.8, 1, 1.2 are used for the recovery (injected EFAC= 1).
g. inTristram et al. (2022);Galloni et al. (2023).Another parameter to consider is the tensor spectral index n T of the tensor perturbations.In the context of slow-roll single-field inflation, these two parameters are linked via the consistency relation r = −8n T .By fixing the consistency relation,Tristram et al. (2022) finds r < 0.032 at 95% CL, while by relaxing it, Planck Collaboration et al. (2020b) finds r < 0.076 and −0.55 < n T < 2.54 at 95% CL.

Fig. 15 :
Fig. 15: 2D posteriors of the tensor-to-scalar ratio (in log 10 ) and the fractional energy density spectral index n T in the PTA frequency range.The 68% and 95% credible regions are displayed.The black dashed line represents the tensor-to-scalar ratio upper bound found in Tristram et al. (2022) assuming single-field slow-roll inflation.

Fig. 16 :
Fig. 16: Comparison of the string tension posteriors for the two string models (BOS and LRS) in case (i), N c = 2 and N k = 0 (Γ = 57).Solid lines assume only a cosmic string background, dashed lines assume both a population of GW-driven circular SMBHBs and cosmic strings.

Fig. 17 :
Fig. 17: SGWB spectra (in terms of log 10 h 2 Ω gw ) for four different early Universe SGWB models considered in this paper.BOS/LRS correspond to a cosmic string background with N c = 2 and N k = 0 (Γ = 57), and log 10 Gµ = −10.1/−10.6.The GWB from turbulence is plotted in solid line for λ * H * = 1, Ω * = 0.3, and T * = 140 MeV.The inflationary spectra is shown for log 10 r = −13.1 and n T = 2.4 (maximum a posteriori value).Power spectrum of the 2nd-order GWB from the scalar curvature perturbations described by the powerlaw model with A 10yr ζ = −2.9 and n s = 2.1 is shown with brown puncture-dot line.The nine first Fourier bins posteriors of the common signal are represented by the grey violin areas.

Fig. 19 :
Fig.19: Results for the monochromatic curvature perturbations described by Eq. 24.Left panel: recovered slopes γ of a simple power-law model as a function of characteristic scale k * of the injected GWB generated by the monochromatic curvature perturbations.The horizontal lines show the theoretical value of γ from a population of circular, GW-driven SMBHBs (grey) and the one obtained in PaperIII (orange).Right panel: 1σ and 2σ contours of the posterior distributions on the amplitude A ζ and characteristic scale of fluctuations k * for DR2new (orange colour).The posterior distribution is overlaid with the current constraints on the primordial power spectrum using Planck data (CMB).The grey colour depicts the 2-σ-confidence intervals.The purple shaded area represents the bounds from spectral distortions(Chluba et al. 2012).For comparison in green we place the prediction of the primordial spectrum of scalar perturbations in the two-field model of inflation described inBraglia et al. (2020) for a range of the model parameters.All three models result in PBH mass functions peaked at ∼ 35 M ⊙ with the brightest line corresponding to the dark matter fraction of PBHs of ∼ 0.01.

Fig. 20 :
Fig. 20: Results for the power-law model of the curvature perturbations described by Eq. (25).Left panel: 1σ and 2σ contours of the posterior

Fig. 21 :Fig. 22 :
Fig. 21: Posterior probabilities for the ULDM amplitude Ψ c and mass m ϕ , from the correlated (top row) and uncorrelated (bottom row) analysis of the DR2new dataset.The pulsar correlated analysis is not shown, but displays the same features.

Figs. 22
and 23 show that the EPTA at current sensitivity is able to constrain the presence of ULDM at the level of the expected DM abundance in the mass range m ϕ ∼ [10 −24 eV, 10 −23.2 eV].The DR2Full analysis performed in Smarra et al. (

Table 1 :
Parameters used for the search and their respective priors.In the correlated limit, φ2 E = φ2 P and can be reabsorbed in a redefinition of Ψ c .
Constraints on the ULDM density ρ ϕ normalised to the DM background value ρ DM = 0.4GeV/cm 3 .The blue, orange and brown lines represent the 95% Bayesian upper limit on ρ ϕ , obtained from the EPTA DR2new dataset with the correlated, uncorrelated and pulsar correlated analysis, respectively.The purple dotted line shows the fiducial local DM density value.