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 (Aghanim et al., 2020), and $t_{0}$ is the cosmic time today. The
sum is over the cusp, kink and kink-kink contributions (labelled with the
index $b$) for which $q_{b}=4/3$, $5/3$ and $2$ respectively.

####  4.2.2 Analysis results

Our analysis follows the one presented in Leclere et al. (2023), which was
based on six pulsars with the best timing precision from the EPTA early Data
Release 2 (Chen et al., 2021). We now analyse the DR2new dataset with $25$
pulsars, using the (computationally heavy) hierarchical likelihood data
analysis method described in PaperIII. We sample the SGWB parameters $(G\mu)$
or $(G\mu,N_{k})$ as well as the individual pulsar noise parameters.

Solid lines in Fig. 16 show the posterior distribution on $\log_{10}(G\mu)$ in
case (i), for which $N_{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 Fig. 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 Fig. 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 (Leclere et al., 2023), this is no longer true
for DR2new. This is because the predicted SGWB PSD is generally steeper than
the measured correlated red signal in the data, as can be seen in Fig. 17.

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

Figure 16: Comparison of the string tension posteriors for the two string
models (BOS and LRS) in case (i), $N_{c}=2$ and $N_{k}=0$ ($\Gamma=57$). Solid
lines assume only a cosmic string background, dashed lines assume both a
population of GW-driven circular SMBHBs and cosmic strings. Figure 17: SGWB
spectra (in terms of $\log_{10}h^{2}\Omega_{gw}$) for four different early
Universe SGWB models considered in this paper. BOS/LRS correspond to a cosmic
string background with $N_{c}=2$ and $N_{k}=0$ ($\Gamma=57$), and
$\log_{10}G\mu=-10.1$/$-10.6$. The GWB from turbulence is plotted in solid
line for $\lambda_{*}\mathcal{H}_{*}=1$, $\Omega_{*}=0.3$, and $T_{*}=140$
MeV. The inflationary spectra is shown for $\log_{10}r=-13.1$ and $n_{T}=2.4$
(maximum a posteriori value). Power spectrum of the 2nd-order GWB from the
scalar curvature perturbations described by the powerlaw model with
$A_{\zeta}^{\text{10yr}}=-2.9$ and $n_{s}=2.1$ is shown with brown puncture-
dot line. The nine first Fourier bins posteriors of the common signal are
represented by the grey violin areas.

###  4.3 Implications on background from turbulence around the QCD energy
scale

Turbulence can arise in the early Universe in the aftermath of a first-order
phase transition (Witten, 1984; Kamionkowski et al., 1994), or can be driven
by pre-existing primordial magnetic fields (Quashnock et al., 1989;
Brandenburg et al., 1996). If the (magneto-)hydrodynamic turbulence were
present around the QCD epoch, when the Universe had a temperature of
$T_{*}\sim 100$ MeV, it would generate a GWB in the PTA band. The
characteristic scale of the turbulence, determining the characteristic GW
frequency, is in fact related to the (comoving) Hubble radius at that epoch
$\lambda_{*}\simeq\mathcal{O}(\mathcal{H}_{*}^{-1})$, where

${\cal H}_{*}\simeq
10^{-8}\frac{T_{*}}{100{\rm\,MeV}}\biggl{(}\frac{g_{*}}{10}\biggr{)}^{1\over
6}\,{\rm\,Hz},$ (20)

and $g_{*}$ denotes the number of relativistic degrees of freedom. If a large
lepton asymmetry and/or primordial magnetic fields were present in the early
Universe, the QCD phase transition might have been of first order (Schwarz &
Stuke, 2009; Wygas et al., 2018; Middeldorf-Wygas et al., 2020, 2022;
Vovchenko et al., 2021; Cao, 2023). In this case, one would expect additional
sources of GWs, from the collision of broken phase bubbles and the subsequent
development of sound waves in the primordial fluid (Kosowsky et al., 1992;
Kosowsky & Turner, 1993; Caprini et al., 2008; Huber & Konstandin, 2008; Jinno
& Takimoto, 2017; Cutting et al., 2018; Hindmarsh et al., 2014, 2015, 2017).
This was analysed for PTAs, e.g. in Moore & Vecchio (2021); Arzoumanian et al.
(2021); Xue et al. (2021). In what follows, we focus on the GWB generated by
decaying (M)HD turbulence.

####  4.3.1 Description of the model

The presence of bulk velocity and magnetic fields produce anisotropic
stresses, which in turn act as a source of GWs (Kamionkowski et al., 1994;
Kosowsky et al., 2002; Dolgov et al., 2002; Caprini & Durrer, 2006;
Gogoberidze et al., 2007; Caprini et al., 2009). This has been recently
studied via numerical simulations in Roper Pol et al. (2020a, b); Brandenburg
et al. (2021); Roper Pol et al. (2022a). In particular, Roper Pol et al.
(2022a) show that the envelope of the GWB produced by decaying MHD turbulence
can be estimated analytically, assuming that the anisotropic stresses from the
velocity and magnetic fields vary more slowly than the dynamical production of
GWs. This was also validated by numerical simulations of purely kinetic
turbulence in Auclair et al. (2022). This assumption leads to the following
GWB signal:

$\Omega_{\rm GW}(f)=3\,{\cal A}\,\Omega_{*}^{2}\,\bigl{(}\lambda_{*}{\cal
H}_{*}\bigr{)}^{2}F_{\rm GW,0}\,S_{\rm turb}(\lambda_{*}f),$ (21)

where $\Omega_{*}$ is the ratio of the (M)HD turbulent energy density to the
radiation one, and $\lambda_{*}{\cal H}_{*}$ is the ratio of the
characteristic length scale of the turbulence, $\lambda_{*}$, to the comoving
Hubble horizon ${\cal H}_{*}^{-1}$ at the QCD epoch. The parameter ${\cal
A}\simeq 1.75\times 10^{-3}$ is the efficiency of GW production,888This
estimate is conservative since it only considers the decaying stage of
turbulence. Numerical simulations find larger values when including a stage of
turbulence production (Roper Pol et al., 2020b; Kahniashvili et al., 2021;
Roper Pol et al., 2022b) estimated in Roper Pol et al. (2022a). 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_{*}$,

$F_{\rm GW,0}\simeq 8\times
10^{-5}\,\biggl{(}\frac{10}{g_{*}}\biggr{)}^{1/3}.$ (22)

The spectral shape of the GWB signal, $S_{\rm turb}(f)$, is

$\displaystyle S_{\rm turb}(\lambda_{*}f)=$ $\displaystyle\,{\cal
B}\,\bigl{(}\lambda_{*}f\bigr{)}^{3}\,p_{\Pi}(\lambda_{*}f)\,$
$\displaystyle\times$
$\displaystyle\,\left\\{\begin{array}[]{ll}\ln^{2}\bigl{[}1+{\cal
H}_{*}\,\delta t_{\rm fin}/(2\pi)\bigr{]},&{\rm\ if\ }f<1/\delta t_{\rm
fin},\\\ \ln^{2}\bigl{[}1+\lambda_{*}{\cal
H}_{*}/(2\pi\lambda_{*}f)\bigr{]},&{\rm\ if\ }f\geq 1/\delta t_{\rm
fin},\end{array}\right.$ (25)

where ${\cal B}\simeq 50\,\bigl{(}\lambda_{*}{\cal H}_{*}\bigr{)}^{-2}$ is a
normalising factor, and $\delta t_{\rm fin}$ denotes the effective duration of
the turbulence. The latter can be estimated, from the numerical simulations
performed in Roper Pol et al. (2022a), to be $\delta t_{\rm fin}\simeq
2\lambda_{*}/\sqrt{1.5\,\Omega_{*}}$. The function $p_{\Pi}(\lambda_{*}f)$ in
Eq. (25) 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\gtrsim 1/\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_{\rm fin}$, $f$
at intermediate frequencies $1/\delta t_{\rm 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_{*}{\cal
H}_{*}$. By causality, $\lambda_{*}{\cal H}_{*}$ is bound to be smaller than
one. In general, also $\Omega_{*}\lesssim 1$, otherwise turbulence would
change the dynamics of the Universe. However, note that the template described
above has been validated in principle only for non-relativistic plasma
motions, for which $\Omega_{*}\lesssim{\cal O}(0.1)$.

####  4.3.2 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\textup{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_{*}\sim 60$ MeV. However, at frequencies
around the peak, the signal corresponds to a power spectral density for the
residuals steeper than $\gamma\sim 4$, which cannot fit the data well. For
this reason, values of $\Omega_{*}\lesssim 0.1$ are disfavoured.

For larger values of $\Omega_{*}$, the $f^{3}$ part of the spectrum at
frequencies below $\delta t_{\rm
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_{\rm fin}^{-1}$ and the spectral peak
at $1/\lambda_{*}$ becomes minimal in the limit $\Omega_{*}\sim 1$. Both of
these characteristics lead to a better fit to the data. This is recovered in
the posteriors of Figure 18, together with the degeneracy between
$\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_{\rm 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.

Figure 18: 2D posteriors for the parameters of the background from turbulence
around the QCD energy scale obtained using a free spectrum fit on DR2new data.
The 68% and 95% credible regions are displayed.

###  4.4 Implications on the 2nd-order GWB produced by primordial curvature
perturbations

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

Figure 20: Results for the power-law model of the curvature perturbations
described by Eq. (27). Left panel: 1$\sigma$ and 2$\sigma$ contours of the
posterior distributions on the amplitude $A_{\zeta}$ and the slope of the
power spectrum $n_{s}$ obtained by the analysis of DR2New. Right panel:
1$\sigma$ and 2$\sigma$ contours of the power spectra inferred from the DR2New
analysis by picking 1000 random samples from the posteriors overlaid with the
current constraints on the primordial power spectrum using the latest Planck
data. The grey colour depicts the 2$\sigma$-confidence intervals. The green
lines and purple shaded areas are the same as in Fig. 19.

It is well-known that scalar, vector and tensor modes of the perturbed metric
do not mix at linear order of the Einstein equations (Lifshitz, 1946; Baumann,
2022). However, scalar curvature perturbations will source propagating
tensorial modes (GWs) at the 2nd order in perturbation theory (Tomita, 1967;
Matarrese et al., 1993, 1998; Noh & Hwang, 2004; Carbone & Matarrese, 2005;
Ananda et al., 2007; Baumann et al., 2007). Such scalar curvature
perturbations and associated primordial density fluctuations inevitably exist
in the Universe and can be directly constrained by observations of the CMB.
The latest Planck data (Planck Collaboration et al., 2020b) suggests that the
power spectrum of the curvature perturbations is nearly scale-invariant with
the amplitude $A_{\zeta}\sim 2\times 10^{-9}$, which implies a marginal
energy-density of the generated GWB. Specifically, when projected to the PTA
sensitivity band, the fractional contribution of the energy density in the
associated GWs becomes $\Omega_{\text{gw}}\sim 10^{-24}$, which is practically
non-detectable by current experiments. On the other hand, some models of
inflation (see, for example, Di & Gong, 2018; Byrnes et al., 2019; Braglia et
al., 2020; Yi & Fei, 2023, and references therein) make it possible to produce
a sharp increase in the power spectrum of primordial curvature perturbations
over many orders of magnitude at small scales.

While the CMB is only capable of directly sampling large cosmological scales
with $k\sim 10^{-3}-10^{-1}$ Mpc-1, small scales stay largely uncovered. PTAs
provide a unique opportunity to complement the CMB measurements by indirectly
probing the scalar curvature perturbations in a scale range $k\sim
10^{6}-10^{8}$ Mpc-1 through the second-order generated GWB, and to place
bounds on the steepest possible growth of the power spectrum as well as
corresponding models of inflation (Saito & Yokoyama, 2009; Bugaev & Klimai,
2011; Chen et al., 2020; Dandoy et al., 2023; Zhao & Wang, 2023).

In this work, we consider two models of the primordial curvature power
spectrum: i) monochromatic and ii) powerlaw. For the former, the primordial
spectrum is modelled as:

$P_{\zeta}=A_{\zeta}\delta(\text{log}k/k^{*})\,$ (26)

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:

$P_{\zeta}=A^{10\text{yr}}_{\zeta}\left(\frac{k}{k_{\text{10yr}}}\right)^{(n_{s}-1)}\,$
(27)

where $n_{s}$ characterises the slope and $k^{\text{10yr}}$ is a normalizing
scale $k_{\text{10yr}}=2\pi/(10yr\times c)$, so that $A^{10\text{yr}}_{\zeta}$
corresponds to dimensionless amplitude at ten years.

In the first scenario, a semi-analytical solution for the induced spectrum of
GWB exists and is given by (Kohri & Terada, 2018; Espinosa et al., 2018):

$\Omega_{\text{GW}}\left(f=\frac{kc}{2\pi}\right)=\frac{3A_{\zeta}^{2}}{64}\left(\frac{4-\tilde{k}^{2}}{4}\right)^{2}\tilde{k}^{2}\left(3\tilde{k}^{2}-2\right)\Theta(2-\tilde{k})\times\\\
\times\left(\pi^{2}(3\tilde{k}^{2}-2)^{2}\Theta(2\sqrt{3}-3\tilde{k})+\left(4+(3\tilde{k}^{2}-2)^{2}\log\left|1-\frac{4}{3\tilde{k}^{2}}\right|\right)^{2}\right)\,$
(28)

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., 2019).

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

$\Omega_{\text{GW}}\left(f=\frac{kc}{2\pi}\right)=Q(n_{s})(A_{\zeta}^{\text{10yr}})^{2}\left(\frac{k}{k_{\text{10yr}}}\right)^{2(n_{s}-1)}\,$
(29)

where $Q(n_{s})$ is the scaling factor which can be evaluated in a range of
$n_{s}$ using interpolation points from Table 1 of Kohri & Terada (2018).

After its production, the GWB is damped due to quantum interactions with the
particles of the primordial plasma at the radiation-dominated epoch, and
redshifted inversely proportionally to the scale factor (as it also occurs to
radiation) starting from the epoch of matter-radiation equality (Saikawa &
Shirai, 2018). The present value of the fractional energy density is then:

$\Omega^{0}_{\text{GW}}=2\Omega_{r}^{0}\left(\frac{g_{*}(T)}{g_{*}(T_{\text{eq}})}\right)\left(\frac{g_{*s}(T)}{g_{*s}(T_{\text{eq}})}\right)^{-\frac{4}{3}}\Omega_{\text{GW}}\,$
(30)

where $T$ is the temperature of the Universe at the moment when structures of
a typical size $1/k$ re-enter the horizon999We conservatively set the
temperature at the epoch of production to 17.35 K., $T_{\text{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:

$S(f)=\frac{H_{0}^{2}}{8\pi^{4}}\frac{\Omega_{\text{GW}}^{0}(f)}{f^{5}},$ (31)

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}^{\text{10yr}}$, uniform in
$[4,12]$ for $\log_{10}(k^{*}/\text{Mpc}^{-1})$, and uniform in $[0.4,2.6]$
for $n_{s}$. Boundaries for the latter are constrained by the limitations of
the numerical approximation of the power law model. For this analysis, we
assumed that the common red noise process detected in the latest EPTA dataset
can be fully explained by the 2nd-order scalar-induced GWs.

Results for monochromatic and power law models are shown in Figures 19 and 20,
respectively. The 2D posterior distribution of the model parameters of the
monochromatic model is depicted with orange contours on the right panel of
Fig. 19. One may notice that the regions of the highest probability are
strongly elongated due to a strong positive correlation between $A_{\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}/\text{Mpc}^{-1})=7.6$ and
$\log_{10}A^{0.05}_{\zeta}=-1.7$ meaning that a whole range of models
predicting Dirac-delta power spectrum can equally good describe the signal of
DR2new. This behaviour is explained in the left panel of Fig. 19, for which we
have simulated GWB signal generated by the monochromatic primordial scalar
perturbations, Eq. (28), 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\sim 3.5$ in the limit of large $k^{*}$, and
becomes consistent with both values obtained with DR2new and the theoretically
predicted 13/3 for the background from circular, GW-driven SMBHBs. This
degeneracy can raise important issues when one tries to disentangle one signal
from another. For the power-law case, the 2D posteriors are shown in the left
panel of Fig. 20 with the following means and 1-$\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 Figs. 19 and 20, we also overplot the inferred
primordial power spectrum with the one obtained from the Planck data. The
orange areas on the right panel of Fig. 20 are 1$\sigma$ and 2$\sigma$
contours of the power-law model described by Eq. 27 reconstructed using 1000
random draws from the posteriors. To explain the observed signatures of the
DR2new in terms of the second-order GWB from the primordial scalar
perturbations, an excess in the primordial spectrum at low scales should be
invoked without violating the CMB inflationary parameters. Such excess has
been proposed in many papers, e.g. in the aforementioned works by Ivanov et
al. (1994); Germani & Prokopec (2017); Di & Gong (2018); Cai et al. (2018);
Biagetti et al. (2018); Byrnes et al. (2019); Motohashi et al. (2020); Braglia
et al. (2020); Yi & Fei (2023). Notably, this power excess would lead to a
copious production of primordial black holes (PBHs) at the radiation-dominated
stage, which is sometimes taken by the authors as the motivation to introduce
them as cold dark matter candidates (Carr et al., 2016, and reference
therein). The PBH formation from cosmological perturbations has been
extensively explored Sasaki et al. (2018). On the radiation-dominated stage,
the PBH mass is related to the mass inside the horizon at the time of the
perturbation entering, $M\sim M_{H}\approx m_{Pl}^{2}t\sim 8M_{\odot}(100{\rm
MeV}/T)^{2}(40/g_{*})^{1/2}$, where $m_{Pl}$ is the Planckian mass, $T$ is the
temperature. In terms of the mode comoving wavenumber at the moment of the
horizon crossing, $k=aH$, the part of the horizon mass collapsing into a PBH
reads

$M(k)\approx
30M_{\odot}\left(\frac{10.75}{g_{*}}\right)^{-1/6}\left(\frac{3\times
10^{5}\rm Mpc^{-1}}{k}\right)^{2}\,.$ (32)

For example, in the two-field inflationary model by Braglia et al. (2020), a
peak around $k\sim 2\times 10^{6}$ Mpc-1 could explain the power
$A_{\zeta}\sim 10^{-2}$ and simultaneously lead to the production of PBHs with
masses peaked at $\sim 35M_{\odot}$ (see also the analysis of the NANOGrav
results in Vaskonen & Veermäe, 2021). Interestingly, such a peak seems to be
observed in the chirp mass distribution of LVK merging binary black holes
(Abbott et al., 2023). A more general list of primordial power spectra, as
well as a careful retranslation of them to the PBH abundance and their mass
function, will be considered in a follow-up paper (Porayko et al., in prep.).

##  5 Implications III: dark matter

Unlike spatially and temporally-correlated stochastic processes discussed in
Secs. 3 and 4, in this Section, we explore a possible deterministic
contribution to the EPTA signal from ultralight scalar-field dark matter
(ULDM). For comparison, the morphology of a putative ULDM signal, predicted by
Khmelnitsky & Rubakov (2014), is similar to a CGW from a SMBHB, that is, it is
prominent only in one frequency bin. Given that the signal observed with
DR2new is mostly apparent in the first two fundamental $T^{-1}$ frequency
bins, it is of interest to consider possible contributions from physical
processes with narrowband spectra. Therefore, the analysis presented here
complements the CGW interpretation of the signal by EPTA and InPTA
Collaborations (2023c). Additionally, the ULDM search with DR2full is
performed in Smarra et al. (2023).

Dark matter currently constitutes approximately 26% of the energy density of
the Universe, as confirmed by, e.g., galactic rotation curves (Rubin & Ford,
1970; Rubin et al., 1980; de Salas et al., 2019), baryonic acoustic
oscillations and cosmic microwave background measurements (Bennett et al.,
2013; Aghanim et al., 2020) as well as galaxy surveys (Escudero et al., 2015).
The standard cold dark matter (CDM) picture, whose leading candidates are the
Weakly Interacting Massive Particles (WIMPs) (Arcadi et al., 2018) and the QCD
axion (Luzio et al., 2020), successfully grasps the large-scale structure of
the Universe. However, it presents some well-known issues, when it comes to
explaining observations at scales smaller than $\mathcal{O}$(kpc). Among
these, the cusp-core problem (Flores & Primack, 1994; Moore, 1994; Karukes, E.
V. et al., 2015) concerns the inconsistency between the observation of a flat
density profile in the centre of galaxies and the power-law-like behavior
predicated by CDM, while the mismatch between the simulated and observed
number of dwarf galaxies in the proximity of our Milky Way is often referred
to as the missing satellite problem (Klypin et al., 1999; Moore, 1994).

While a thorough understanding of baryonic physics feedback mechanisms
(Navarro et al., 1996; Governato et al., 2012; Brooks et al., 2013; Chan et
al., 2015; Oñorbe et al., 2015; Read et al., 2016; Morganti, 2017) might help
to alleviate some of these issues, they can be more easily disposed of
assuming that DM is an ultralight ($m_{\phi}\sim 10^{-22}$eV)
scalar/pseudoscalar or axion-like field, whose astrophysically large
($\mathcal{O}$(kpc)) de Broglie wavelength suppresses power on small scales.
Moreover, ultralight scalars are also generally present in string theory
compactifications (Green et al., 1988; Svrcek & Witten, 2006; Arvanitaki et
al., 2010), which makes them interesting candidates for new physics as well.
CMB-based arguments are used to constrain $m_{\phi}\gtrsim 10^{-24}$eV (Hlozek
et al., 2015), while Lyman-$\alpha$ bounds push the limit up to
$m_{\phi}\gtrsim 10^{-21}$eV, provided that the ultralight particles account
for more than 30% of the full dark matter budget (Iršič et al., 2017;
Armengaud et al., 2017; Kobayashi et al., 2017; Nori et al., 2018; Rogers &
Peiris, 2021). A lower limit of $m_{\phi}\sim 10^{-19}$eV is claimed by
studies of ultra-faint dwarf (UFD) galaxies (Hayashi et al., 2021; Dalal &
Kravtsov, 2022), but a wide consensus is yet to be reached. In a seminal work,
Khmelnitsky & Rubakov (2014) showed that the travel time of pulsar radio beams
is affected by the gravitational potential induced by ULDM particles, making
thus PTAs excellent facilities to investigate the existence of ULDM particles.
Moreover, they represent complementary probes which do not suffer from the
small-scale structure modelling uncertainties that affect non-CMB bounds
(Schive et al., 2014; Zhang et al., 2019), as the aforementioned
Lyman-$\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 pulses as described in Khmelnitsky & Rubakov (2014).

Being non-relativistic and with a very large characteristic occupation number,
the ULDM scalar field can be described as a classical wave whose oscillation
frequency is twice its mass $m_{\phi}$ (Khmelnitsky & Rubakov, 2014). The
periodic displacement that ULDM induces on the TOAs of signal from a pulsar
$P$ can be written as follows (Porayko et al., 2018):

$\delta
t(t)=\frac{\Psi_{c}(\@vec{x})}{2m_{\phi}}[\hat{\phi}^{2}_{E}\sin{(2m_{\phi}+\gamma_{E})}-\hat{\phi}^{2}_{P}\sin{(2m_{\phi}+\gamma_{P})}],$
(33)

where

$\Psi_{c}(\@vec{x})=\frac{\pi G\rho_{\phi}}{m^{2}_{\phi}}\approx 6.52\cdot
10^{-18}\left(\frac{10^{-22}\text{eV}}{m_{\phi}}\right)^{2}\left(\frac{\rho_{\phi}}{0.4\text{GeV}/\text{cm}^{3}}\right),$
(34)

where $\gamma_{P}$ ($\gamma_{E}$) is a pulsar (Earth) dependent phase and
$\hat{\phi}^{2}_{P}$ ($\hat{\phi}^{2}_{E}$) takes into account the stochastic
nature of the axion-like field near the pulsar (Earth). The parameters and
their priors are summarised in Table 1. Considering a typical value of
$v_{\phi}\sim 10^{-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:

$l_{c}\approx\frac{2\pi}{m_{\phi}v_{\phi}}\approx
0.4\text{kpc}\left(\frac{10^{-22}\text{eV}}{m_{\phi}}\right).$ (35)

In particular, notice that $\hat{\phi}^{2}_{E}$ and $\hat{\phi}^{2}_{P}$
should be taken as:

  * •

different parameters when the average pulsar-Earth and pulsar-pulsar distance
is larger than the coherence length.

  * •

the same parameter when the average pulsar-Earth and pulsar-pulsar distance is
smaller than the coherence length.

Following the procedure in Smarra et al. (2023), we analyze three separate
cases, which we refer to as the uncorrelated, the pulsar correlated and the
correlated limit. As the average inter-pulsar and Earth-pulsar separation is
$\sim$ 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
$\sim 20$ kpc), but larger than the average inter-pulsar and pulsar-Earth
distance. More specifically, the correlated regime holds for masses lower than
$m_{\phi}\sim 2\times 10^{-24}$ eV; the pulsar correlated regime for $2\times
10^{-24}\,\text{eV}\lesssim m_{\phi}\lesssim 5\times 10^{-23}\,\text{eV}$ and
the uncorrelated limit for $m_{\phi}\gtrsim 5\times
10^{-23}\leavevmode\nobreak\ \text{eV}$. We defer a more detailed study to
future analysis.

Table 1: Parameters used for the search and their respective priors. In the correlated limit, $\hat{\phi}^{2}_{E}=\hat{\phi}^{2}_{P}$ and can be reabsorbed in a redefinition of $\Psi_{c}$. Parameter | Description | Prior | Occurrence
---|---|---|---
White Noise $\left(\sigma=E_{f}^{2}\sigma^{2}_{ToA}+E_{q}^{2}\right)$
$E_{f}$ | EFAC per backend/receiver system | Uniform $[0,10]$ | 1 per pulsar
$E_{q}$ | EQUAD per backend/receiver system | Log10-Uniform $[-10,-5]$ | 1 per pulsar
Red Noise
$A_{\text{red }}$ | Red noise power-law amplitude | Log10-Uniform $[-20,-6]$ | 1 per pulsar
$\gamma_{\text{red }}$ | Red noise power-law spectral index | Uniform $[0,10]$ | 1 per pulsar
ULDM
$\Psi_{\text{c}}$ | ULDM amplitude | Log10-Uniform $[-20,-12]$ | 1 for PTA
$m_{\phi}\leavevmode\nobreak\ [\mathrm{eV}]$ | ULDM mass | Log10-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 per 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 | Log10-Uniform $[-20,-6]$ | 1 for PTA

Figure 21: Posterior probabilities for the ULDM amplitude $\Psi_{c}$ and mass
$m_{\phi}$, from the correlated (top row) and uncorrelated (bottom row)
analysis of the DR2new dataset. The pulsar correlated analysis is not shown,
but displays the same features.

Based on the above, we fit the model from Eq. (33) to the DR2new as in Smarra
et al. (2023). As an example, marginalised posterior distributions for
$\Psi_{c}$ and $m_{\phi}$ are plotted in Fig. 21. All the three limits peak at
a ULDM mass $m_{\phi}\sim 10^{-23}\text{eV}$ with amplitude $\Psi_{c}\sim
10^{-15}$, which translates into a density $\rho_{\phi}\simeq
0.6\,\text{GeV}/\text{cm}^{3}$ of the scalar field. While this value is higher
than the fiducial local DM density $\rho_{\text{DM}}\approx
0.4\text{GeV}/\text{cm}^{3}$, it is well within the measurement uncertainties
(Bovy & Tremaine, 2012; Read, 2014; Sivertsson et al., 2018; de Salas, 2020).

Additionally, we search for a potential ULDM signature alongside the SMBHB
gravitational-wave background. Thus, we introduce a new model that contains
ULDM contributions alongside a common red signal to account for gravitational
wave contributions. We find no ULDM signal under this hypothesis, in 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: Constraints on $\Psi_{c}$ as a function of $m_{\phi}$ using the
EPTA DR2new dataset from PaperIII. Previous analyses are shown for comparison,
cf. Porayko & Postnov (2014); Porayko et al. (2018) for further details. The
blue, orange and brown lines represent the 95% Bayesian upper limit on
$\Psi_{c}$ obtained from the EPTA DR2new dataset with the correlated,
uncorrelated and pulsar correlated analysis, respectively. The purple line
shows the expected ULDM abundance computed from Eq. (34). Figure 23:
Constraints on the ULDM density $\rho_{\phi}$ normalised to the DM background
value $\rho_{\text{DM}}=0.4\text{GeV}/\text{cm}^{3}$. The blue, orange and
brown lines represent the 95% Bayesian upper limit on $\rho_{\phi}$, obtained
from the EPTA DR2new dataset with the correlated, uncorrelated and pulsar
correlated analysis, respectively. The purple dotted line shows the fiducial
local DM density value.

Figs. 22 and 23 show that the EPTA at current sensitivity is able to constrain
the presence of ULDM at the level of the expected DM abundance in the mass
range $m_{\phi}\sim[10^{-24}\text{eV},10^{-23.2}\text{eV}]$. The DR2Full
analysis performed in Smarra et al. (2023) pushes these limits down, thus
implying that such ULDM candidates cannot constitute the entire DM. We
highlight that our bounds extend below the DR2new sampling frequency
$f=1/T\approx 3\cdot 10^{-9}\text{Hz}$, with $T\approx\text{10yr}$. In fact,
while it might naively be thought that the ULDM field needs to complete an
oscillation in the timescale $\tau_{\text{PTA}}$ of the PTA experiment to
produce a detectable effect, we point out that an ULDM wave with frequency
$m_{\phi}<1/\tau_{\text{PTA}}$ can still be approximated by an expansion in
powers of $m_{\phi}t$ (Kaplan et al., 2022), though the sensitivity of PTAs
will be reduced due to degeneracies with the timing model (Ramani et al.,
2020). Relying on the robust CMB bounds mentioned before, we fix the lower end
of our search at $m_{\phi}=10^{-24}\text{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 Fig. 22, we report a bump in upper
limits at $m_{\phi}\sim 10^{-23}\text{eV}$, which is at around the maximum-a-
posteriori boson mass in Fig. 21. This mass further corresponds to the
frequency of the CGW analyzed in EPTA and InPTA Collaborations (2023c).
Looking at the posteriors in this mass range, we also find an additional
contribution, on top of the CURN process. A similar bump in upper limits is
also present in the DR2Full dataset, as discussed in Smarra et al. (2023).
However, the Bayes factor of ln$\mathcal{B}\sim 0.3$ we find in favour of the
presence of ULDM signal is still inconclusive. We recommend following up on
this bump in future work.

##  6 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 $\approx
60$ or, equivalently, at a $p-$value of $\approx 0.001$, indicative of a
$\gtrsim 3\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_{f=1{\rm
yr}^{-1}}=-13.94_{-0.48}^{+0.23}$ and $\gamma=2.71_{-0.71}^{+1.18}$
(PaperIII).

We considered several physical processes separately, investigating the
implications of the detected signal under the hypothesis that it is generated
by that specific process. Our main findings can be summarised as follows.

SMBHBs. The signal is consistent with a cosmic population of merging SMBHBs.
Phenomenological models based on galaxy pairs observations can account for the
power spectral distribution of the correlated signal. Those models can also be
used to predict the chance to detect CGWs in the current data, for which the
search has given inconclusive results so far (EPTA and InPTA Collaborations,
2023c). According to those models, there is roughly a 50% chance to detect a
CGW in DR2new at S/N$>3$. The relatively high amplitude of the signal can be
used to place constraints on key properties of the cosmic SMBHB population. By
exploiting the inference framework developed in Middleton et al. (2016); Chen
et al. (2019), we can infer that SMBHBs merge in less than 1Gyr following
galaxy mergers, and that the SMBH-stellar bulge relation has a normalization
${\rm log}_{10}(M_{*}/{\rm M}_{\odot})\approx 8.4$ at a reference stellar mass
bulge of $10^{11}$M⊙, in line with recent compiled results (e.g. Kormendy &
Ho, 2013). Finally, we investigated how the detected signal compares with
predictions from state-of-the-art galaxy formation SAMs (Barausse, 2012;
Izquierdo-Villalba et al., 2022). We found that the high amplitude of the
signal imposes 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\approx 2.7$ implies $n_{T}\approx 2.3$ for a radiation-
dominated universe with $w=1/3$. A cosmic string origin would allow narrowing
down the string tension to values of $-11\lesssim{\rm
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 ${\rm log}_{10}(m_{\phi}/{\rm eV})\approx-23$. This corresponds to a
field oscillation frequency of ${\rm log}_{10}(f/{\rm Hz})\approx-8.3$, which
is consistent with the CGW candidate examined in EPTA and InPTA Collaborations
(2023c). When a joint ULDM+GWB search is performed, however, the peak in the
ULDM mass posterior distribution vanishes, as the data strongly prefer the
presence of an HD correlation in the common power. We therefore conclude that
an ULDM origin of the detected signal is disfavoured, placing a direct
constraint on the abundance of ULDM in our Galaxy. The non-detection in DR2new
implies that only about 80% of the DM density in the solar neighbourhood can
be attributed to ULDM with $-24<{\rm log}_{10}(m_{\phi}/{\rm eV})<-23.7$. More
stringent constraints are obtained from DR2full and are presented in a
separate paper (Smarra et al., 2023).

It is interesting to remark that the best fit to the measured power-law slope
is $\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 Sec. 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 Section 4. Another possibility, proposed by
Lieu et al. (2022) is that GWs can lose energy while propagating through the
intergalactic medium via the acceleration of charged particles. Owing to the
much higher thermal velocity of electrons over protons and the presence of an
intergalactic medium magnetic field, a significant fraction of the GWs can be
converted to EM radiation. The ensuing electromagnetic intensity is inversely
proportional to frequency. Partial GW loss via this mechanism could explain a
flatter-than-expected GWB spectrum. The efficiency of such a mechanism depends
on the strength of the intergalactic medium magnetic field.

It is also possible that the observed signal is coming from multiple
overlapping GWBs. In our analyses, we have assumed a single origin for all the
observed signal power in the EPTA+InPTA data, however, an overlap of GWBs of
different origins can cause the spectral shape of the recovered signal to
deviate from the expected value of any single GWB. Searching, disentangling
and identifying the underlying physical processes will be part of the spectral
characterisation moving forward (Moore & Vecchio, 2021; Kaiser et al., 2022).

As time goes by and the PTA experiment improves, the low-frequency GWB will
leave an increasingly distinctive signature in the data. Interpreting all of
the details of this signature will be necessary to understand the nature of
the signal and exploit the full potential of this new window into the
Universe. As mentioned in the introduction, a population of SMBHBs is expected
to generate a highly non-Gaussian, partially anisotropic and perhaps non-
stationary signal. As we increase the timespan of our data, improve our
instrumentation and combine more pulsars, the detailed properties of the
signal will eventually reveal themselves in the data (e.g. Cornish & Sesana,
2013; Taylor & Gair, 2013; Taylor et al., 2020; Pol et al., 2022). While many
early Universe signals are expected to be isotropic and Gaussian, noticeable
exceptions exist, such as GWBs from bursts of cosmic strings. Cross-
correlating the power distribution of the nanohertz GW sky to the distribution
of massive galaxies and large-scale structures in the low-$z$ universe can
eventually provide the key to determining the true nature – astrophysical vs.
early Universe – of this signal (Rosado & Sesana, 2014; Mingarelli et al.,
2017). In this respect, combining all of the available high-quality datasets
within the IPTA framework is the next step towards the fulfilment of the
promises of nanohertz GW science.

###### 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 Paris (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 and 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 H2020 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. JPWV acknowledges support by the Deutsche
Forschungsgemeinschaft (DFG) through thew Heisenberg programme (Project No.
433075039) and by the NSF through AccelNet award #2114721. NKP is funded by
the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) –
Projektnummer PO 2758/1–1, through the Walter–Benjamin programme. ASa thanks
the Alexander von Humboldt foundation in Germany for a Humboldt fellowship for
postdoctoral researchers. APo, DP and MBu acknowledge support from the
research grant “iPeska” (P.I. Andrea Possenti) funded under the INAF national
call Prin-SKA/CTA 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, 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 Starobinskiy, Sergei
Blinnikov and Alexander Dolgov for discussions on the early Universe physics.
HM 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
CRediT101010https://credit.niso.org/ 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.

## References

  * Abbott et al. (2017) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017, Phys. Rev. Lett., 118, 121101
  * Abbott et al. (2016) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016, Physical Review Letters, 116, 061102
  * Abbott et al. (2018) Abbott, B. P. et al. 2018, Phys. Rev. D, 97, 102002
  * Abbott & Wise (1984) Abbott, L. F. & Wise, M. B. 1984, Nuclear Physics B, 244, 541
  * Abbott et al. (2023) Abbott, R., Abbott, T. D., Acernese, F., et al. 2023, Physical Review X, 13, 011048
  * Abbott et al. (2021) Abbott, R. et al. 2021, Phys. Rev. Lett., 126, 241102
  * Aghanim et al. (2020) Aghanim, N. et al. 2020, Astron. Astrophys., 641, A6, [Erratum: Astron.Astrophys. 652, C4 (2021)]
  * Allen (1996) Allen, B. 1996, in Les Houches School of Physics: Astrophysical Sources of Gravitational Radiation, 373–417
  * Allen & Koranda (1994) Allen, B. & Koranda, S. 1994, Phys. Rev. D, 50, 3713
  * Amaro-Seoane et al. (2010) Amaro-Seoane, P., Sesana, A., Hoffman, L., et al. 2010, MNRAS, 402, 2308
  * Ananda et al. (2007) Ananda, K. N., Clarkson, C., & Wands, D. 2007, Phys. Rev. D, 75, 123518
  * Anber & Sorbo (2012) Anber, M. M. & Sorbo, L. 2012, Phys. Rev. D, 85, 123537
  * Anber & Sorbo (2012) Anber, M. M. & Sorbo, L. 2012, Phys. Rev. D, 85, 123537
  * Antonini et al. (2015) Antonini, F., Barausse, E., & Silk, J. 2015, ApJ, 806, L8
  * Arcadi et al. (2018) Arcadi, G., Dutra, M., Ghosh, P., et al. 2018, The European Physical Journal C, 78
  * Armengaud et al. (2017) Armengaud, E., Palanque-Delabrouille, N., Yèche, C., Marsh, D. J. E., & Baur, J. 2017, Monthly Notices of the Royal Astronomical Society, 471, 4606
  * Arvanitaki et al. (2010) Arvanitaki, A., Dimopoulos, S., Dubovsky, S., Kaloper, N., & March-Russell, J. 2010, Physical Review D, 81
  * Arzoumanian et al. (2020) Arzoumanian, Z., Baker, P. T., Blumer, H., et al. 2020, ApJ, 905, L34
  * Arzoumanian et al. (2016) Arzoumanian, Z., Brazier, A., Burke-Spolaor, S., et al. 2016, ApJ, 821, 13
  * Arzoumanian et al. (2021) Arzoumanian, Z. et al. 2021, Phys. Rev. Lett., 127, 251302
  * Auclair et al. (2023a) Auclair, P., Babak, S., Quelquejay Leclere, H., & Steer, D. A. 2023a [arXiv:2305.11653]
  * Auclair et al. (2022) Auclair, P., Caprini, C., Cutting, D., et al. 2022, JCAP, 09, 029
  * Auclair & Ringeval (2022) Auclair, P. & Ringeval, C. 2022, Phys. Rev. D, 106, 063512
  * Auclair et al. (2023b) Auclair, P., Steer, D. A., & Vachaspati, T. 2023b [arXiv:2306.08331]
  * Auclair et al. (2020) Auclair, P. et al. 2020, JCAP, 04, 034
  * Auclair (2020) Auclair, P. G. 2020, JCAP, 11, 050
  * Babak & Sesana (2012) Babak, S. & Sesana, A. 2012, Phys. Rev. D, 85, 044034
  * Barausse (2012) Barausse, E. 2012, MNRAS, 423, 2533
  * Barausse et al. (2020) Barausse, E., Dvorkin, I., Tremmel, M., Volonteri, M., & Bonetti, M. 2020, ApJ, 904, 16
  * Barnaby et al. (2012) Barnaby, N., Moxon, J., Namba, R., et al. 2012, Phys. Rev. D, 86, 103508
  * Bartolo et al. (2016) Bartolo, N., Caprini, C., Domcke, V., et al. 2016, J. Cosmology Astropart. Phys., 2016, 026
  * Bartolo et al. (2007) Bartolo, N., Matarrese, S., Riotto, A., & Väihkönen, A. 2007, Phys. Rev. D, 76, 061302
  * Bassa et al. (2016) Bassa, C. G., Janssen, G. H., Karuppusamy, R., et al. 2016, 456, 2196
  * Baumann (2022) Baumann, D. 2022, Cosmology (Cambridge University Press)
  * Baumann et al. (2007) Baumann, D., Steinhardt, P., Takahashi, K., & Ichiki, K. 2007, Phys. Rev. D, 76, 084019
  * Bécsy et al. (2022) Bécsy, B., Cornish, N. J., & Kelley, L. Z. 2022, ApJ, 941, 119
  * Bennett et al. (2013) Bennett, C. L., Larson, D., Weiland, J. L., et al. 2013, The Astrophysical Journal Supplement Series, 208, 20
  * Biagetti et al. (2015) Biagetti, M., Dimastrogiovanni, E., Fasiello, M., & Peloso, M. 2015, J. Cosmology Astropart. Phys., 2015, 011
  * Biagetti et al. (2013) Biagetti, M., Fasiello, M., & Riotto, A. 2013, Phys. Rev. D, 88, 103518
  * Biagetti et al. (2018) Biagetti, M., Franciolini, G., Kehagias, A., & Riotto, A. 2018, J. Cosmology Astropart. Phys., 2018, 032
  * Bian et al. (2022) Bian, L., Shu, J., Wang, B., Yuan, Q., & Zong, J. 2022, Phys. Rev. D, 106, L101301
  * Blanco-Pillado & Olum (2017) Blanco-Pillado, J. J. & Olum, K. D. 2017, Phys. Rev. D, 96, 104046
  * Blanco-Pillado et al. (2014) Blanco-Pillado, J. J., Olum, K. D., & Shlaer, B. 2014, Phys. Rev. D, 89, 023512
  * Blanco-Pillado et al. (2015) Blanco-Pillado, J. J., Olum, K. D., & Shlaer, B. 2015, Phys. Rev. D, 92, 063528
  * Blasi et al. (2021) Blasi, S., Brdar, V., & Schmitz, K. 2021, Phys. Rev. Lett., 126, 041305
  * Bonetti & Sesana (2020) Bonetti, M. & Sesana, A. 2020, Phys. Rev. D, 102, 103023
  * Bonetti et al. (2018a) Bonetti, M., Sesana, A., Barausse, E., & Haardt, F. 2018a, MNRAS, 477, 2599
  * Bonetti et al. (2018b) Bonetti, M., Sesana, A., Barausse, E., & Haardt, F. 2018b, MNRAS, 477, 2599
  * Bovy & Tremaine (2012) Bovy, J. & Tremaine, S. 2012, The Astrophysical Journal, 756, 89
  * Boylan-Kolchin et al. (2009) Boylan-Kolchin, M., Springel, V., White, S. D. M., Jenkins, A., & Lemson, G. 2009, MNRAS, 398, 1150
  * Boyle & Buonanno (2008) Boyle, L. A. & Buonanno, A. 2008, Phys. Rev. D, 78, 043531
  * Boyle & Steinhardt (2008) Boyle, L. A. & Steinhardt, P. J. 2008, Phys. Rev. D, 77, 063504
  * Braglia et al. (2020) Braglia, M., Hazra, D. K., Finelli, F., et al. 2020, J. Cosmology Astropart. Phys., 2020, 001
  * Brandenburg et al. (2021) Brandenburg, A., Clarke, E., He, Y., & Kahniashvili, T. 2021, Phys. Rev. D, 104, 043513
  * Brandenburg et al. (1996) Brandenburg, A., Enqvist, K., & Olesen, P. 1996, Phys. Rev. D, 54, 1291
  * Brooks et al. (2013) Brooks, A. M., Kuhlen, M., Zolotov, A., & Hooper, D. 2013, The Astrophysical Journal, 765, 22
  * Bugaev & Klimai (2011) Bugaev, E. & Klimai, P. 2011, Phys. Rev. D, 83, 083521
  * Byrnes et al. (2019) Byrnes, C. T., Cole, P. S., & Patil, S. P. 2019, J. Cosmology Astropart. Phys., 2019, 028
  * Cai et al. (2018) Cai, Y.-F., Tong, X., Wang, D.-G., & Yan, S.-F. 2018, Phys. Rev. Lett., 121, 081306
  * Cao (2023) Cao, G. 2023, Phys. Rev. D, 107, 014021
  * Capelo et al. (2017) Capelo, P. R., Dotti, M., Volonteri, M., et al. 2017, MNRAS, 469, 4437
  * Caprini & Durrer (2006) Caprini, C. & Durrer, R. 2006, Phys. Rev. D, 74, 063521
  * Caprini et al. (2008) Caprini, C., Durrer, R., & Servant, G. 2008, Phys. Rev. D, 77, 124015
  * Caprini et al. (2009) Caprini, C., Durrer, R., & Servant, G. 2009, J. Cosmology Astropart. Phys., 2009, 024
  * Caprini & Figueroa (2018) Caprini, C. & Figueroa, D. G. 2018, Classical and Quantum Gravity, 35, 163001
  * Carbone & Matarrese (2005) Carbone, C. & Matarrese, S. 2005, Phys. Rev. D, 71, 043508
  * Carr et al. (2016) Carr, B., Kühnel, F., & Sandstad, M. 2016, Phys. Rev. D, 94, 083504
  * Chan et al. (2015) Chan, T. K., Kereš, D., Oñorbe, J., et al. 2015, Monthly Notices of the Royal Astronomical Society, 454, 2981
  * Chen et al. (2021) Chen, S., Caballero, R. N., Guo, Y. J., et al. 2021, MNRAS, 508, 4970
  * Chen et al. (2017a) Chen, S., Middleton, H., Sesana, A., Del Pozzo, W., & Vecchio, A. 2017a, MNRAS, 468, 404
  * Chen et al. (2019) Chen, S., Sesana, A., & Conselice, C. J. 2019, MNRAS, 488, 401
  * Chen et al. (2017b) Chen, S., Sesana, A., & Del Pozzo, W. 2017b, MNRAS, 470, 1738
  * Chen et al. (2022) Chen, Z.-C., Wu, Y.-M., & Huang, Q.-G. 2022, Astrophys. J., 936, 20
  * Chen et al. (2020) Chen, Z.-C., Yuan, C., & Huang, Q.-G. 2020, Phys. Rev. Lett., 124, 251101
  * Chluba et al. (2012) Chluba, J., Erickcek, A. L., & Ben-Dayan, I. 2012, ApJ, 758, 76
  * Cook & Sorbo (2013) Cook, J. L. & Sorbo, L. 2013, J. Cosmology Astropart. Phys., 2013, 047
  * Cook & Sorbo (2013) Cook, J. L. & Sorbo, L. 2013, JCAP, 11, 047
  * Cornish & Sesana (2013) Cornish, N. J. & Sesana, A. 2013, Classical and Quantum Gravity, 30, 224005
  * Cutting et al. (2018) Cutting, D., Hindmarsh, M., & Weir, D. J. 2018, Phys. Rev. D, 97, 123513
  * Dalal & Kravtsov (2022) Dalal, N. & Kravtsov, A. 2022, Not so fuzzy: excluding FDM with sizes and stellar kinematics of ultra-faint dwarf galaxies
  * Damour & Vilenkin (2000) Damour, T. & Vilenkin, A. 2000, Phys. Rev. Lett., 85, 3761
  * Damour & Vilenkin (2001) Damour, T. & Vilenkin, A. 2001, Phys. Rev. D, 64, 064008
  * Damour & Vilenkin (2001) Damour, T. & Vilenkin, A. 2001, Phys. Rev. D, 64, 064008
  * Damour & Vilenkin (2005) Damour, T. & Vilenkin, A. 2005, Phys. Rev. D, 71, 063510
  * Dandoy et al. (2023) Dandoy, V., Domcke, V., & Rompineve, F. 2023, arXiv e-prints, arXiv:2302.07901
  * de Salas et al. (2019) de Salas, P., Malhan, K., Freese, K., Hattori, K., & Valluri, M. 2019, Journal of Cosmology and Astroparticle Physics, 2019, 037
  * de Salas (2020) de Salas, P. F. 2020, Journal of Physics: Conference Series, 1468, 012020
  * Dehnen (2014) Dehnen, W. 2014, Computational Astrophysics and Cosmology, 1, 1
  * Desvignes et al. (2016) Desvignes, G., Caballero, R. N., Lentati, L., et al. 2016, MNRAS, 458, 3341
  * Di & Gong (2018) Di, H. & Gong, Y. 2018, J. Cosmology Astropart. Phys., 2018, 007
  * Dolgov et al. (2002) Dolgov, A. D., Grasso, D., & Nicolis, A. 2002, Phys. Rev. D, 66, 103505
  * D’Orazio & Duffell (2021) D’Orazio, D. J. & Duffell, P. C. 2021, ApJ, 914, L21
  * Dvali & Vilenkin (2004) Dvali, G. & Vilenkin, A. 2004, J. Cosmology Astropart. Phys., 2004, 010
  * Ellis & Lewicki (2021) Ellis, J. & Lewicki, M. 2021, Phys. Rev. Lett., 126, 041304
  * Ellis et al. (2020) Ellis, J. A., Vallisneri, M., Taylor, S. R., & Baker, P. T. 2020, ENTERPRISE: Enhanced Numerical Toolbox Enabling a Robust PulsaR Inference SuitE, Zenodo
  * Enoki & Nagashima (2007) Enoki, M. & Nagashima, M. 2007, Progress of Theoretical Physics, 117, 241
  * EPTA and InPTA Collaboration (2023) EPTA and InPTA Collaboration. 2023, submitted to A&A
  * EPTA and InPTA Collaborations (2023a) EPTA and InPTA Collaborations. 2023a, submitted to A&A
  * EPTA and InPTA Collaborations (2023b) EPTA and InPTA Collaborations. 2023b, A&A, 678, A50
  * EPTA and InPTA Collaborations (2023c) EPTA and InPTA Collaborations. 2023c, submitted to A&A
  * Escudero et al. (2015) Escudero, M., Mena, O., Vincent, A. C., Wilkinson, R. J., & Bœhm, C. 2015, Journal of Cosmology and Astroparticle Physics, 2015, 034
  * Espinosa et al. (2018) Espinosa, J. R., Racco, D., & Riotto, A. 2018, J. Cosmology Astropart. Phys., 2018, 012
  * Fabbri & Pollock (1983) Fabbri, R. & Pollock, M. D. 1983, Physics Letters B, 125, 445
  * Farris et al. (2014) Farris, B. D., Duffell, P., MacFadyen, A. I., & Haiman, Z. 2014, ApJ, 783, 134
  * Flores & Primack (1994) Flores, R. A. & Primack, J. R. 1994, ApJ, 427, L1
  * Foster & Backer (1990) Foster, R. S. & Backer, D. C. 1990, ApJ, 361, 300
  * Fujita et al. (2015) Fujita, T., Yokoyama, J., & Yokoyama, S. 2015, Progress of Theoretical and Experimental Physics, 2015, 043E01
  * Galloni et al. (2023) Galloni, G., Bartolo, N., Matarrese, S., et al. 2023, J. Cosmology Astropart. Phys., 2023, 062
  * Germani & Prokopec (2017) Germani, C. & Prokopec, T. 2017, Physics of the Dark Universe, 18, 6
  * Giarè et al. (2023) Giarè, W., Forconi, M., Di Valentino, E., & Melchiorri, A. 2023, MNRAS, 520, 1757
  * Giarè & Melchiorri (2021) Giarè, W. & Melchiorri, A. 2021, Physics Letters B, 815, 136137
  * Giovannini (1998) Giovannini, M. 1998, Phys. Rev. D, 58, 083504
  * Gogoberidze et al. (2007) Gogoberidze, G., Kahniashvili, T., & Kosowsky, A. 2007, Phys. Rev. D, 76, 083002
  * Goncharov et al. (2021) Goncharov, B., Shannon, R. M., Reardon, D. J., et al. 2021, On the Evidence for a Common-spectrum Process in the Search for the Nanohertz Gravitational-wave Background with the Parkes Pulsar Timing Array
  * Goncharov et al. (2022) Goncharov, B., Thrane, E., Shannon, R. M., et al. 2022, ApJ, 932, L22
  * Governato et al. (2012) Governato, F., Zolotov, A., Pontzen, A., et al. 2012, Monthly Notices of the Royal Astronomical Society, 422, 1231
  * Graham & Scott (2013) Graham, A. W. & Scott, N. 2013, ApJ, 764, 151
  * Green et al. (1988) Green, M. B., Schwarz, J. H., & Witten, E. 1988, SUPERSTRING THEORY. VOL. 1: INTRODUCTION, Cambridge Monographs on Mathematical Physics
  * Grishchuk (1975) Grishchuk, L. P. 1975, Soviet Journal of Experimental and Theoretical Physics, 40, 409
  * Gualandris et al. (2022) Gualandris, A., Khan, F. M., Bortolas, E., et al. 2022, MNRAS, 511, 4753
  * Hayashi et al. (2021) Hayashi, K., Ferreira, E. G. M., & Chan, H. Y. J. 2021, The Astrophysical Journal Letters, 912, L3
  * Hellings & Downs (1983) Hellings, R. W. & Downs, G. S. 1983, ApJ, 265, L39
  * Henriques et al. (2015) Henriques, B. M. B., White, S. D. M., Thomas, P. A., et al. 2015, MNRAS, 451, 2663
  * Hindmarsh et al. (2014) Hindmarsh, M., Huber, S. J., Rummukainen, K., & Weir, D. J. 2014, Phys. Rev. Lett., 112, 041301
  * Hindmarsh et al. (2015) Hindmarsh, M., Huber, S. J., Rummukainen, K., & Weir, D. J. 2015, Phys. Rev. D, 92, 123009
  * Hindmarsh et al. (2017) Hindmarsh, M., Huber, S. J., Rummukainen, K., & Weir, D. J. 2017, Phys. Rev. D, 96, 103520
  * Hindmarsh & Kibble (1995) Hindmarsh, M. B. & Kibble, T. W. B. 1995, Rept. Prog. Phys., 58, 477
  * Hlozek et al. (2015) Hlozek, R., Grin, D., Marsh, D. J. E., & Ferreira, P. G. 2015, Phys. Rev. D, 91, 103512
  * Horndeski (1974) Horndeski, G. W. 1974, International Journal of Theoretical Physics, 10, 363
  * Huber & Konstandin (2008) Huber, S. J. & Konstandin, T. 2008, J. Cosmology Astropart. Phys., 2008, 022
  * Iršič et al. (2017) Iršič, V., Viel, M., Haehnelt, M. G., Bolton, J. S., & Becker, G. D. 2017, Phys. Rev. Lett., 119, 031302
  * Ivanov et al. (1994) Ivanov, P., Naselsky, P., & Novikov, I. 1994, Phys. Rev. D, 50, 7173
  * Izquierdo-Villalba et al. (2022) Izquierdo-Villalba, D., Sesana, A., Bonoli, S., & Colpi, M. 2022, MNRAS, 509, 3488
  * Jaffe & Backer (2003) Jaffe, A. H. & Backer, D. C. 2003, ApJ, 583, 616
  * Jeannerot et al. (2003) Jeannerot, R., Rocher, J., & Sakellariadou, M. 2003, Phys. Rev. D, 68, 103514
  * Jinno & Takimoto (2017) Jinno, R. & Takimoto, M. 2017, Phys. Rev. D, 95, 024009
  * Jones et al. (2003) Jones, N. T., Stoica, H., & Tye, S. H. H. 2003, Physics Letters B, 563, 6
  * Joshi et al. (2022) Joshi, B. C., Gopakumar, A., Pandian, A., et al. 2022, Journal of Astrophysics and Astronomy, 43, 98
  * Kahniashvili et al. (2021) Kahniashvili, T., Brandenburg, A., Gogoberidze, G., Mandal, S., & Roper Pol, A. 2021, Phys. Rev. Res., 3, 013193
  * Kaiser et al. (2022) Kaiser, A. R., Pol, N. S., McLaughlin, M. A., et al. 2022, ApJ, 938, 115
  * Kamionkowski et al. (1994) Kamionkowski, M., Kosowsky, A., & Turner, M. S. 1994, Phys. Rev. D, 49, 2837
  * Kaplan et al. (2022) Kaplan, D. E., Mitridate, A., & Trickle, T. 2022, Physical Review D, 106
  * Karukes, E. V. et al. (2015) Karukes, E. V., Salucci, P., & Gentile, G. 2015, A&A, 578, A13
  * Kelley et al. (2017) Kelley, L. Z., Blecha, L., Hernquist, L., Sesana, A., & Taylor, S. R. 2017, MNRAS, 471, 4508
  * Kelley et al. (2018) Kelley, L. Z., Blecha, L., Hernquist, L., Sesana, A., & Taylor, S. R. 2018, MNRAS, 477, 964
  * Khan et al. (2016) Khan, F. M., Fiacconi, D., Mayer, L., Berczik, P., & Just, A. 2016, ApJ, 828, 73
  * Khan et al. (2012) Khan, F. M., Preto, M., Berczik, P., et al. 2012, ApJ, 749, 147
  * Khmelnitsky & Rubakov (2014) Khmelnitsky, A. & Rubakov, V. 2014, Journal of Cosmology and Astroparticle Physics, 2014, 019
  * Kibble (1976) Kibble, T. 1976, J. Phys. A, 9, 1387
  * Klein et al. (2016) Klein, A., Barausse, E., Sesana, A., et al. 2016, Phys. Rev. D, 93, 024003
  * Klypin et al. (1999) Klypin, A., Kravtsov, A. V., Valenzuela, O., & Prada, F. 1999, The Astrophysical Journal, 522, 82
  * Kobayashi et al. (2017) Kobayashi, T., Murgia, R., Simone, A. D., Irš ič, V., & Viel, M. 2017, Physical Review D, 96
  * Kocsis & Sesana (2011) Kocsis, B. & Sesana, A. 2011, MNRAS, 411, 1467
  * Kohri & Terada (2018) Kohri, K. & Terada, T. 2018, Phys. Rev. D, 97, 123532
  * Kormendy & Ho (2013) Kormendy, J. & Ho, L. C. 2013, ARA&A, 51, 511
  * Kosowsky (1996) Kosowsky, A. 1996, Annals of Physics, 246, 49
  * Kosowsky et al. (2002) Kosowsky, A., Mack, A., & Kahniashvili, T. 2002, Phys. Rev. D, 66, 024030
  * Kosowsky & Turner (1993) Kosowsky, A. & Turner, M. S. 1993, Phys. Rev. D, 47, 4372
  * Kosowsky et al. (1992) Kosowsky, A., Turner, M. S., & Watkins, R. 1992, Phys. Rev. D, 45, 4514
  * Koss et al. (2018) Koss, M. J., Blecha, L., Bernhard, P., et al. 2018, Nature, 563, 214
  * Kulier et al. (2015) Kulier, A., Ostriker, J. P., Natarajan, P., Lackner, C. N., & Cen, R. 2015, ApJ, 799, 178
  * Lamb et al. (2023) Lamb, W. G., Taylor, S. R., & van Haasteren, R. 2023, The Need For Speed: Rapid Refitting Techniques for Bayesian Spectral Characterization of the Gravitational Wave Background Using PTAs
  * Lasky et al. (2016) Lasky, P. D., Mingarelli, C. M. F., Smith, T. L., et al. 2016, Physical Review X, 6, 011035
  * Leclere et al. (2023) Leclere, H. Q. et al. 2023 [arXiv:2306.12234]
  * Lee (2016) Lee, K. J. 2016, in Astronomical Society of the Pacific Conference Series, Vol. 502, Frontiers in Radio Astronomy and FAST Early Sciences Symposium 2015, ed. L. Qain & D. Li, 19
  * Lentati et al. (2015) Lentati, L., Taylor, S. R., Mingarelli, C. M. F., et al. 2015, MNRAS, 453, 2576
  * Lieu et al. (2022) Lieu, R., Lackeos, K., & Zhang, B. 2022, Classical and Quantum Gravity, 39, 075014
  * Lifshitz (1946) Lifshitz, E. M. 1946, Zhurnal Eksperimentalnoi i Teoreticheskoi Fiziki, 16, 587
  * Lorenz et al. (2010) Lorenz, L., Ringeval, C., & Sakellariadou, M. 2010, JCAP, 10, 003
  * Luzio et al. (2020) Luzio, L. D., Giannotti, M., Nardi, E., & Visinelli, L. 2020, Physics Reports, 870, 1
  * Maggiore (2000) Maggiore, M. 2000, Phys. Rept., 331, 283
  * Manchester et al. (2013) Manchester, R. N., Hobbs, G., Bailes, M., et al. 2013, PASA, 30, e017
  * Matarrese et al. (1998) Matarrese, S., Mollerach, S., & Bruni, M. 1998, Phys. Rev. D, 58, 043504
  * Matarrese et al. (1993) Matarrese, S., Pantano, O., & Saez, D. 1993, Phys. Rev. D, 47, 1311
  * McConnell & Ma (2013) McConnell, N. J. & Ma, C.-P. 2013, ApJ, 764, 184
  * McConnell et al. (2011) McConnell, N. J., Ma, C.-P., Gebhardt, K., et al. 2011, Nature, 480, 215
  * McLaughlin (2013) McLaughlin, M. A. 2013, Classical and Quantum Gravity, 30, 224008
  * McWilliams et al. (2014) McWilliams, S. T., Ostriker, J. P., & Pretorius, F. 2014, ApJ, 789, 156
  * Middeldorf-Wygas et al. (2020) Middeldorf-Wygas, M. M., Oldengott, I. M., Bödeker, D., & Schwarz, D. J. 2020, arXiv e-prints, arXiv:2009.00036
  * Middeldorf-Wygas et al. (2022) Middeldorf-Wygas, M. M., Oldengott, I. M., Bödeker, D., & Schwarz, D. J. 2022, Phys. Rev. D, 105, 123533
  * Middleton et al. (2016) Middleton, H., Del Pozzo, W., Farr, W. M., Sesana, A., & Vecchio, A. 2016, MNRAS, 455, L72
  * Middleton et al. (2021) Middleton, H., Sesana, A., Chen, S., et al. 2021, MNRAS, 502, L99
  * Miles et al. (2023) Miles, M. T., Shannon, R. M., Bailes, M., et al. 2023, MNRAS, 519, 3976
  * Milosavljević & Merritt (2003) Milosavljević, M. & Merritt, D. 2003, ApJ, 596, 860
  * Mingarelli et al. (2017) Mingarelli, C. M. F., Lazio, T. J. W., Sesana, A., et al. 2017, Nature Astronomy, 1, 886
  * Moore (1994) Moore, B. 1994, Nature, 370, 629
  * Moore et al. (1999) Moore, B., Ghigna, S., Governato, F., et al. 1999, ApJ, 524, L19
  * Moore & Vecchio (2021) Moore, C. J. & Vecchio, A. 2021, Nature Astronomy, 5, 1268
  * Morganti (2017) Morganti, R. 2017, Frontiers in Astronomy and Space Sciences, 4
  * Motohashi et al. (2020) Motohashi, H., Mukohyama, S., & Oliosi, M. 2020, J. Cosmology Astropart. Phys., 2020, 002
  * Nasim et al. (2020) Nasim, I., Gualandris, A., Read, J., et al. 2020, MNRAS, 497, 739
  * Nasim et al. (2021) Nasim, I. T., Gualandris, A., Read, J. I., et al. 2021, MNRAS, 502, 4794
  * Navarro et al. (1996) Navarro, J. F., Eke, V. R., & Frenk, C. S. 1996, Monthly Notices of the Royal Astronomical Society, 283, L72
  * Noh & Hwang (2004) Noh, H. & Hwang, J.-C. 2004, Phys. Rev. D, 69, 104011
  * Nori et al. (2018) Nori, M., Murgia, R., Iršič, V., Baldi, M., & Viel, M. 2018, Monthly Notices of the Royal Astronomical Society, 482, 3227
  * Oñorbe et al. (2015) Oñorbe, J., Boylan-Kolchin, M., Bullock, J. S., et al. 2015, Monthly Notices of the Royal Astronomical Society, 454, 2092
  * Perera et al. (2019) Perera, B. B. P., DeCesar, M. E., Demorest, P. B., et al. 2019, MNRAS, 490, 4666
  * Phinney (2001) Phinney, E. S. 2001, arXiv e-prints, astro
  * Pillepich et al. (2018) Pillepich, A., Springel, V., Nelson, D., et al. 2018, MNRAS, 473, 4077
  * Planck Collaboration et al. (2014) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2014, A&A, 571, A16
  * Planck Collaboration et al. (2020a) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020a, A&A, 641, A6
  * Planck Collaboration et al. (2020b) Planck Collaboration, Akrami, Y., Arroja, F., et al. 2020b, A&A, 641, A10
  * Pol et al. (2022) Pol, N., Taylor, S. R., & Romano, J. D. 2022, ApJ, 940, 173
  * Porayko & Postnov (2014) Porayko, N. & Postnov, K. 2014, Physical Review D, 90
  * Porayko et al. (2018) Porayko, N. K., Zhu, X., Levin, Y., et al. 2018, Physical Review D, 98
  * Preto et al. (2011) Preto, M., Berentzen, I., Berczik, P., & Spurzem, R. 2011, ApJ, 732, L26
  * Quashnock et al. (1989) Quashnock, J. M., Loeb, A., & Spergel, D. N. 1989, Astrophys. J. Lett., 344, L49
  * Quinlan (1996) Quinlan, G. D. 1996, New A, 1, 35
  * Rajagopal & Romani (1995) Rajagopal, M. & Romani, R. W. 1995, ApJ, 446, 543
  * Ramani et al. (2020) Ramani, H., Trickle, T., & Zurek, K. M. 2020, Journal of Cosmology and Astroparticle Physics, 2020, 033
  * Ravi et al. (2012) Ravi, V., Wyithe, J. S. B., Hobbs, G., et al. 2012, ApJ, 761, 84
  * Ravi et al. (2015) Ravi, V., Wyithe, J. S. B., Shannon, R. M., & Hobbs, G. 2015, MNRAS, 447, 2772
  * Ravi et al. (2014) Ravi, V., Wyithe, J. S. B., Shannon, R. M., Hobbs, G., & Manchester, R. N. 2014, MNRAS, 442, 56
  * Read (2014) Read, J. I. 2014, Journal of Physics G: Nuclear and Particle Physics, 41, 063101
  * Read et al. (2016) Read, J. I., Agertz, O., & Collins, M. L. M. 2016, Monthly Notices of the Royal Astronomical Society, 459, 2573
  * Reardon et al. (2023) Reardon, D. J., Zic, A., Shannon, R. M., et al. 2023, ApJ, 951, L6
  * Ringeval & Suyama (2017) Ringeval, C. & Suyama, T. 2017, JCAP, 12, 027
  * Roedig et al. (2011) Roedig, C., Dotti, M., Sesana, A., Cuadra, J., & Colpi, M. 2011, MNRAS, 415, 3033
  * Roedig & Sesana (2012) Roedig, C. & Sesana, A. 2012, in Journal of Physics Conference Series, Vol. 363, Journal of Physics Conference Series, 012035
  * Rogers & Peiris (2021) Rogers, K. K. & Peiris, H. V. 2021, Phys. Rev. Lett., 126, 071302
  * Roper Pol et al. (2020a) Roper Pol, A., Brandenburg, A., Kahniashvili, T., Kosowsky, A., & Mandal, S. 2020a, Geophys. Astrophys. Fluid Dynamics, 114, 130
  * Roper Pol et al. (2022a) Roper Pol, A., Caprini, C., Neronov, A., & Semikoz, D. 2022a, Phys. Rev. D, 105, 123502
  * Roper Pol et al. (2022b) Roper Pol, A., Mandal, S., Brandenburg, A., & Kahniashvili, T. 2022b, JCAP, 04, 019
  * Roper Pol et al. (2020b) Roper Pol, A., Mandal, S., Brandenburg, A., Kahniashvili, T., & Kosowsky, A. 2020b, Phys. Rev. D, 102, 083512
  * Rosado & Sesana (2014) Rosado, P. A. & Sesana, A. 2014, MNRAS, 439, 3986
  * Rosado et al. (2015) Rosado, P. A., Sesana, A., & Gair, J. 2015, MNRAS, 451, 2417
  * Rubakov et al. (1982) Rubakov, V. A., Sazhin, M. V., & Veryaskin, A. V. 1982, Physics Letters B, 115, 189
  * Rubin et al. (1980) Rubin, V. C., Ford, W. K., J., & Thonnard, N. 1980, ApJ, 238, 471
  * Rubin & Ford (1970) Rubin, V. C. & Ford, W. Kent, J. 1970, ApJ, 159, 379
  * Sachs & Wolfe (1967) Sachs, R. K. & Wolfe, A. M. 1967, ApJ, 147, 73
  * Saikawa & Shirai (2018) Saikawa, K. & Shirai, S. 2018, J. Cosmology Astropart. Phys., 2018, 035
  * Saito & Yokoyama (2009) Saito, R. & Yokoyama, J. 2009, Phys. Rev. Lett., 102, 161101
  * Sanidas et al. (2012) Sanidas, S. A., Battye, R. A., & Stappers, B. W. 2012, Phys. Rev. D, 85, 122003
  * Sasaki et al. (2018) Sasaki, M., Suyama, T., Tanaka, T., & Yokoyama, S. 2018, Classical and Quantum Gravity, 35, 063001
  * Schechter (1976) Schechter, P. 1976, ApJ, 203, 297
  * Schive et al. (2014) Schive, H.-Y., Liao, M.-H., Woo, T.-P., et al. 2014, Phys. Rev. Lett., 113, 261302
  * Schwarz & Stuke (2009) Schwarz, D. J. & Stuke, M. 2009, J. Cosmology Astropart. Phys., 2009, 025
  * Sesana (2010) Sesana, A. 2010, ApJ, 719, 851
  * Sesana (2013a) Sesana, A. 2013a, Classical and Quantum Gravity, 30, 224014
  * Sesana (2013b) Sesana, A. 2013b, MNRAS, 433, L1
  * Sesana (2015) Sesana, A. 2015, in Astrophysics and Space Science Proceedings, Vol. 40, Gravitational Wave Astrophysics, 147
  * Sesana et al. (2014) Sesana, A., Barausse, E., Dotti, M., & Rossi, E. M. 2014, ApJ, 794, 104
  * Sesana et al. (2004) Sesana, A., Haardt, F., Madau, P., & Volonteri, M. 2004, ApJ, 611, 623
  * Sesana et al. (2008) Sesana, A., Vecchio, A., & Colacino, C. N. 2008, MNRAS, 390, 192
  * Sesana et al. (2009) Sesana, A., Vecchio, A., & Volonteri, M. 2009, MNRAS, 394, 2255
  * Siemens et al. (2007) Siemens, X., Mandic, V., & Creighton, J. 2007, Phys. Rev. Lett., 98, 111101
  * Simon (2023) Simon, J. 2023, ApJ, 949, L24
  * Sivertsson et al. (2018) Sivertsson, S., Silverwood, H., Read, J. I., Bertone, G., & Steger, P. 2018, Monthly Notices of the Royal Astronomical Society, 478, 1677
  * Siwek et al. (2020) Siwek, M. S., Kelley, L. Z., & Hernquist, L. 2020, MNRAS, 498, 537
  * Smarra et al. (2023) Smarra, C., Goncharov, B., Barausse, E., & EPTA and InPTA. 2023, Phys. Rev. Lett., 131, 171001
  * Sorbo (2011) Sorbo, L. 2011, J. Cosmology Astropart. Phys., 2011, 003
  * Sorbo (2011) Sorbo, L. 2011, JCAP, 06, 003
  * Sotiriou & Faraoni (2010) Sotiriou, T. P. & Faraoni, V. 2010, Reviews of Modern Physics, 82, 451
  * Springel et al. (2005) Springel, V., White, S. D. M., Jenkins, A., et al. 2005, Nature, 435, 629
  * Starobinskii (1985) Starobinskii, A. A. 1985, Soviet Astronomy Letters, 11, 133
  * Svrcek & Witten (2006) Svrcek, P. & Witten, E. 2006, Journal of High Energy Physics, 2006, 051
  * Taylor & Gair (2013) Taylor, S. R. & Gair, J. R. 2013, Phys. Rev. D, 88, 084001
  * Taylor et al. (2020) Taylor, S. R., van Haasteren, R., & Sesana, A. 2020, Phys. Rev. D, 102, 084039
  * The Nanograv Collaboration et al. (2023) The Nanograv Collaboration, Agazie, G., Anumarlapudi, A., et al. 2023, ApJ, submitted
  * Tiburzi et al. (2016) Tiburzi, C., Hobbs, G., Kerr, M., et al. 2016, MNRAS, 455, 4339
  * Tomita (1967) Tomita, K. 1967, Progress of Theoretical Physics, 37, 831
  * Tristram et al. (2022) Tristram, M., Banday, A. J., Górski, K. M., et al. 2022, Phys. Rev. D, 105, 083524
  * Vachaspati & Vilenkin (1985) Vachaspati, T. & Vilenkin, A. 1985, Phys. Rev. D, 31, 3052
  * Vallisneri (2020) Vallisneri, M. 2020, libstempo: Python wrapper for Tempo2, Astrophysics Source Code Library, record ascl:2002.017
  * Valtolina et al. (2024) Valtolina, S., Shaifullah, G., Samajdar, A., & Sesana, A. 2024, A&A, 683, A201
  * Vaskonen & Veermäe (2021) Vaskonen, V. & Veermäe, H. 2021, Phys. Rev. Lett., 126, 051303
  * Verbiest et al. (2016) Verbiest, J. P. W., Lentati, L., Hobbs, G., et al. 2016, MNRAS, 458, 1267
  * Vilenkin & Shellard (2000) Vilenkin, A. & Shellard, E. P. S. 2000, Cosmic Strings and Other Topological Defects
  * Vovchenko et al. (2021) Vovchenko, V., Brandt, B. B., Cuteri, F., et al. 2021, Phys. Rev. Lett., 126, 012701
  * Witten (1984) Witten, E. 1984, Phys. Rev. D, 30, 272
  * Wygas et al. (2018) Wygas, M. M., Oldengott, I. M., Bödeker, D., & Schwarz, D. J. 2018, Phys. Rev. Lett., 121, 201302
  * Wyithe & Loeb (2003) Wyithe, J. S. B. & Loeb, A. 2003, ApJ, 590, 691
  * Xu et al. (2023) Xu, H., Chen, S., Guo, Y., et al. 2023, Research in Astronomy and Astrophysics, 23, 075024
  * Xue et al. (2021) Xue, X. et al. 2021, Phys. Rev. Lett., 127, 251303
  * Yi & Fei (2023) Yi, Z. & Fei, Q. 2023, European Physical Journal C, 83, 82
  * Zhang et al. (2019) Zhang, J., Liu, H., & Chu, M.-C. 2019, Frontiers in Astronomy and Space Sciences, 5
  * Zhao & Wang (2023) Zhao, Z.-C. & Wang, S. 2023, Universe, 9, 157

##  Appendix A Supermassive black hole binaries - full corner plots

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

Figure 24:  Marginalised posterior distributions for all $18$ parameters of
the astrophysically-informed model. The posterior and prior are shown in grey
and green, respectively.  Figure 25:  Marginalised posteriors for all five
parameters of the agnostic SMBHB model.