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

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 correlations¹ 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. [87]), the Indian PTA (InPTA Joshi et al. [136]), the North American Nanohertz Observatory for Gravitational Waves (NANOGrav, McLaughlin [175]), the Parkes PTA (PPTA, Manchester et al. [170]), and the Chinese PTA (CPTA, Lee [163]). Although the significance of the signal does not yet reach the $5\sigma$ mark, which is generally accepted as the ‘golden rule’ for a firm detection claim, the evidence reported by the different collaborations ranges between $2\sigma$ and $4\sigma$ (EPTA and InPTA Collaborations 2023b; The Nanograv Collaboration et al. [95]; Reardon et al. [213]; Xu et al. 2023), strongly suggesting a genuine gravitational wave (GW) origin. Awaiting decisive confirmation within the International PTA (IPTA, Verbiest et al. [263]; Perera et al. [194]) framework, with the additional contribution of the MeerKAT PTA (Miles et al. [181]), 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 ground-based interferometers (Abbott et al. [2]). 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 [104]). 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 [206]; Jaffe & Backer [132]). The signal manifests as a stochastic Gaussian process characterised by a power-law Fourier spectrum of delays-advances to pulse arrival times, with characteristic inter-pulsar correlations of general relativity identified by Hellings & Downs [120]. 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. [241]; Ravi et al. [208]). 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. [242]; Babak & Sesana [26]; Kelley et al. [143]) 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 [94]; Kocsis & Sesana [151]; Sesana [236]; Ravi et al. [210]). 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 [182]), 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. [224]; Starobinskii [252]; Fabbri & Pollock [101]; Abbott & Wise [4]). 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 [264]; Damour & Vilenkin 2001; Jones et al. [135]; Dvali & Vilenkin [91]; Damour & Vilenkin [82]), particle production during inflation (Sorbo 2011; Barnaby et al. [28]; Cook & Sorbo 2013; Anber & Sorbo [12]), (magneto-)hydrodynamic turbulence ((M)HD, Kamionkowski et al. [139]; Kosowsky et al. [155]; Dolgov et al. [89]; Caprini & Durrer [60]; Gogoberidze et al. [111]; Caprini et al. [62]), the collision of bubble walls during a first-order primordial phase transition (Kosowsky et al. [157]; Kosowsky & Turner [156]; Caprini et al. [61]; Huber & Konstandin [128]; Jinno & Takimoto [134]; Cutting et al. [77]), sound waves in the aftermath of a first-order phase transition (Hindmarsh et al. [122, 123, 124]), as well as scalar perturbations at second order in cosmological perturbation theory (Baumann et al. [33]; Ananda et al. [11]), 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 [63]).

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 [8]; Maggiore [169]; Caprini & Figueroa [63]). 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, sub-horizon 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 [79] (2001)), or the GWB from particle production during inflation (Cook & Sorbo 2013; Sorbo 2011; Anber & Sorbo [12]). 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. [256]) 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 oscillation of the gravitational potential associated with the presence of ultralight dark matter (ULDM, Smarra et al. [247]⁶, also known as fuzzy dark matter (FDM), in the Galactic halo Khmelnitsky & Rubakov [146]. The existence of ultralight scalars, generally motivated by string-theoretical frameworks [116, 253, 16], 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 [103, 184, 141] or the missing satellite [149, 185] problems, could be disposed of or, at least, mitigated assuming that dark matter is made of ultralight particles. As predicated by Khmelnitsky & Rubakov [146], 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 [186]. 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 $\approx3\sigma$, 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 § 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 (§ 3) or by a number of processes occurring in the early Universe (§ 4). In § 5, we investigate the compatibility of the observed signal with a DM origin and place constraints on ULDM candidates. Finally, in § 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 [31]. 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 [95] and the noise properties of each pulsar are then extracted according to the optimised custom noise models presented in EPTA and InPTA Collaborations [96]. 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 new-generation 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 [120]) 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 free spectra 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 $\alpha=-2/3$ and amplitude $A\approx10^{-15}$, assuming $f_0=1\,\mathrm{yr}^{-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 $\Omega(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. [164]:

where $\gamma=3-2\alpha$. PTAs search for HD correlated time delays with such a power spectrum in the data, and measure the parameters $A$ and $\gamma$. For an observation timespan $T$, measurements are discretised in frequency bins $\Delta 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 Figure 1. The length of the dataset is $T = 10.3$ years, and excess common correlated power is detected in several frequency bins up to $\approx30\,\mathrm{nHz}$ (Figure 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 = 1\,\mathrm{yr}^{-1}$. The data are, however, most informative at the lowest frequencies, while the common power at $1\,\mathrm{yr}^{-1}$ is essentially unconstrained. This naturally leads to a strong degeneracy of the $A-\gamma$ 2D posterior, as shown for example in Figure 1 of PaperIII. Therefore, unless otherwise stated, we change the reference frequency to $f_0 = 10\,\mathrm{yr}^{-1}$, where the data are actually constraining, which results in a weaker dependence of $A$ upon $\gamma$, as shown in the right panel of Figure 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 SMBHBs (at redshifts $z \lesssim1$), 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 $\mathrm{d}^{5}N/(\mathrm{d}z\mathrm{d}m_1\mathrm{d}q\mathrm{d}e\mathrm{d}t_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 [236]

Here, $\mathrm{d}t_r/\mathrm{d}\ln 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. [223], 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 $\mathcal{M}=(m_1m_2)^{3/5}/(m_1+m_2)^{1/5}$, and Eq. (5) takes the familiar form [241]

This can be recast in terms of the comoving number density of merging binaries $\mathrm{d}^{2}n/(\mathrm{d}z\mathrm{d}\mathcal{M})$ [195]

which highlights that, in this case, the expected spectrum follows a power law $h_c \propto f^{-2/3}$, and the only free parameter is its overall amplitude. The latter is set by the function $\mathrm{d}^{2}n/(\mathrm{d}z\mathrm{d}\mathcal{M})$, 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 $d t_{\rm r}/d\ln f_{\rm 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. [210, 238, 142, 70]).

The literature investigating the GWB produced by a population of SMBHBs is vast, dating back to the mid-nineties and early 2000s [206, 132, 268, 240], 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 [241, 242, 208, 159, 142, 45, 246, 131], and empirical models based on observed properties of galaxy pairs coupled to SMBH-host galaxy relations [237, 223, 209, 244], or on the evolution of the SMBH mass function inferred from observations [176]. 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 § 3.1, we perform a semi-quantitative comparison between the measured signal and predictions of an extended version of the Rosado et al. [223] models, hereafter RSG15 including binary eccentricity and environmental coupling. In § 3.2, we exploit the framework developed in Middleton et al. [179]; Chen et al. [68, 69] 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 § 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 [121] and the model developed by Barausse and collaborators [46, 27] (Barausse 2012). We discuss caveats and future directions of improvement in § 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 [237] (S13 hereinafter) and RSG15.

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, $d^3 n/(d m_1 d q d z)$ 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 Figure 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 activity onto the individual SMBHs prior to merger [e.g. 59, 158]. Moreover, hydrodynamical simulations of sub-pc scale binaries, have found most of the accretion to occur on the secondary (i.e. less massive) SMBH [102]. 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. Second, observations of overmassive black holes in the centre of large ellipticals [174] has led to an upward revision of the SMBH-galaxy relations. Contrary to S13 and RSG15, here we consider only those revised relations, namely the ones reported by Kormendy & Ho [153]; McConnell & Ma [173]; Graham & Scott [115]. 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 Figure 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.⁸ This is done by employing the hardening models of Sesana [235] that self-consistently 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 $\mathrm{d}t/\mathrm{d}\ln f_{\mathrm{K,r}}$ in Eq. (5), and thus to reconstruct from the density distribution of merging binaries, $\mathrm{d}^3n/(\mathrm{d}m_1\mathrm{d}q\mathrm{d}z)$, 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, $\mathrm{d}^5N/(\mathrm{d}m_1\mathrm{d}q\mathrm{d}z\mathrm{d}f\mathrm{d}e)$. For each model, we consider 10 values of $e_0 = 0, 0.1, \ldots, 0.9$ and three different normalizations of the stellar density profile, described as $\rho= C \times\rho_0(r/r_0)^{-1.5}$, with $C = 0.1, 1, 10$ (details in Sesana 2010).

In total, we explore $108 \times10 \times3 = 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 $\Delta f = 10.3\ \mathrm{yr}^{-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. [10]; Bonetti & Sesana [44]. 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 Figure 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 GWB 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 $1\ \mathrm{yr}^{-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 thick 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$.⁹ 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 [96], 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 Figure 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 $\approx3.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_{\mathrm{c}} \propto 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-\prod(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. [34], 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-\gamma$ 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 $\log A-\log f$ plane. As described in Sec. 2, we normalise the measurement to $f_0 = 10\,\mathrm{yr}^{-1}$, where our data are informative, which alleviates the $A-\gamma$ degeneracy in the posterior. Results are shown in Figure 5. Although the measured spectrum tends to be shallower than the theoretical one (see also Sec. 3.4), the contours overlap at $2\sigma$ 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. [17]), exploiting the framework laid out in Middleton et al. [179]; ); Chen et al. [68, 69]. The SMBHB population of a given model $\mathbf{M}$ is described by a set of parameters $\{\theta\}$, we then use Bayesian inference to find the posterior distribution $p(\theta|d,\mathbf{M})$ for the parameters $\theta$ of model $\mathbf{M}$ given the observation data $d$:

Here, $p(\theta|\mathbf{M})$ is the prior distribution on the model parameters, $p(d|\theta,\mathbf{M})$ is the likelihood of model $\mathbf{M}$ with parameters $\theta$ of producing the data, and $p(d|\mathbf{M})$ is the evidence. For each set of $\theta$, the likelihood is computed by comparing the signal amplitude produced by $\mathbf{M}$ at frequencies $f_i = i/(10.3\,\mathrm{yr})$, $i = 1,\ldots,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_{\mathbf{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 [186].

Note that the models we use are deterministic in the sense that they have a 1:1 correspondence between $\theta$ and the predicted $h_{\mathrm{c}}(f)$. In reality, given $\theta$, the ensemble of SMBHBs generating the signal is statistically drawn from the deterministic distribution function, which results in a significant scatter of $h_{\mathrm{c}}(f)$, as demonstrated by the individual spectra shown in Figure 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. [179], makes minimal assumptions about the underlying population of SMBHBs. 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 [232] in $z$ and $\mathcal{M}$,

where $t_{\mathrm{R}}$ is the time in the source frame and we assume cosmological parameters from Planck18 [198]. The five model parameters are $\theta=\{\dot{n}_{0},\alpha_{\mathcal{M}},\mathcal{M}_{\star},\beta_{z},z_{0}\}$, where $\dot{n}_{0}$ is the merger rate per unit rest-frame time, co-moving volume, and logarithmic $\mathcal{M}$ interval, and the parameter pairs $\{\alpha_{\mathcal{M}},\mathcal{M}_{\star}\}$ and $\{\beta_{z},z_{0}\}$ control the shape of the $\mathcal{M}$ and $z$ distributions, respectively. The integration limits in $\mathcal{M}$ and $z$ are $10^{6}\leq\mathcal{M}/\mathcal{M}_{\odot}\leq10^{11}$ and $0\leq z\leq5$, 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. [69]. 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: $\{\Phi_{0},\Phi_{I},M_{0},\alpha_{0},\alpha_{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},\alpha_{f},\beta_{f},\gamma_{f}\}$ for the pair fraction, and $\{\tau_{0},\alpha_{\tau},\beta_{\tau},\gamma_{\tau}\}$ 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 $\mathcal{N}\{x,y\}$ is a log normal distribution with mean value $x$ and standard deviation $y$, which adds three further parameters, $\{M_{\star},\alpha_{\star},\epsilon\}$, to the model. The final two parameters describe the eccentricity at SMBHB pairing $\{e_{0}\}$ and the density of the stellar environment $\{\zeta_{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. [69].

Results of the inference

The main results of the inference are shown in Figure 6, Figure 7 and Figure 8. Figure 6 shows the marginalised posterior distribution for the normalization of the merger rate density $\dot{n}_{0}$ from the agnostic model compared to an equivalent parameter derived from the astrophysically-informed model. The constraint on the amplitude of the signal imposes an informative constraint on the normalization of the SMBHB merger density. Using nine frequency bins, the median and central 90% credible regions for $\log_{10}\dot{n}_{0}/[\mathrm{Mpc}^{3}\mathrm{Gyr}]$ are $-1.95^{+2.91}_{-5.87}$ and $-3.51^{+2.90}_{-0.62}$ for the agnostic and astrophysically-informed models, respectively. The measurement essentially constrains the amplitude of the signal, which imposes an informative constraint on $\dot{n}_{0}$. The astrophysically-informed model clearly shows that the signal favours an $\dot{n}_{0}$ at the upper edge of the astrophysical prior. All other parameters of the agnostic model are unconstrained and

the posterior is very similar to the prior (see Appendix A for full posterior distributions for both models).

Figure 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 $\mathcal{M} \approx10^{7}M_{\odot}$ 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 Figure 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 $\mathcal{M} \approx10^{9}M_{\odot}$ 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 Figure 8, the SMBH merger timescale $\tau_0$ following galaxy pairing must be shorter than $\approx1$ Gyr (90% confidence), with the data mildly favouring shorter merger times for massive galaxies at low redshifts (i.e. $\alpha_\tau,\beta_\tau< 0$). Moreover, the data favour a high normalization of the SMBH-bulge mass relation $\log_{10}M_* \approx8.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 § 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, $\beta_f > 1$. Despite the low $\gamma$ 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 $\zeta_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. [239]; Antonini et al. [13]; Klein et al. [148]; Bonetti et al. [46]; Barausse et al. [27]) and L-Galaxies (Henriques et al. [121]; Izquierdo-Villalba et al. [131]). 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 $\alpha= -2/3$.

SAMs and SMBHB delays

In Figure 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. [148] (K+16), Bonetti et al. [46] (B+18) and Barausse et al. [27] (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)”, “Q3-nod (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)¹⁰; (ii) models “popIII-d (K+16)”, “Q3-d (K+16)”, “LS-d (B+18)”, “HS-d (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-noSN-d (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 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 $\gamma= 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 extrapolation – performed as described in Figure 4 of Klein et al. [148] – 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 $\sim4$) 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 Figure 10. Unlike in the case of Figure 9, these predictions were obtained for multiple specific realizations of the SMBHB population, following Sesana et al. (2008)¹¹. 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 Figure 11 we show the theoretical forecasts for $A(f = 1/10\mathrm{yr})$ from a subset of the models presented above (namely those in qualitative agreement with the data in Figure 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 $\gamma= 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/10\mathrm{yr})$ (assuming again $\gamma= 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 $\sim100$ 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 [121, 131], a sophisticated SAM constructed on the dark matter merger trees extracted from the Millennium simulation suite [251, 48]. 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

array of observations including, among others, the galaxy mass function and morphological distribution, the quasar luminosity function and the SMBH-host galaxy relations. Izquierdo-Villalba et al. [131] found that the standard L-Galaxies tuning results in a GWB with $\log_{10} A = -14.9$ at $f_0 = 1\,\mathrm{yr}^{-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:

• std: standard configuration (Izquierdo-Villalba et al. [131]);

• 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 Figure 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 $\approx1.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. [131]). 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).

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 Figure 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 [216]). While gas-driven dynamics is expected to saturate the binary eccentricity at a value $e \approx0.4-0.6$ (Roedig et al. [215]; D’Orazio & Duffell [90]), stellar hardening is known to statistically increase eccentricity without any saturation point (Quinlan [205]), potentially leading to extremely eccentric systems (Sesana [235]). 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 $\tau_0$ derived in Sec. 3.2.

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. [203]; Khan et al. [145, 144]; Nasim et al. [189]; Gualandris et al. [118]). In particular, the binary eccentricity is very sensitive to the initial orbits of the merging galaxies (Gualandris et al. [118]). 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. [196]). Preliminary results show that mergers of massive galaxies occur preferentially on nearly radial orbits, potentially resulting in highly eccentric binaries. Figure 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 inspiral, pairing and hardening of the SMBHs via a hybrid numerical scheme able to model the evolution self-consistently from kpc to mpc scales [86]. Despite the intrinsic stochasticity of the processes driving binary formation and hardening [188], binaries tend to form with a large eccentricity, which then further grows due to encounters with background stars, as also found by scattering experiments [235]. 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 [261]. 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 multi-frequency observations and added a GWB spectrum from an astrophysical population of circular SMBHBs producing a nominal GWB with $A_{\mathrm{1yr}^{-1}} = 2.5 \times10^{-15}$, consistent with the DR2new measurement at $\gamma= 13/3$. We performed 100 experiments by changing the specific noise realization and the sampling of the injected SMBHB population. The analysis was carried out using the ENTERPRISE software package [93].

Two illustrative examples of injection-recovery mismatch are shown in Figure 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 $\gamma= 3.10$. Similar cases have been observed when injecting a pure $\gamma= 13/3$ power law with the createGWB function of libstempo [260]. 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 $\simeq$ 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 Figure 14, we show how the presence of some extra high-frequency 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, $\gamma$ can shift significantly depending on whether the high-frequency 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 common-red 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 [63] 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. [113, 138]) 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 ([224, 252, 101, 4] (Grishchuk 1975)). 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 [63], and references therein). Furthermore, tensor metric perturbations lead to CMB temperature anisotropies and polarisation at large angular scales ([227, 252, 154, 9]). 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.g. in [258, 106]. Another parameter to consider is the tensor spectral index $n_{\mathrm{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_{\mathrm{T}}$. By fixing the consistency relation, Tristram et al. [258] finds $r<0.032$ at 95% CL, while by relaxing it, Planck Collaboration et al. [199] finds $r<0.076$ and $-0.55<n_{\mathrm{T}}<2.54$ at 95% CL.

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 ([161]). 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_{\mathrm{T}}$ (see e.g. [1]). 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$ ([49, 110, 50]). 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. [28, 12, 29] (Sorbo 2011; Cook & Sorbo 2013)) (ii) enhancement of tensor perturbations for example by spectator fields, or space-dependent inflation (see e.g. [30, 37, 36, 105]) (iii) modified gravity theories such as $f(R)$ or Horndeski gravity ([127, 250]) and iv) enhanced scalar perturbations at small scales and/or primordial black holes, which are treated in § 4.4.

Analysis

Similarly to what was done in [161], we constrain the key parameters defining the GWB, $r$ and $n_{\mathrm{T}}$, while being agnostic on the underlying process generating the blue-tilted 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_{\mathrm{T}}$ is expected to vary with the frequency scale considered (see e.g. [109, 108, 22]). However, for simplicity, here we consider a constant $n_{\mathrm{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 [161, 63]

where the second line is valid in the PTA frequency band, and has been obtained by setting $h^{2}\Omega_{\mathrm{rad}}=2.47\cdot10^{-5}$ with $h=0.67$, the amplitude of the scalar spectrum $\mathcal{P}_{\mathcal{R}}^{*}=2\cdot10^{-9}$, and $f_{*}\approx7.7\times10^{-17}$ Hz related to the CMB pivot scale $k_{*}=0.05/\mathrm{Mpc}$ ([197]). $f_{\mathrm{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 Figure 1 (see [186, 160]; for details on the method) to fit the inflationary spectrum of Eq. (14) and obtain posteriors on $\log_{10}r$ and $n_{\mathrm{T}}$. Results are reported in Figure 15. Note that, since $\gamma=5-n_{\mathrm{T}}$, the correlation between the amplitude and spectral index of the signal is compatible with Figure 5. The 90% credible (symmetric) intervals are $\log_{10}r=-12.18^{+8.81}_{-7.00}$ and $n_{\mathrm{T}}=2.29^{+0.87}_{-1.11}$. The obtained value of $n_{\mathrm{T}}$ corresponds to a PSD spectral index of $\gamma\approx2.7$, as in Figure 1. The excessively small value of $r$ is a consequence of the simplistic parameterisation of Eq. (14), which assumes a constant $n_{\mathrm{T}}$ at all scales. The fractional energy density spectrum obtained from the maximum a posteriori parameter values is plotted in Figure 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. [138] 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 $\gamma=4$ as opposed to $\gamma=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 $\gamma< 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 $\gamma= 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. [161], $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. [161] (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 - \gamma\simeq2.3$, which is directly linked to the best-fit PSD spectral index $\gamma\simeq2.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 \ne1/3$, the relation between the PSD spectral index $\gamma$ and the primordial tensor spectral index $n_T$ changes to ([18, 63])

again with $\gamma\simeq2.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 \gtrsim0.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 \simeq0$) for the best fit value $\gamma\simeq2.7$. However, by broadening the range of possible values to $\gamma\ge3$, $n_T \simeq0$ 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 ([147, 125]); they are generic predictions of most Grand Unification Theories scenarios ([133]). These 1D objects are characterised by the string tension $G\mu$ (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 ([3, 6]). 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\mu\sim10^{-11}$ ([20]). 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 ([23]).

The incoherent sum of these GW bursts also produces a stochastic GWB (SGWB), which has been looked for by LVK ([3, 6]). LISA is expected to be able to detect a SGWB from cosmic strings with tensions as low as $G\mu\sim10^{-17}$ ([24]). The strings background has already been looked for in EPTA ([230, 162]), in NANOGrav ([43, 92]) and PPTA ([39, 71]). 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\mu\sim10^{-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(\ell,t)$, where $\ell$ is the invariant length of the loop (defined by its total energy divided by $\mu$) and $t$ is cosmic time. We focus on the most up-to-date loop distributions, calibrated on large scales with Nambu-Goto simulations ([167, 41]). Loops are formed from the intercommutation of super-Horizon infinite strings, with a maximal size $\ell\simeq0.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. [42]) and LRS (Lorenz et al. [167]), 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 [25]. The explicit expressions for the two loop distributions are given in Abbott et al. [3], and both have been considered by the LVK [6] and LISA [24] collaborations. The third observing run of LVK placed constraints on $G\mu$, based on the non-detection of a SGWB. These are $G\mu\lesssim10^{-8}$ for BOS and $G\mu\lesssim10^{-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 [214, 25, 162].

Other than $G\mu$, a further parameter appears in $n(\ell,t)$, namely $\Gamma$, which describes the rate at which loops lose energy through gravitational radiation: $\dot{\ell}=-\Gamma G\mu$. If $\Gamma$ 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 [25]. The value of $\Gamma$ 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_{\mathrm{c,k}}$ is the number of cusps/kink events per oscillation period of the loop, $N_{\mathrm{k}}^2/4$ the number of kink-kink collisions, and

where $g_2=\sqrt{3}/4$, $g_1^{\mathrm{c}}\approx0.85$, $g_1^{\mathrm{k}}\approx0.29$, $g_1^{\mathrm{kk}}\approx0.1$.

In this paper, we consider 2 cases. For the first model (i) $N_{\mathrm{c}}=2$ and $N_{\mathrm{k}}=0$, for which $\Gamma=57$ [162], a value close to that observed in numerical simulations of individual loops [259, 40]. Therefore, the only free parameter for this model is $G\mu$. As for the second model (ii), we fix $N_{\mathrm{c}}=1$ and allow $N_{\mathrm{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\mu,N_{\mathrm{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 $\Lambda$CDM cosmology with the Planck-2018 fiducial parameters [7], 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. [162], which was based on six pulsars with the best timing precision from the EPTA early Data Release 2 [67]. We now analyse the DR2new dataset with 25 pulsars, using the (computationally heavy) hierarchical likelihood data analysis method described in Paper III. We sample the SGWB parameters ($G\mu$) or $(G\mu,N_{\mathrm{k}})$ as well as the individual pulsar noise parameters.

Solid lines in Figure 16 show the posterior distribution on $\log_{10}(G\mu)$ in case (i), for which $N_{\mathrm{k}}=0$ and $\Gamma=57$. The string tension 90% credible (symmetric) interval is $\log_{10}G\mu=-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 Figure 17 where, for each model, we set $G\mu$ 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\mu$) are highly correlated, since both provide a possible explanation for the detected signal. As a result, the posterior on $\log_{10}G\mu$ no longer has compact support, but a tail to lower values (see the dashed lines in Figure 17). We therefore extract the 95-quantile of the string tension posterior to obtain an upper bound of $\log_{10}G\mu<-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 [162], 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 Figure 17.

For case (ii), we obtain very similar results to those discussed in Leclere et al. [162]. Namely, we obtain quasi non-informative posteriors for $N_{\mathrm{k}}$, showing that the data can be equally explained by a population of kinky loops with $N_{\mathrm{k}}\gtrsim120$. In other words, we cannot extract any upper bound on the number of kinks, since this quantity is degenerate with $G\mu$.

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 [266, 139], or can be driven by pre-existing primordial magnetic fields [204, 53]. If the (magneto-)hydrodynamic turbulence were present around the QCD epoch, when the Universe had a temperature of $T_* \sim100$ 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 $\lambda_* \simeq O(\mathcal{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 [234, 267, 265, 58] (Middeldorf-Wygas et al. 2020; Middeldorf-Wygas et al. 2022). 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 [157, 156, 61, 128, 134, 77, 122, 123, 124]. This was analysed for PTAs, e.g. in [186, 270] (); Arzoumanian et al. 2021). In what follows, we focus on the GWB generated by decaying (M)HD turbulence [139, 155, 89, 60, 111, 62]. This has been recently studied via numerical simulations in Roper Pol et al. (2020a,b); Brandenburg et al. [52]; Roper Pol et al. [219]. In particular, Roper Pol et al. [219] 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. [21]. This assumption leads to the following GWB signal:

where $\Omega_*$ is the ratio of the (M)HD turbulent energy density to the radiation one, and $\lambda_* \mathcal{H}_*$ is the ratio of the characteristic length scale of the turbulence, $\lambda_*$, to the comoving Hubble horizon $\mathcal{H}_*^{-1}$ at the QCD epoch. The parameter $\mathcal{A} \sim1.75 \times10^{-3}$ is the efficiency of GW production, estimated in Roper Pol et al. [219]. The function $F_{\rm 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_{\rm turb}(f)$, is

Description of the model

The presence of bulk velocity and magnetic fields produce anisotropic stresses, which in turn act as a source of GWs where $\mathcal{B} \simeq50 (\lambda_{*} \mathcal{H}_{*})^{-2}$ is a normalising factor, and $\delta t_{\mathrm{fin}}$ denotes the effective duration of the turbulence. The latter can be estimated, from the numerical simulations performed in Roper Pol et al. [219], to be $\delta t_{\mathrm{fin}} \simeq2\lambda_{*}/\sqrt{1.5}\,\Omega_{*}$. The function $p_{\Pi}(\lambda_{*}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 $\lambda_{*}$, it is constant for $f < 1/\lambda_{*}$. Furthermore, it decays as $f^{-11/3}$ for $f \gtrsim1/\lambda_{*}$, 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/\delta t_{\mathrm{fin}}$, $f$ at intermediate frequencies $1/\delta t_{\mathrm{fin}} < f < 1/\lambda_{*}$, and $f^{-8/3}$ at large frequencies $f > 1/\lambda_{*}$.

The GWB produced from vortical (M)HD turbulence is therefore determined by three parameters: the temperature scale $T_{*}$, the turbulence strength $\Omega_{*}$, and the turbulence characteristic length scale $\lambda_{*}\mathcal{H}_{*}$. By causality, $\lambda_{*}\mathcal{H}_{*}$ is bound to be smaller than one. In general, also $\Omega_{*} \lesssim1$, 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 $\Omega_{*} \lesssim 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 Figure 1. We use $\log_{10}$-uniform priors for the model parameters, choosing $\log_{10}(\lambda_{*}\mathcal{H}_{*}) \in[-3,0]$, $\log_{10}\Omega_{*} \in[-2,0]$, and $\log_{10}(T_{*}/1\,\mathrm{MeV}) \in[1,3]$. The 2D posteriors obtained are shown in Figure 18.

For values of $\Omega_{*}$ 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 $\lambda_{*}\mathcal{H}_{*}$ is close to 1 and the peak frequency lies within the PTA frequency range, implying $T_{*} \sim60$ MeV. However, at frequencies around the peak, the signal corresponds to a power spectral density for the residuals steeper than $\gamma\sim4$, which cannot fit the data well. For this reason, values of $\Omega_{*} \lesssim0.1$ are disfavoured.

For larger values of $\Omega_{*}$, the $f^{3}$ part of the spectrum at frequencies below $\delta t_{\mathrm{fin}}^{-1} \propto\sqrt{\Omega_{*}/(\lambda_{*}\mathcal{H}_{*})}\times\mathcal{H}_{*}(T_{*})$ can enter the PTA band with a sufficiently high amplitude. Furthermore, the distance between the break at $\delta t_{\mathrm{fin}}^{-1}$ and the spectral peak at $1/\lambda_{*}$ becomes minimal in the limit $\Omega_{*} \sim1$. 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 $\lambda_{*}\mathcal{H}_{*}$ and $\Omega_{*}$ from the signal amplitude (see Equation (21)), and the degeneracy between $\lambda_{*}\mathcal{H}_{*}$ and $T_{*}$ from the break at $1/\delta t_{\mathrm{fin}}$ (note that the dependence of the latter on $\sqrt{\Omega_{*}}$ is subdominant).

The model therefore provides a good fit to the data in the limit of large $\Omega_{*}$, 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 [166]; Baumann [32]). However, scalar curvature perturbations will source propagating tensorial modes (GWs) at the 2nd order in perturbation theory (Tomita [257]; Matarrese et al. [172, 171]; Noh & Hwang [191]; Carbone & Matarrese [64]; Ananda et al. [11]; Baumann et al. [33]). 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. [199]) suggests that the power spectrum of the curvature perturbations is nearly scale-invariant with the amplitude $A_{\zeta} \sim2 \times10^{-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 $\Omega_{\mathrm{gw}} \sim10^{-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. [56]; Braglia et al. [51]; 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 \sim10^{-3} - 10^{-1}\,\mathrm{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 \sim10^{6} - 10^{8}\,\mathrm{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 [229]; Bugaev & Klimai [55]; Chen et al. [72]; Dandoy et al. [83]; Zhao & Wang [273]).

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_{\zeta}$ is a dimensionless amplitude and $k^{*}$ is a wavenumber at which the monochromatic power spectrum has a Dirac-delta peak. For the second case:

where $n_s$ characterises the slope and $k^{10\mathrm{yr}}$ is a normalizing scale $k_{10\mathrm{yr}}=2\pi/(10\mathrm{yr}\times c)$, so that $A_{\zeta}^{10\mathrm{yr}}$ 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 ([152, 100]):

where $\tilde{k}=k/k^{*}$ and $\Theta$ is the Heaviside theta function. In spite of being nonphysical, the $\delta$-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. [56]).

In the second case of a more general (and more realistic) power-law spectrum, the result can only be obtained numerically (Kohri & Terada [152]):

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 [152].

Parameter	Description	Prior	Occurrence
	White Noise ($\sigma= E_f^2\sigma_{ToA}^2 + E_q^2$)
$E_f$	EFAC per backend/receiver system	Uniform [0, 10]	1 per pulsar
$E_q$	EQUAD per backend/receiver system	Log$_{10}$-Uniform [$-10$, $-5$]	1 per pulsar
	Red Noise
$A_{\mathrm{red}}$	Red noise power-law amplitude	Log$_{10}$-Uniform [$-20$, $-6$]	1 per pulsar
$\gamma_{\mathrm{red}}$	Red noise power-law spectral index	Uniform [0, 10]	1 per pulsar
	ULDM
$\Psi_c$	ULDM amplitude	Log$_{10}$-Uniform [$-20$, $-12$]	1 for PTA
$m_\phi$ [eV]	ULDM mass	Log$_{10}$-Uniform [$-24$, $-22$]	1 for PTA
$\hat{\phi}_E^2$	Earth factor	$e^{-x}$	1 for PTA
$\hat{\phi}_P^2$	Pulsar factor	$e^{-x}$	1 per pulsar
$\gamma_E$	Earth signal phase	Uniform [0, $2\pi$]	1 for PTA
$\gamma_P$	Pulsar signal phase	Uniform [0, $2\pi$]	1 per pulsar
	Common spatially Uncorrelated Red Noise (CURN)
$A_{\mathrm{GWB}}$	Common process strain amplitude	Log$_{10}$-Uniform [$-20$, $-6$]	1 for PTA

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 [228]). 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¹³, $T_{\mathrm{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_{\zeta}$ and $\log_{10} A_{\zeta}^{10\mathrm{yr}}$, uniform in $[4,12]$ for $\log_{10}(k^{*}/\mathrm{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 Figure 19 and Figure 20, respectively. The 2D posterior distribution of the model parameters of the monochromatic model is depicted with orange contours on the right panel of Figure 19. One may notice that the region of the highest probability are strongly elongated due to a strong positive correlation between $A_{\zeta}$ 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}/\mathrm{Mpc}^{-1})=7.6$ and $\log_{10} A_{\zeta}^{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 Figure 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)\sim f^{-\gamma}$ used in PaperIII to model an arbitrary common red noise process. After a rapid decrease, the recovered slope stabilises at $\gamma\sim3.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 Figure 20 with the following means and 1-$\sigma$ uncertainties: $\log_{10} A_{\zeta}=-2.94^{+0.42}_{-0.46}$ and $n_s=2.11^{+0.25}_{-0.32}$.

On the right panels of Figure 19 and Figure 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 Figure 20 are 1$\sigma$ and 2$\sigma$ 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. [130]; ); Germani & Prokopec [107]; ); Di & Gong (2018); ); Cai et al. [57]; ); Bian et al. (2018); ); Byrnes et al. [56]; ); Motohashi et al. (2020); ); Braglia et al. [51]; ); 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. [65], and reference therein). The PBH formation from cosmological perturbations has been extensively explored Sasaki et al. [231]. On the radiation-dominated stage, the PBH mass is related to the mass inside the horizon at the time of the perturbation entering, $M\sim M_H\approx m_{\mathrm{Pl}}^2 t=8M_{\odot}(100\mathrm{MeV}/T)^2(40/g_*)^{1/2}$, where $m_{\mathrm{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. [51], a peak around $k\sim2\times10^{6}\,\mathrm{Mpc}^{-1}$ could explain the power $A_{\zeta}\sim10^{-2}$ and simultaneously lead to the production of PBHs with masses peaked at $\sim35M_{\odot}$ (see also the analysis of the NANOGrav results in Vaskonen & Veermäe [262]). Interestingly, such a peak seems to be observed in the chirp mass distribution of LVK merging binary black holes (Abbott et al. [5]). 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 [146], 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 [98]. Additionally, the ULDM search with DR2full is performed in Smarra et al. [247].

Dark matter currently constitutes approximately 26% of the energy density of the Universe, as confirmed by, e.g., galactic rotation curves [226, 225, 84], baryonic acoustic oscillations and cosmic microwave background measurements [35, 7] as well as galaxy surveys [99]. The standard cold dark matter (CDM) picture, whose leading candidates are the Weakly Interacting Massive Particles (WIMPs) [14] and the QCD axion [168], 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(\mathrm{kpc})$. Among these, the cusp-core problem [103, 184, 141] 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 [149, 184].

While a thorough understanding of baryonic physics feedback mechanisms [190, 114, 54, 66, 193, 212, 187] might help to alleviate some of these issues, they can be more easily disposed of assuming that DM is an ultralight ($m_{\phi} \sim10^{-22}\,\mathrm{eV}$) scalar/pseudoscalar or axion-like field, whose astrophysically large ($O(\mathrm{kpc})$) de Broglie wavelength suppresses power on small scales. Moreover, ultralight scalars are also generally present in string theory compactifications [116, 253, 16], which makes them interesting candidates for new physics as well. CMB-based arguments are used to constrain $m_{\phi} \gtrsim10^{-24}\,\mathrm{eV}$ [126], while Lyman-$\alpha$ bounds push the limit up to $m_{\phi} \gtrsim10^{-21}\,\mathrm{eV}$, provided that the ultralight particles account for more than 30% of the full dark matter budget [15, 150, 192, 217] (Iršić et al. 2017). A lower limit of $m_{\phi} \sim10^{-19}\,\mathrm{eV}$ is claimed by studies of ultra-faint dwarf (UFD) galaxies [119, 78], but a wide consensus is yet to be reached. In a seminal work, Khmelnitsky & Rubakov [146] 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 [233, 272], as the aforementioned Lyman-$\alpha$ 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 pulsars as described in Khmelnitsky & Rubakov [146].

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_{\phi}$ [146]. The periodic displacement that ULDM induces on the TOAs of signal from a pulsar $P$ can be written as follows [202]:

where

where $\gamma_{P}$ ($\gamma_{E}$) is a pulsar (Earth) dependent phase and $\hat{\phi}_{P}^{2}$ ($\hat{\phi}_{E}^{2}$) 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_{\phi} \sim10^{-3}$ for the ULDM velocity, the region in which the scalar field oscillates coherently, i.e. with the same value of $\hat{\phi}$, is spanned by the coherence length:

In particular, notice that $\hat{\phi}_{E}^{2}$ and $\hat{\phi}_{P}^{2}$ 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. [247], 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 $\sim\mathrm{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 $\sim20$ kpc), but larger than the average inter-pulsar and pulsar-Earth distance. More specifically, the correlated regime holds for masses lower than $m_{\phi} \sim2 \times10^{-24}\,\mathrm{eV}$; the pulsar correlated regime for $2 \times10^{-24}\,\mathrm{eV} \lesssim m_{\phi} \lesssim5 \times10^{-23}\,\mathrm{eV}$ and the uncorrelated limit for $m_{\phi} \gtrsim5 \times10^{-23}\,\mathrm{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. [247]. As an example, marginalised posterior distributions for $\Psi_{\mathrm{c}}$ and $m_{\phi}$ are plotted in Figure 21. All the three limits peak at a ULDM mass $m_{\phi} \sim10^{-23}\,\mathrm{eV}$ with amplitude $\Psi_{\mathrm{c}} \sim10^{-15}$, which translates into a density $\rho_{\phi} \sim0.6\,\mathrm{GeV}/\mathrm{cm}^{3}$ of the scalar field. While this value is higher than the fiducial local DM density $\rho_{\mathrm{DM}} \approx0.4\,\mathrm{GeV}/\mathrm{cm}^{3}$, it is well within the measurement uncertainties [47, 211, 245, 85].

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

Table 1: Parameters used for the search and their respective priors. In the correlated limit, $\hat{\phi}_E^2 = \hat{\phi}_P^2$ and can be reabsorbed in a redefinition of $\Psi_c$.

agreement with the fact that the data support HD correlation, which naturally favours GWs over ULDM in a joint search. Thus

we put 95% upper limit on the amplitude $\Psi_c$ and the density $\rho_\phi$ of the scalar field.

Figure 22 and Figure 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_\phi\sim[10^{-24}\,\mathrm{eV}, 10^{-23.2}\,\mathrm{eV}]$. The DR2Full analysis performed in Smarra et al. [247] 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 \approx3 \cdot10^{-9}\,\mathrm{Hz}$, with $T \approx10\,\mathrm{yr}$. In fact, while it might naively be thought that the ULDM field needs to complete an oscillation in the timescale $\tau_{\rm PTA}$ of the PTA experiment to produce a detectable effect, we point out that an ULDM wave with frequency $m_\phi< 1/\tau_{\rm PTA}$ can still be approximated by an expansion in powers of $m_\phi t$ [140], though the sensitivity of PTAs will be reduced due to degeneracies with the timing model [207]. Relying on the robust CMB bounds mentioned before, we fix the lower end of our search at $m_\phi= 10^{-24}\,\mathrm{eV}$. Importantly, we remove the $\hat{\phi}_{E} = \hat{\phi}_{P}$ parameter from the search in the correlated limit, as it is fully degenerate with $\Psi_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, $\Psi_c^{0} = \Psi_c\phi^2$, which represents the instance of DM in the Milky Way. Finally, as shown by Figure 22, we report a bump in upper limits at $m_\phi\sim10^{-23}\,\mathrm{eV}$, which is at around the maximum-a-posteriori boson mass in Figure 21. This mass further corresponds to the frequency of the CGW analyzed in EPTA and InPTA Collaborations [98]. 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. [247]. However, the Bayes factor of $\ln B \sim0.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 $\approx60$ or, equivalently, at a $p$-value of $\approx0.001$, indicative of a $\gtrsim3\sigma$ confidence. The signal can be modelled by a single power law spectrum $S_h(f)$ as in Eq. (3), with best-fit parameters $A_{\rm yr^{-1}}=-13.94^{+0.24}_{-0.48}$ and $\gamma=2.71^{+1.18}_{-0.71}$ (Paper III).

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.

SMBHs. The signal is consistent with a cosmic population of merging SMBHBs. Phenomenological models based on galaxy 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 $\geq3$. 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. [179]; Chen et al. [69], we can infer that SMBHBs merge in less than 1 Gyr following galaxy mergers, and that the SMBH–stellar bulge relation has a normalization $\log_{10}(M_*/M_\odot) \approx8.4$ at a reference stellar mass bulge of $10^{11}M_\odot$, in line with recent compiled results [e.g. 153]. 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 non-trivial 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 $\gamma\approx2.7$ implies $n_T \approx2.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 \lesssim\log_{10}G\mu\lesssim-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_\phi/\mathrm{eV}) \approx-23$. This corresponds to a field oscillation frequency of $\log_{10}(f/\mathrm{Hz}) \approx-8.3$, which is consistent with the CGW candidate examined in EPTA and InPTA Collaborations [98]. 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_\phi/\mathrm{eV}) < -23.7$. More stringent constraints are obtained from DR2full and are presented in a separate paper ([247]).

It is interesting to remark that the best fit to the measured power-law slope is $\gamma= 2.71$, well below the predicted $\gamma= 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 § 3.4. First, $\gamma= 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 $\gamma= 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 $\gamma= 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 $\gamma$, 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 § 4. Another possibility, proposed by Lieu et al. [165] 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 GWB 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 ([186, 138]).

As time goes by and the PTA experiment improves, the low-frequency GWB will leave an increasingly distinctive signature in the data. Interpreting all 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 GWB 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 Universe – of this signal [222, 183]. 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.

Acknowledgements

The European Pulsar Timing Array (EPTA) is a collaboration between European and partner institutes, namely ASTRON (NL), INAF/Osservatorio di Cagliari (IT), Max-Planck-Institut für Radioastronomie (GER), Nançay/Paris Observatory (FRA), the University of Manchester (UK), the University of Birmingham (UK), the University of East Anglia (UK), the University of Bielefeld (GER), the University of Pisa (FRA), the University of Milan-Bicocca (IT), the Foundation for Research and Technology, Hellas (GR), and Peking University (CHN), with the aim to provide high-precision pulsar timing to work towards the direct detection of low-frequency gravitational waves. An Advanced Grant of the European Research Council allowed to implement the Large European Array for Pulsars (LEAP) under Grant Agreement Number 227947 (PI M. Kramer). The EPTA is part of the International Pulsar Timing Array (IPTA); we thank our IPTA colleagues for their support and help with this paper and the external Detection Committee members for their work on the Detection Checklist. Part of this work is based on observations with the 100-m telescope of the Max-Planck-Institut für Radioastronomie (MPIfR) at Effelsberg in Germany. Pulsar research at the Jodrell Bank Centre for Astrophysics and the observations using the Lovell Telescope are supported by a Consolidated Grant (ST/T000414/1) from the UK’s Science and Technology Facilities Council (STFC). ICN is also supported by the STFC doctoral training grant ST/T506291/1. The Nançay radio Observatory is operated by the Paris Observatory, associated with the French Centre National de la Recherche Scientifique (CNRS), and partially supported by the Region Centre in France. We acknowledge financial support from “Programme National de Cosmologie et Galaxies” (PNCG), and “Programme National Hautes Energies” (PNHE) funded by CNRS/INSU-IN2P3-INP, CEA and CNES, France. We acknowledge financial support from Agence Nationale de la Recherche (ANR-18-CE31-0015), France. The Westerbork Synthesis Radio Telescope is operated by the Netherlands Institute for Radio Astronomy (ASTRON) with support from the Netherlands Foundation for Scientific Research (NWO). The Sardinia Radio Telescope (SRT) is funded by the Department of University and Research (MIUR), the Italian Space Agency (ASI), and the Autonomous Region of Sardinia (RAS) and is operated as a National Facility by the National Institute for Astrophysics (INAF). The work is supported by the National SKA programme of China (2020SKA0120100), Max-Planck Partner Group, NSFC 11690024, CAS Cultivation Project for FAST Scientific. This work is also supported as part of the “LEGACY” MPG-CAS collaboration on low-frequency gravitational wave astronomy. JA acknowledges support from the European Commission (Grant Agreement number: 101094354). JA and SCha were partially supported by the Stavros Niarchos Foundation (SNF) and the Hellenic Foundation for Research and Innovation (H.F.R.I.) under the 2nd Call of the “Science and Society – Action Always Strive for excellence – Theodoros Papazoglou” (Project Number: 01431). AC acknowledges support from the Paris Île-de-France Region. AC, AF, ASe, ASa, EB, DI, GMS, MBo acknowledge financial support provided under the European Union’s Horizon 2020 ERC Consolidator Grant “Binary Massive Black Hole Astrophysics” (B Massive, Grant Agreement: 818691). GD, KLi, RK and MK acknowledge support from European Research Council (ERC) Synergy Grant “BlackHoleCam”, Grant Agreement Number 610058. This work is supported by the ERC Advanced Grant “LEAP”, Grant Agreement Number 227947 (PI M. Kramer). AV and PRB are supported by the UK’s Science and Technology Facilities Council (STFC; grant ST/W000946/1). AV also acknowledges the support of the Royal Society and Wolfson Foundation. JPW acknowledges support by the Deutsche Forschungsgemeinschaft (DFG) through the Heisenberg programme (Project No. 433075039) and by the NSF through AccelNet award #2114721. NKP is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – Projekt-nummer PO 2758/1–1, through the Walter–Benjamin programme. ASa thanks the Alexander von Humboldt foundation for a Humboldt fellowship for postdoctoral researchers. AP, DP and MBu acknowledge support from the research grant “iPeska” (P.I. Andrea Possenti) funded under the INAF national call PRIN-SKA/CIA approved with the Presidential Decree 70/2016 (Italy). RNC acknowledges financial support from the Special Account for Research Funds of the Hellenic Open University (ELKE-HOU) under the research programme “GRAVPUL” (grant agreement 319/10-10-2022). EvdW, CGB and GHJ acknowledge support from the Dutch National Science Agenda, NWA Startimpuls – 400.17.608. BG is supported by the Italian Ministry of Education, University and Research within the PRIN 2017 Research Program Framework, n. 2017SYRTCN. LS acknowledges the use of the HPC system Cobra at the Max Planck Computing and Data Facility. The Indian Pulsar Timing Array (InPTA) is an Indo-Japanese collaboration that routinely employs TIFR’s upgraded Giant Metrewave Radio Telescope for monitoring a set of IPTA pulsars. BCJ, YG, YM, SD, AG and PR acknowledge the support of the Department of Atomic Energy, Government of India, under Project Identification # RTI 4002. BCJ, YG and YM acknowledge support of the Department of Atomic Energy, Government of India, under project No. 12-R&D-TFR-5.02-0700 while SD, AG and PR acknowledge support of the Department of Atomic Energy, Government of India, under project no. 12-R&D-TFR-5.02-0200. KT is partially supported by JSPS KAKENHI Grant Numbers 20H00180, 21H01130, and 21H04467, Bilateral Joint Research Projects of JSPS, and the ISM Cooperative Research Program (2021-ISMCRP-2017). AS is supported by the NANOGrav NSF Physics Frontiers Center (awards #1430284 and 2020265). AKP is supported by CSIR fellowship Grant number 09/0079(15784)/2022-EMR-I. SH is supported by JSPS KAKENHI Grant Number 20J20509. KN is supported by the Birla Institute of Technology & Science Institute fellowship. AmS is supported by CSIR fellowship Grant number 09/1001(12656)/2021-EMR-I and T-641 (DST-ICPS). TK is partially supported by the JSPS Overseas Challenge Program for Young Researchers. We acknowledge the National Supercomputing Mission (NSM) for providing computing resources of ‘PARAM Ganga’ at the Indian Institute of Technology Roorkee as well as ‘PARAM Seva’ at IIT Hyderabad, which is implemented by C-DAC and supported by the Ministry of Electronics and Information Technology (MeitY) and Department of Science and Technology (DST), Government of India. DD acknowledges the support from the Department of Atomic Energy, Government of India through Apex Project - Advance Research and Education in Mathematical Sciences at IMSc. The work presented here is a culmination of many years of data analysis as well as software and instrument development. In particular, we thank Drs. N. D’Amico, P. C. C. Freire, R. van Haasteren, C. Jordan, K. Lazaridis, P. Lazarus, R. L. Lentati, O. Lö hmer and R. Smits for their past contributions. We also thank Dr. N. Wex for supporting the calculations of the galactic acceleration as well as the related discussions. We would like to thank Prof. Drs. Alexey Starobinsky, Sergei Blinnikov and Alexander Dolgov for discussions on the early Universe physics. MHM acknowledges the support of the UK Space Agency, Grant No. ST/V002813/1 and ST/X002071/1. Some of the computations described in this paper were performed using the University of Birmingham’s BlueBEAR HPC service, which provides a High Performance Computing service to the University’s research community. See http://www.birmingham.ac.uk/bear for more details. The EPTA is also grateful to staff at its observatories and telescopes who have made the continued observations possible.

Author contributions. The EPTA is a multi-decade effort and all authors have contributed through conceptualisation, funding acquisition, data-curation, methodology, software and hardware developments as well as (aspects of) the continued running of the observational campaigns, which includes writing and proofreading observing proposals, evaluating observations and observing systems, mentoring students, developing science cases. All authors also helped in (aspects of) verification of the data, analysis and results as well as in finalising the paper draft. Specific contributions from individual EPTA members are listed in the CRediT¹⁴ format below. InPTA members contributed in uGMRT observations and data reduction to create the InPTA data set which is employed while assembling the DR2full+ and DR2new+ data sets.

{sup}43` Department of Astronomy,School of Physics, Peking University, Beijing 100871, P. R. China

{sup}44` National Astronomical Observatories, Chinese Academy of Sciences, Beijing 100101, P. R. China

{sup}45` E.A. Milne Centre for Astrophysics, University of Hull, Cottingham Road, Kingston-upon-Hull, HU6 7RX, UK

{sup}46` Centre of Excellence for Data Science, Artificial Intelligence and Modelling (DAIM), University of Hull, Cottingham Road, Kingston-upon-Hull, HU6 7RX, UK

{sup}47` Laboratory of Astrophysics, École Polytechnique Fédérale de Lausanne, CH-1015 Lausanne, Switzerland

{sup}48` Department of Physics, BITS Pilani Hyderabad Campus, Hyderabad 500078, Telangana, India

{sup}49` Joint Astronomy Programme, Indian Institute of Science, Bengaluru, Karnataka, 560012, India

{sup}50` Arecibo Observatory, HC3 Box 53995, Arecibo, PR 00612, USA

{sup}51` IRFU, CEA, Université Paris-Saclay, F-91191 Gif-sur-Yvette, France

{sup}52` Raman Research Institute India, Bengaluru, Karnataka, 560080, India

{sup}53` Institut für Physik und Astronomie, Universität Potsdam, Haus 28, Karl-Liebknecht-Str. 24/25, 14476, Potsdam, Germany

{sup}54` Kazan Federal University, 18 Kremlyovskaya, 420008 Kazan, Russia

{sup}55` Department of Physics, IISER Bhopal, Bhopal Bypass Road, Bhauri, Bhopal 462066, Madhya Pradesh, India

{sup}56` Ollscoil na Gaillimhe — University of Galway, University Road, Galway, H91 TK33, Ireland

{sup}57` Center for Gravitation, Cosmology, and Astrophysics, University of Wisconsin-Milwaukee, Milwaukee, WI 53211, USA

{sup}58` Division of Natural Science, Faculty of Advanced Science and Technology, Kumamoto University, 2-39-1 Kurokami, Kumamoto 860-8555, Japan

{sup}59` International Research Organization for Advanced Science and Technology, Kumamoto University, 2-39-1 Kurokami, Kumamoto 860-8555, Japan

{sup}60` Laboratoire Univers et Théories LUTh, Observatoire de Paris, Université PSL, CNRS, Université de Paris, 92190 Meudon, France

{sup}61` Max Planck Institute for Gravitational Physics (Albert Einstein Institute), Leibniz Universität Hannover, Callinstrasse 38, D-30167, Hannover, Germany

{sup}62` Florida Space Institute, University of Central Florida, 12354 Research Parkway, Partnership 1 Building, Suite 214, Orlando, 32826-0650, FL, USA

{sup}63` Ruhr University Bochum, Faculty of Physics and Astronomy, Astronomical Institute (AIRUB), 44780 Bochum, Germany

{sup}64` Advanced Institute of Natural Sciences, Beijing Normal University, Zhuhai 519087, China

I : The European Pulsar Timing Array

II : The Indian Pulsar Timing Array

Supermassive black hole binaries - full corner plots

Figs. A.1 and A.2 show the full posterior results for the astrophysically-informed and agnostic mode, respectively. Individual parameters are listed in the main text.