# Practical approaches to analyzing PTA data: Cosmic strings with six pulsars

Hippolyte Quelquejay Leclere, Pierre Auclair, Stanislav Babak, Aurélien Chalumeau, Danièle A Steer, J Antoniadis, A.-S. Bak Nielsen, C.G Bassa, A Berthereau, M Bonetti, et al.

► To cite this version:

Hippolyte Quelquejay Leclere, Pierre Auclair, Stanislav Babak, Aurélien Chalumeau, Danièle A Steer, et al.. Practical approaches to analyzing PTA data: Cosmic strings with six pulsars. Physical Review D, 2023, 108 (12), pp.123527. 10.1103/PhysRevD.108.123527. hal-04152091

HAL Id: hal-04152091

https://hal.science/hal-04152091v1

Submitted on 15 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Practical approaches to analyzing PTA data: Cosmic strings with six pulsars

Hippolyte Quelquejay Leclere,¹, * Pierre Auclair,²,† Stanislav Babak,¹, ‡ Aurélien Chalumeau,³

Danièle A. Steer,¹, § J. Antoniadis,⁴,⁵ A.-S. Bak Nielsen,⁵,⁶ C. G. Bassa,⁷ A. Berthereau,⁸,⁹

M. Bonetti,³,¹⁰,¹¹ E. Bortolas,³,¹⁰,¹¹ P. R. Brook,¹² M. Burgay,¹³ R. N. Caballero,¹⁴ D. J. Champion,⁵

S. Chanlaridis,⁴ S. Chen,¹⁵ I. Cognard,⁸,⁹ G. Desvignes,⁵ M. Falxa,⁸,¹ R. D. Ferdman,¹⁶

A. Franchini,³,¹⁰ J. R. Gair,¹⁷ B. Goncharov,¹⁸,¹⁹ E. Graikou,⁵ J.-M. Grießmeier,⁸,⁹ L. Guillemot,⁸,⁹

Y. J. Guo,⁵ H. Hu,⁵ F. Iraci,¹³,²⁰ D. Izquierdo-Villalba,³,¹⁰ J. Jang,⁵ J. Jawor,⁵ G. H. Janssen,⁷,²¹

A. Jessner,⁵ R. Karuppusamy,⁵ E. F. Keane,²² M. J. Keith,²³ M. Kramer,⁵,²³ M. A. Krishnakumar,⁵,⁶

K. Lackeos,⁵ K. J. Lee,⁴,⁵,⁹ K. Liu,⁵ Y. Liu,⁶,²⁴ A. G. Lyne,²³ J. W. McKee,²⁵,²⁶ R. A. Main,⁵

M. B. Mickaliger,²³ I. C. Nițu,²³ A. Parthasarathy,⁵ B. B. P. Perera,²⁷ D. Perrodin,¹³ A. Petiteau,²⁸,¹

N. K. Porayko,⁵,³ A. Possenti,¹³ A. Samajdar,²⁹ S. A. Sanidas,²³ A. Sesana,³,¹⁰,¹¹ G. Shaifullah,³,¹⁰,¹³

L. Speri,¹⁷ R. Spiewak,²³ B. W. Stappers,²³ S. C. Susarla,³⁰ G. Theureau,⁸,⁹,³¹ C. Tiburzi,¹³

E. van der Wateren,²¹,⁷ A. Vecchio,¹² V. Venkatraman Krishnan,⁵ J. Wang,⁶,³²,³³ L. Wang,²³ and Z. Wu²⁴

(The European Pulsar Timing Array)

¹Université Paris Cité CNRS, Astroparticule et Cosmologie, 75013 Paris, France

²Cosmology, Universe and Relativity at Louvain (CURL),
Institute of Mathematics and Physics, University of Louvain,
2 Chemin du Cyclotron, 1348 Louvain-la-Neuve, Belgium

³Dipartimento di Fisica “G. Occhialini”, Università degli Studi di Milano-Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy

⁴Institute of Astrophysics, FORTH, N. Plastira 100, 70013, Heraklion, Greece

⁵Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany

⁶Fakultät für Physik, Universität Bielefeld, Postfach 100131, 33501 Bielefeld, Germany

⁷ASTRON, Netherlands Institute for Radio Astronomy,
Oude Hoogeveensdijk 4, 7991 PD, Dwingeloo, The Netherlands

⁸Laboratoire de Physique et Chimie de l'Environnement et de l'Espace,
Université d'Orléans / CNRS, 45071 Orléans Cedex 02, France

⁹Observatoire Radioastronomique de Nançay, Observatoire de Paris,
Université PSL, Université d'Orléans, CNRS, 18330 Nançay, France

¹⁰INFN, Sezione di Milano-Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy

¹¹INAF - Osservatorio Astronomico di Brera, via Brera 20, I-20121 Milano, Italy

¹²Institute for Gravitational Wave Astronomy and School of Physics and Astronomy,
University of Birmingham, Edgbaston, Birmingham B15 2TT, UK

¹³INAF - Osservatorio Astronomico di Cagliari, via della Scienza 5, 09047 Selargius (CA), Italy

¹⁴Hellenic Open University, School of Science and Technology, 26335 Patras, Greece

¹⁵Kavli Institute for Astronomy and Astrophysics,
Peking University, Beijing 100871, P. R. China

¹⁶School of Physics, Faculty of Science, University of East Anglia, Norwich NR4 7TJ, UK

¹⁷Max Planck Institute for Gravitational Physics (Albert Einstein Institute), Am Muühlenberg 1, 14476 Potsdam, Germany

¹⁸Gran Sasso Science Institute (GSSI), I-67100 L'Aquila, Italy

¹⁹INFN, Laboratori Nazionali del Gran Sasso, I-67100 Assergi, Italy

²⁰Universitá di Cagliari, Dipartimento di Fisica,
S.P. Monserrato-Sestu Km 0,700 - 09042 Monserrato (CA), Italy

²¹Department of Astrophysics/IMAPP, Radboud University Nijmegen,
P.O. Box 9010, 6500 GL Nijmegen, The Netherlands

²²School of Physics, Trinity College Dublin, College Green, Dublin 2, D02 PN40, Ireland

²³Jodrell Bank Centre for Astrophysics, Department of Physics and Astronomy,
University of Manchester, Manchester M13 9PL, UK

²⁴National Astronomical Observatories, Chinese Academy of Sciences, Beijing 100101, P. R. China

²⁵E.A. Milne Centre for Astrophysics, University of Hull,
Cottingham Road, Kingston-upon-Hull, HU6 7RX, UK

²⁶Centre of Excellence for Data Science, Artificial Intelligence and Modelling (DAIM),
University of Hull, Cottingham Road, Kingston-upon-Hull, HU6 7RX, UK

²⁷Arecibo Observatory, HC3 Box 53995, Arecibo, PR 00612, USA

²⁸IRFU, CEA, Université Paris-Saclay, F-91191 Gif-sur-Yvette, France

²⁹Institut für Physik und Astronomie, Universität Potsdam,
Haus 28, Karl-Liebknecht-Str. 24/25, 14476, Potsdam, Germany

³⁰Ollscoil na Gaillimhe — University of Galway, University Road, Galway, H91 TK33, Ireland

³¹Laboratoire Univers et Théories LUTH, Observatoire de Paris,
Université PSL, CNRS, Université de Paris, 92190 Meudon, France
³²Ruhr University Bochum, Faculty of Physics and Astronomy, Astronomical Institute (AIRUB), 44780 Bochum, Germany

³³Advanced Institute of Natural Sciences, Beijing Normal University, Zhuhai 519087, China

(Dated: May 6, 2024)

We search for a stochastic gravitational wave background (SGWB) generated by a network of cosmic strings using six millisecond pulsars from Data Release 2 (DR2) of the European Pulsar Timing Array (EPTA). We perform a Bayesian analysis considering two models for the network of cosmic string loops, and compare it to a simple power-law model which is expected from the population of supermassive black hole binaries. Our main strong assumption is that the previously reported common red noise process is a SGWB. We find that the one-parameter cosmic string model is slightly favored over a power-law model thanks to its simplicity. If we assume a two-component stochastic signal in the data (supermassive black hole binary population and the signal from cosmic strings), we get a 95% upper limit on the string tension of $log_{10}(G\mu) < -9.9$ ($-10.5$) for the two cosmic string models we consider. In extended two-parameter string models, we were unable to constrain the number of kinks. We test two approximate and fast Bayesian data analysis methods against the most rigorous analysis and find consistent results. These two fast and efficient methods are applicable to all SGWBs, independent of their source, and will be crucial for analysis of extended data sets.

# I. INTRODUCTION

All regional Pulsar Timing Array (PTA) collaborations have recently, and independently, reported the presence of a common red noise process in the observed data [1–3]. Furthermore, the combined dataset shows even stronger evidence for its presence [4]. The exact nature of this signal is as yet very uncertain. While its spectral properties are consistent with an expected stochastic gravitational wave background (SGWB), the data is not sensitive enough to make any informative statement for or against its gravitational wave (GW) nature <fnref n="1" />.

In this paper, we make a strong assumption that it is, in fact, a SGWB and work towards its possible interpretation. The most favorable model for an anticipated SGWB is a superposition of monochromatic GW signals from a population of supermassive black hole binaries (SMBHBs) in the local Universe [1, 10]. This signal, assuming SMBHBs in circular orbits, is a power-law with the theoretical spectral index $\gamma = 13/3$ in the power spectral density of the residuals [11]. However, realistic astrophysical simulations suggest that the spectral index could vary based on the realization of the observed Universe [12]. For this reason, the spectral index is usually inferred from the data, giving a two-parameter model: A, $\gamma$. The amplitude A is referenced at the frequency 1/year.

We propose an alternative to SMBHBs, namely a SGWB from an early universe source, and in particular a

network of cosmic string loops. (See [13–17] for previous work on cosmic strings in the nHz band.) Cosmic strings are topological defects that could have emerged from symmetry-breaking phase transitions in the early Universe [18–21]. These quasi one-dimensional objects are characterized by their dimensionless tension $G\mu$, where G is Newton’s constant. Numerical simulations show that on large scales, cosmic string networks reach an attractor ‘scaling’ regime in which all the characteristic length scales in the network grow as $t$ [22–25]. Cosmic string loops are formed at all times by the scaling infinite string network, oscillate with period $\ell/2$, and decay into GWs<fnref n="2" />. The superposition of the GW emitted generates a SGWB that depends on the cosmic string loop distribution. While the loop distribution is well-understood on large scales close to the Hubble radius, there is an active debate regarding the distribution of smaller loops. To account for this theoretical uncertainty, in the following we consider two loop distribution models that have been used by the LIGO-Virgo-Kagra [30, 31] and LISA [32, 33] collaborations.

In this paper we use the six pulsars (see Table I) from the early Data Release 2 (DR2) of European PTA (EPTA) collaboration [1] <fnref n="3" /> Six pulsars, at the current sensitivity and observational time span, are not sufficient to detect the stochastic GW signal, however the presence of the common red process is evident. More pulsars are required to populate the expected Hellings-Downs (HD) [34] correlation curve and to get a high statistical significance for a SGWB. Nevertheless, in our study, we still utilize Hellings-Downs correlations to analyze the com-

<footnote>* quelquejay@apc.in2p3.fr</footnote>

<footnote>† pierre.auclair@uclouvain.be</footnote>

<footnote>‡ stas@apc.in2p3.fr</footnote>

<footnote>§ steer@apc.in2p3.fr</footnote>

<footnote n="1">¹ After submission of this paper on the ArXiv and during the peer review process, the results from the regional PTA collaborations became public [5–9]. These papers have demonstrated evidence of GW-induced correlations in their respective datasets with a considerable statistical significance.</footnote>

<footnote n="2">² Field theory simulations show that cosmic string loops also decay into particles [26, 27]. Although the balance between GW and particle emission is still under debate, preliminary quantitative studies tend to show that the SGWB from cosmic string is unaffected at the frequency of PTAs [28, 29]</footnote>

<footnote n="3">³ This dataset is a subset of the full DR2 dataset recently used in [5] which includes a total of 25 pulsars.</footnote>
mon red noise among pulsars. This approach allows us to account for the correlations among pulsars when characterizing the common red noise (CRN), resulting in a distinct interaction with the individual red noises. Given that we consider only 6 pulsars, the diagonal noise component (autocorrelation part of the noise matrix) will dominate and we expect only a small difference from the inclusion of Hellings-Downs correlation. Here we describe in detail the methodology of inferring parameters of the string network. This same methodology will later be applied on the extended dataset [5].

Assuming that the observed common red noise is a SGWB, we use Bayesian methods to infer the characteristics of the string network using a one-parameter model ($G\mu$) and a two-parameter model ($G\mu, N_k$) where $N_k$ is the average number of kinks on a loop per oscillation period, see section II. In addition to the SGWB signal, we also model the three components of the noise: white noise, pulsar red noise and chromatic noises (dispersion measurement variations and scattering variations) (see [35] for details).

We carry the Bayesian analysis through three different approaches with varying accuracy of description and computational cost. The first approach, dubbed in this paper as “Full method” [36, 37], is the standard approach in which the parameters of the SGWB and of the noise are explored simultaneously. This method is the most accurate but also the most computationally expensive due to the high-dimensionality of the analysis and because we consider Hellings-Downs correlations for the SGWB. The second approach, dubbed as “Resampling method” [38], neglects at first the spatial correlation between pulsars considering a common uncorrelated signal. The resulting posterior distributions are then resampled by taking into account GW-induced correlation. In this approach, the likelihood is factorized into a product of likelihoods for each pulsar [39], and the computational cost is significantly reduced with respect to the Full method as the number of pulsars increases. Finally, the third method, dubbed as “Free spectrum” [40, 41], amounts to obtaining the correlated power of the common noise for each frequency bin independently, marginalizing over the single pulsar noise parameters, before inferring the SGWB parameters. This method drastically reduces the dimensionality of the analysis and allows for a very fast parameter estimation of the SGWB models.

The paper is organized as follows. We start with a description of the SGWB produced by cosmic string loops in section II. In section III we explain the three data analysis approaches used. In section IV, we present the constraints on single parameter models ($G\mu$) and two parameter models ($G\mu, N_k$). We perform a Bayesian model comparison between the different cosmic string models and a simple power-law SGWB. Finally, we conclude in section V.

## II. STOCHASTIC GRAVITATIONAL WAVE BACKGROUND FROM COSMIC STRINGS

The SGWB generated by a cosmic string network has been studied in depth, see e.g. [32, 42]. It includes a contribution from the uncorrelated superposition GW bursts emitted from cusps, kinks and kink-kink collisions on oscillating loops [43]. Cusps, which are points at which the loop reaches ultra relativistic velocities, emit a short beam of GWs.

Kinks, which are discontinuities in the tangent vector of the loop, propagate along the string at the speed of light and emit a one-dimensional “fan-like” GW burst with the same beaming angle as for cusps [43]. Finally, when left- and right-moving kinks collide (kink-kink collisions), an isotropic burst of GW emission is emitted [44].

A loop of length $\ell$ oscillates periodically with corresponding fundamental frequency $f_0 = 2/\ell$. For the n-th oscillatory mode of the loop, the power (in units of $G\mu^2$) emitted in GWs by each of these burst events (labelled by b) is given by 4 [43]

$$P_n^{(b)} = \frac{\Gamma^{(b)}}{\zeta(q_b)} n^{-q_b}. \quad (1)$$

Here $\zeta(q_b)$ is the Riemann zeta function, $q_b = 4/3, 5/3$ and 2 for cusps, kinks and kink-kink collisions respectively and [43, 44]

$$\Gamma^c = \frac{3(\pi g_1^c)^2}{2^{1/3} g_2^{2/3}}, \quad \Gamma^k = \frac{3(\pi g_1^k)^2}{2^{2/3} g_2^{1/3}}, \quad \Gamma^{kk} = 2(\pi g_1^{kk})^2, \quad (2)$$

where $g_2 = \sqrt{3}/4, g_1^c \approx 0.85, g_1^k \approx 0.29, g_1^{kk} \approx 0.1$.

Given the exponents $q_b$ in Eq. (1), one might expect that, compared to the cusp contribution, the kink and kink-kink collision contributions could be neglected. However, as already noted in [44], if many kinks are present they can dominate the power emitted in GWs. Indeed, since two kinks (right and left moving) are formed when two pieces of strings intercommute, the number of kink-kink collisions on a loop with $N_k$ kinks per oscillation period should scale as $N_k^2/4$, which can then dominate for high $N_k$. Furthermore, kinks can proliferate on loops with junctions [45, 46]. To summarize, for a loop containing $N_c$ cusps and $N_k$ kinks per oscillation period, the total power emitted, in units of $G\mu^2$, is

$$\Gamma = N_c \Gamma^c + N_k \Gamma^k + \frac{N_k^2}{4} \Gamma^{kk}. \quad (3)$$

The fractional energy density of the SGWB per logarithmic interval of frequency is

$$\Omega_{gw}(t_0, f) \equiv \frac{8\pi G}{3H_0^2} f \frac{d\rho_{gw}}{df}, \quad (4)$$

<footnote n="4">4 Formally this expression is only valid for large $n$, but as elsewhere we assume that extrapolating to $n \ge 1$ is a good approximation to describe low-harmonic modes.</footnote>
where $H_0$ is the Hubble constant, and $d\rho_{gw}/df$ is the energy density in gravitational waves per unit frequency $f$, observed today (at $t = t_0$). Following [31, 47], for a network of cosmic string loops with distribution $n(\ell, t)$, this is given by

$$ \frac{d\rho_{gw}}{d\ln f} = \frac{2G\mu^2}{f} \sum_b \frac{N_b \Gamma^{(b)}}{\zeta(q_b)} \times \sum_{n=1}^{+\infty} \int \frac{n^{1-q_b} dz}{(1+z)^6 H(z)} \mathbf{n} \left[ \frac{2n}{(1+z)f}, t(z) \right], \quad (5) $$

where we have summed over the cusp, kink and kink-kink contributions, and $H(z)$ is the Hubble parameter. In the following we consider standard $\Lambda$CDM cosmology with the Planck-2018 fiducial parameters [48].

We consider the most updated loop distribution models $n(\ell, t)$ available, calibrated with numerical simulations, that have also been used in LVK collaboration papers [30, 31] and LISA [33, 49]. They are denoted in the following as BOS [42, 50] and LRS [51]. Each has an intersection probability $p = 1$, and the largest loops produced are of size $0.1t$. The main difference between these models is that, relative to the BOS model, the LRS model has an additional population of very small loops which aim to account for physics at the gravitational backreaction scale. These small loops emit high frequency GWs, and lead to modifications of the SGWB at high frequencies $f \gg H_0(\Gamma G\mu)^{-1}$ [52]. This can be seen in Figure 1 which plots the numerical evaluation of Eq. (5) for $N_c = 2$, $N_k = 0$ for the BOS model (continuous lines) and LRS model (dashed lines) for different values of $G\mu$.

Notice that as a consequence, the two models considered give a significantly different spectrum in PTA frequency range only for $\Gamma G\mu \gtrsim 2 \times 10^{-11}$ (as seen for $G\mu = 10^{-14}$ and $\Gamma = 57$ in Figure 1).

In this work, we consider two cases (i) smooth cosmic string loops with two cusps only (as in Figure 1), (ii) fixed number of cusps, $N_c = 1$, but a varying number of kinks $N_k$ with $N_k^2/4$ kink-kink collisions. Our aim is to quantify how well the SGWB from cosmic string loops can explain the common red process that is already seen with strong evidence in three PTAs, namely the EPTA, NANOGRAV and PPTA consortiums [1–3].

### III. DATA ANALYSIS METHOD

In this section, we use a Bayesian analysis to perform parameter estimation and model selection. We will not repeat the basis of Bayesian inference in the context of PTA, and instead refer the reader to [1, 35] for a more detailed description. In the following subsection, we briefly describe the noise and SGWB models which we test.

We have applied three Bayesian approaches to the six pulsars of the early second EPTA data release (DR2) [1]:

1. We directly infer the parameters of each model using the spectral shape of the SGWB as a prior,

<figure><img src="image_2.png" /><figcaption>FIG. 1. SGWB from a network of cosmic string loops, expressed in terms of characteristic energy density. The spectrum is computed using the BOS (resp. LRS) loop number density model for the solid (resp. dashed) lines. Here we have taken $N_c = 2$, $N_k = 0$ leading to $\Gamma = 57$. For each model, computations using three different tension values $G\mu$ are represented. The sensitivity frequency range of EPTA corresponds to the yellow band.</figcaption></figure>

considering HD correlation between pulsars. This is the most computationally expensive path, which is referred to as “full”.

2. Similar to 1, but initially we neglect all correlations between pulsars, assuming an uncorrelated common process. Then, we resample the posterior points using “correct” likelihood with HD correlation following the procedure outlined in [38]. This path is significantly faster, we refer to it as “resampling” (RS).

3. Finally, we also use the method suggested in [41] and further developed in [40]. There, we first estimate the correlated power at each Fourier bin with a relatively large log-uniform prior: we refer to this as the “free spectrum”. Then, we perform Bayesian analysis using the free spectrum as observational data. This is the computationally cheapest path, which we refer to as the “free spectrum” (FS) method.

#### A. PTA data analysis principles

In PTA, we work with time residuals, obtained as the difference between observed and predicted time of arrivals of the pulsars’ radio pulses

$$ \delta t_i = t_i^{\text{obs}} - t_i^{\text{TM}}(\vec{\beta}_b), \quad (6) $$

where $i$ is the index of each observation and $t^{\text{TM}}$ is the timing model attempting to explain the deviations from
the observed time-of-arrival using a set of parameters $\vec{\beta}_b$. Following the usual path [36, 37] we marginalize analytically over timing model parameters.

We assume that the noise in each pulsar consists of the radiometer white noise, intrinsic pulsar spin noise, dispersion measurement and finally scattering variations. These last two components are referred to as chromatic noises, since they are due to the propagation of the radio pulses through the interstellar medium and depend on the observational radio frequency. We follow the standard procedure [35], initially we infer the noise model in each pulsar individually, and then we use these results as a starting point in the search for the SGWB common to all pulsars.

Besides these individual noise components, we also introduce SGWB into the data model. It is characterized by the spectral shape (subject to its source) and exhibits a particular correlation across pulsar pairs which depends on the angular separation of pulsars in the sky and is described by the Hellings-Downs (HD) curve [34].

We assume that the data does not contain any deterministic signals, with the exception of the two deterministic chromatic signals, referred to as exponential dips, of PSR J1713+0747 (see Table I).

We can summarize the data model in the form of the total covariance matrix

$$C_{(ai)(bj)} = \mathcal{N}_{a,(ij)}\delta_{ab} + (C_{a,(ij)}^{\text{RN}} + C_{a,(ij)}^{\text{ChRN}})\delta_{ab} + \Gamma_{ab}C_{(ij)}^{\text{CRN}}, \quad (7)$$

where $\delta_{(..)}$ is the Kronecker delta function, $a, b$ are indexing the pulsars, $i, j$ the timing residuals, $\mathcal{N}_a$ is the white noise covariance matrix of pulsar $a$, $C_a^{\text{RN}}$ its intrinsic red-noise covariance matrix and $C_a^{\text{ChRN}}$ its chromatic red-noise covariance matrix. $C_a^{\text{CRN}}$ is the common red noise covariance matrix, and $\Gamma_{ab}$ is the overlap reduction function giving the correlations between pulsars $a$ and $b$. The function $\Gamma_{ab}$ reduces to $\delta_{ab}$ for an uncorrelated common red noise and is described by the HD correlations for a SGWB,

$$\Gamma_{ab} = \frac{3}{2}x \ln x - \frac{1}{4}x + \frac{1}{2}(1 + \delta_{ab}), \quad (8)$$

where $x = \frac{1-\cos\xi_{ab}}{2}$, $\xi_{ab}$ being the angular separation between pulsars $a$ and $b$.

In PTA data analysis, red stochastic processes are often described via Gaussian Processes with truncated Fourier basis (see [37] for details)

$$\delta t_{\text{red}}(t_i) = \sum_{n=1}^{N_f} \left[ a_n \cos \left( \frac{2\pi nt_i}{T} \right) + b_n \sin \left( \frac{2\pi nt_i}{T} \right) \right], \quad (9)$$

where $T$ is the observational time span and $a_n, b_n$ are random variables (weights) defined by the spectral shape of the underlying process. In this work, we use $N_f = 30$ for the common red process while for the individual pulsar noises, we employ a custom number of Fourier bins determined through single pulsar analysis, following the

method detailed in [35] (see Table I for an overview of the different $N_f$ values used for each pulsar noise). For the chromatic noises, we use a similar approach but introduce a scaling with the observational radio frequency $\nu$, $\delta t_{\text{red}} \propto \nu^{-\chi}$, where $\chi$ is a chromatic index equals to 2 for dispersion measure variations and 4 for scattering variations.

For all intrinsic red noise components attached to each pulsar, we assume a power-law shape for the one-sided power spectral density (PSD) of the residuals

$$S(f; A, \gamma) = \frac{A^2}{12\pi^2} \left( \frac{f}{\text{yr}^{-1}} \right)^{-\gamma} \text{yr}^3, \quad (10)$$

characterized by the spectral index $\gamma$ and the amplitude $A$ defined at the reference frequency of 1/year.

Here we assume that the CRN is a SGWB with the one-sided power spectral density given in terms of the fractional energy density by

$$S(f; \vec{\eta}) = \frac{H_0^2}{8\pi^4} \frac{1}{f^5} \Omega_{\text{SGWB}}(f; \vec{\eta}), \quad (11)$$

where $\vec{\eta}$ are the hyper-parameters of the SGWB model considered. The SGWB from a population of SMBHB is well approximated by a power-law model, and takes the same form as in (10), see [10]. Moreover, for the SMBHBs in circular orbits with GW-driven evolution, we expect $\gamma \approx 13/3$. The spectral shape for the network of cosmic strings is obtained by integration of Eq. (5) which adds computational time.

In the following analysis, we fix the parameters (and, therefore the level) of the white noise based on each pulsar investigation [35]. We vary the parameters of all red noise components together with the parameters of SGWB. In total, there are $30 + N_{\text{CRN}}$ parameters, where $N_{\text{CRN}}$ is the number of hyper-parameters of the CRN (varies from model to model).

## B. “Full” Method

This is the standard path one would follow to infer the parameters of each model (pulsars and SGWB). The likelihood for the concatenated array of observations (residuals) $\vec{\delta}t$ is given by the usual Gaussian form

$$p(\vec{\delta}t|\vec{\eta}) = \frac{1}{\sqrt{|2\pi C|}} \exp\left(-\frac{1}{2}\vec{\delta}t^T C^{-1} \vec{\delta}t\right), \quad (12)$$

where $C$ is the full covariance matrix of Eq. (7). It is parametrized by $A_a, \gamma_a$ for each red noise component for each pulsar $a$ and by parameters describing the SGWB spectrum. We use the python package ENTERPRISE [53] for its computation (with an extension to include cosmic string models). We sample parameters using PTMCMC sampler [54]. Note that we do not employ the “parallel tempering” features due to technical reasons.
<table><thead><tr><th>Pulsar</th><th>Time Span (years)</th><th>Number of ToAs</th><th>Red Noise</th><th>DM Variations</th><th>Scattering Variations</th><th>Deterministic Signals</th></tr></thead><tbody><tr><td>J0613-0200</td><td>22.89</td><td>2909</td><td>10</td><td>144</td><td>—</td><td>—</td></tr><tr><td>J1012+5307</td><td>23.68</td><td>5325</td><td>149</td><td>45</td><td>—</td><td>—</td></tr><tr><td>J1600-3053</td><td>14.32</td><td>2982</td><td>—</td><td>26</td><td>137</td><td>—</td></tr><tr><td>J1713+0747</td><td>24.46</td><td>5003</td><td>11</td><td>148</td><td>—</td><td>Two exponential dips<sup>a</sup></td></tr><tr><td>J1744-1134</td><td>24.01</td><td>1946</td><td>9</td><td>151</td><td>—</td><td>—</td></tr><tr><td>J1909-3744</td><td>15.74</td><td>2503</td><td>20</td><td>151</td><td>—</td><td>—</td></tr></tbody></table>

a Interested readers are directed to section 5.2.2 of [35] for further details regarding these events.

<caption>TABLE I. Noise model used for each of six pulsars. For each pulsar, we also consider white noise whose parameters are set to the maximum likelihood values obtained through single pulsar analysis. We indicate the number of Fourier bins ($N_f$ in Equation 9) employed for each noise type and pulsar. In cases where a particular noise type is not included in the pulsar's noise model, we denote it with a – sign.</caption>

For the SGWB we have used either the power-law model, the BOS or the LRS models with fixed and variable number of kinks. This is computationally the most expensive (but most rigorous) approach. There are two bottlenecks in the analysis: (i) integration of Eq. (5), and (ii) accounting for the correlation between pulsars (inverting the $\Gamma_{ab}$ part in the covariance matrix). However, thanks to the low dimensionality of the parameter space of the cosmic string spectrum ($f, G\mu, N_k$), we were able to address point (i) by using numerical interpolation method over a precomputed grid to speed-up the computation of its power spectrum.

We can also perform Bayesian model selection. We use the Bayes Factor (BF) as a way to measure which model is preferred by the data. The Bayes Factor is equivalent to the posterior odds ratio $\mathcal{O}$ if we assume (and we do) equal model priors

$$B_{M_2}^{M_1} = \frac{p(d|M_1)}{p(d|M_2)} = \frac{p(M_1|d)p(M_2)}{p(M_2|d)p(M_1)} = O_{M_2}^{M_1} \frac{p(M_2)}{p(M_1)}. \quad (13)$$

It is computationally prohibitive to compute the evidence for each model directly using the Nested sampling method. Instead, we used the product-space sampling approach described in [55] and already used in several PTA analysis (e.g. in [1, 56]). The artificially introduced hyperparameter (model index) allows the sampler also to jump between several models and the BF is given by a ratio of the number of samples accumulated in each model.

### C. Resampling method

This approach is based on importance sampling and described in detail in [38]. The main idea is to perform sampling like in the "full" method neglecting the correlation between pulsars, taking $\Gamma_{ab} = \delta_{ab}$ in (7). This way the likelihood is factorized into a product of likelihoods for each pulsar [39]. The likelihood computational time is reduced by a factor 6 for six pulsars, but this gain

increases as the number of pulsars analyzed increases. Once it is done, we resample (reweigh) the posteriors by recomputing the likelihood for each sample using $\Gamma_{ab}$ defined by the HD curve. As a result, it is possible to parallelize the calculation and therefore gain even more speed. This method works well if the difference in the posteriors obtained with and without HD correlation is not very large (not disjoint distributions).

We can use the same approach to get a cheap estimation of the Bayes factors. To the best of our knowledge, it was not used before, so we provide here some detailed computation.

We refer the reader to section 2 of [38] for details on the resampling of the posterior from an approximate model "A" to a target posterior "T". We apply this method to the product space approach for computing the Bayes Factor. We consider two hyper-models $\mathcal{M}, \mathcal{M}^{\text{HD}}$ with the same model parameters $\theta$ and hyperparameter $n$ indexing the models describing the common red noise we want to compare. The difference between the two hyper-models is in the likelihood computation: we consider uncorrelated common process in the hyper-model $\mathcal{M}$. Computations of Bayes factor with $\mathcal{M}$ are significantly faster and more reliable, we want to evaluate the Bayes factor of models in $\mathcal{M}^{\text{HD}}$ by resampling the hyper-chain (product space) $\mathcal{M}$. We repeat the steps outlined in section 3 of [55] but introducing the approximate likelihood in the posterior probability of $n$ for $\mathcal{M}^{\text{HD}}$

$$ \begin{aligned} P(n|\mathcal{D}, \mathcal{M}^{\text{HD}}) &= \frac{\pi(n)}{Z_{\mathcal{M}^{\text{HD}}}} \int \mathcal{L}^{\text{HD}}(\theta_n) \pi(\theta_n|n) d\theta_n \\ &= \frac{\pi(n)}{Z_{\mathcal{M}^{\text{HD}}}} \int \frac{\mathcal{L}^{\text{HD}}(\theta_n)}{\mathcal{L}(\theta_n)} \mathcal{L}(\theta_n) \pi(\theta_n|n) d\theta_n \\ &\approx \frac{\pi(n)}{Z_{\mathcal{M}^{\text{HD}}}} \frac{1}{N_s^{(n)}} \sum_{i=1}^{N_s^{(n)}} w_i^{(n)} Z_n. \end{aligned} \quad (14) $$

Here we introduced $\pi(n)$ the prior on the model index (chosen to be uniform), $\mathcal{L}^{\text{HD}}/\mathcal{L}$ are likelihood functions with/without HD spatial correlations and $Z_{\mathcal{M}^{\text{HD}}}$ is the
evidence of the hyper-model with HD (in our case just normalization factor which cancels out). In the last line, we have approximated the integral using Monte Carlo approach, with $N_s^{(n)}$ being the number of samples corresponding to the model associated to $n$ in the hyper-model $\mathcal{M}$ posterior. The weights $w_i^{(n)} = \mathcal{L}^{\text{HD}}(\theta_n^{(i)})/\mathcal{L}(\theta_n^{(i)})$ are computed at each posterior points (indexed by $i$) of $\mathcal{M}$, the $\theta_n$ being the $\theta$ parameters attached to the model $n$.

The Bayes Factor (assuming uniform prior across the models $\pi(n)$) between the models within each hyper-model posterior is given as

$$B_{M_b}^{M_a} \equiv B_{a/b} = \frac{P(\lfloor n \rfloor = n_a | \mathcal{D}, \mathcal{M})}{P(\lfloor n \rfloor = n_b | \mathcal{D}, \mathcal{M})}, \quad (15)$$

where $\lfloor \cdot \rfloor$ is the floor function.

Applying this expression to the resampled hyper-model posterior, we obtain

$$B_{a/b}^{\mathcal{M}} = \frac{Z_a^{\mathcal{H}\mathcal{D}}}{Z_b^{\mathcal{H}\mathcal{D}}} = \frac{\bar{w}_{(a)}}{\bar{w}_{(b)}} \times B_{a/b}^{\mathcal{M}}, \quad (16)$$

where $\bar{w}_{(n)} = \frac{1}{N_s^{(n)}} \sum_{i=1}^{N_s^{(n)}} w_i^{(n)}$ is a posterior average weight of the model indexed by $n$. As previously mentioned, the Bayes factor for the hyper-model $\mathcal{M}$ (without HD correlations) is simply given by $B_{a/b}^{\mathcal{M}} \approx N_s^{(n_a)}/N_s^{(n_b)}$. The performance of this method depends on the sampling efficiency, as described in [38]. Note that the novelty of this method is in the efficient evaluation of Bayes factor between two or more GWB models via resampling.

### D. “Free spectrum” method

As mentioned, this method was suggested in [41] and further developed in [40]. The main idea is to get rid of the high dimensionality of the problem by marginalizing over the single pulsar noise parameters. We want to estimate the PSD assuming the HD correlation at each frequency bin $f_k = k/T$, where $T$ is the total time span of the array. Using MCMC we get a posterior distribution for the amplitude $\rho_k$ at each frequency bin, where

$$\rho_k^2 = \frac{S(f_k)}{T}, \quad (17)$$

where $S(f)$ is the one-sided PSD of the SGWB we want to characterize, see Eq. (11). We assume that all bins are independent and impose log-uniform prior for each amplitude, $\log_{10} \rho_k \in [-12, -4]$, a justification for this range is given in section IV B.

We use a kernel density estimator to obtain a smooth probability distribution function $p(\rho_k|\vec{\delta}t)$. In this study, we use Gaussian kernels, and their bandwidth is selected using Scott's rule of thumb [57]. We use the estimated PSD of a common process with HD correlation as our new data. We can then build a factorized likelihood that

depends only on the parameters of the SGWB $\vec{\eta}$, following [40], we get

$$\begin{align} p(\vec{\delta}t|\vec{\eta}) &\approx \prod_{k=1}^{N_f} \int d\rho_k p(\vec{\delta}t|\rho_k) p(\rho_k|\vec{\eta}) \nonumber \\ &\propto \prod_{k=1}^{N_f} \int d\rho_k \frac{p(\rho_k|\vec{\delta}t)}{p(\rho_k)} p(\rho_k|\vec{\eta}). \tag{18} \end{align}$$

We have used the following notation: $p(\rho_k|\vec{\delta}t)$ are the posterior distributions of the $\rho_k$ describing our observations in Fourier domain, $p(\rho_k)$ is the prior on the $\rho_k$ and $p(\rho_k|\vec{\eta})$ is the probability to have $\rho_k$ given the CRN parameters $\vec{\eta}$. For the latter, we consider that the PSD is perfectly deterministic such that

$$p(\rho_k|\vec{\eta}) = \delta(\rho_k - \rho_k^{\text{SGWB}}(\vec{\eta})), \quad (19)$$

where $\rho_k^{\text{SGWB}}(\vec{\eta})$ is derived from the PSD of a corresponding SGWB model computed at $f_k = k/T$ for a set of background parameters $\vec{\eta}$.

Note that the Fourier frequency bins are not entirely independent as the PTA data are not evenly sampled and pulsars do not have the same observational duration. We have computed the associated Pearson correlation coefficient matrix $P_{i,j} = C_{i,j}/\sqrt{C_{i,i}C_{j,j}}$ where $C_{i,j}$ is the covariance matrix of the $\rho$'s. We have found that the frequency bins are approximately independent: the average of the non-diagonal coefficients being $\langle |P_{i,j}| \rangle_{i<j} \approx 0.02$.

Taking all into account, the likelihood of Eq. (18) takes the form

$$p(\vec{\delta}t|\vec{\eta}) \propto \prod_{k=1}^{N_c} p(\rho_k^{\text{SGWB}}(\vec{\eta})|\vec{\delta}t). \quad (20)$$

This method is extremely fast and drastically reduces the dimensionality of the problem, making it possible to use the Nested sampling [58, 59] algorithm to obtain the evidence for each model and posterior distributions of the SGWB parameters $\vec{\eta}$. We used a particular implementation of nested sampling, dynasty, described in [60, 61].

## IV. RESULTS

Before discussing our results, we should mention that the six pulsars dataset considered here is non-informative about the presence of the HD correlations, as was shown in [1, 35]. Consequently, we expect that the results of the analysis will be similar (though not identical) using CURN or GWB assumptions for the CRN. The main purpose of this paper is to introduce the methodology which will be applied to the extended (25-pulsar) EPTA dataset.

We start by considering the cosmic string models with the fixed average number of cosmic string bursts per oscillation to $N_c = 2$ and $N_k = N_{kk} = 0$, leading to $\Gamma = 57$, which is expected by a population of smooth loops [62]. Later we extend our analysis by allowing the number of kinks to vary (and so the value of $\Gamma$).
# A. Parameter estimation

## 1. Smooth loops

When we infer a SGWB from a network of smooth loops with $N_k = 0$ and $N_c = 2$ (alongside the individual pulsar noise models), we find a very constrained distribution for the tension $G\mu$ for both loop density distribution models, see Figure 2. As one can see the posterior distributions given by the Full, the RS and the FS methods agree very well, though we find that the FS method gives slightly broader posteriors. The string tension 90% credible (symmetric) intervals for each Bayesian method can be found in the two first lines of Table II. These intervals show nice consistency across all three Bayesian methods: $\log_{10} G\mu \sim -10.1$ (resp. $-10.6$) for the BOS (resp. LRS) model.

These results can be understood as follows. For the BOS model, the SGWB exhibits a peak around $f \sim F_0 \equiv 3H_0(\Gamma G\mu)^{-1}$ before decreasing to reach a plateau at higher frequencies (see Figure 1) [52]. According to the posterior for $G\mu$, the peak is in the middle of the PTA band, thus the spectral index of $\Omega_{\text{SGWB}}$ transits from 1 to 0 (flat spectrum) throughout the PTA frequency range. This is comparable to the value of $2/3$ (corresponding to $\gamma = 13/3$) expected for a SMBHB background. As can be seen in Figure 5, the BOS model slope is very compatible with the spectral estimated from the data.

Regarding the LRS model, the extra population of small loops (see [52]) is responsible for the fact that instead of decreasing to a plateau for $f > F_0$, the characteristic energy density rather increases again up to a second peak at very high frequencies compared to the PTA band, see Figure 1. The fractional energy density follows there a power-law of spectral index $\approx 0.4$, which also fits well the spectral shape of the data (see Figure 5).

As an alternative, we also consider a two-component SGWB model: a SGWB generated by a population of circular binary black holes with $\gamma = 13/3$ plus a stochastic signal from the cosmic strings.

We find a very strong correlation between those two components as one can see in Figure 3. The 2D-posteriors exhibit a pronounced L-shape, meaning that both model can explain the data. Based on the obtained posterior distribution, we can extract an upper limit which we evaluate to be $\log_{10} G\mu < -9.9$ (resp. $-10.5$) for the BOS (resp. LRS) model using the Full method, and it corresponds to 95% confidence. In addition, the evaluation obtained with the two fast methods (RS, [FS]) give very consistent upper bounds: $\log_{10} G\mu < -9.8[-9.5]$ (resp. $-10.5[-10.4]$) for the BOS (resp. LRS model). The slightly broader posterior obtained with FS earlier leads to a less stringent upper bound for the BOS model. One can notice that we constrain better the tension, $G\mu$, for the LRS model, this is due to the excess of power that the small loop population introduces.

## 2. Kinky loops

Next we also vary the number of kinks (and thus kink-kink collisions) in addition to the string tension, considering two-dimensional models ($\log_{10} G\mu, N_k$). The PSD of the SGWB created by such a model is calculated from Eq. (5) by setting $N_c = 1$ (in order to still have some GW power at the low number of kinks) and $N_{kk} = N_k^2/4$.

We took a uniform prior for $N_k$ letting it vary between 0 and 200, the same prior was used in the LVK analysis [31]. This large prior on $N_k$ accounts for theoretical uncertainties on the initial number of kinks at loop creation and on the efficiency of the gravitational backreaction that is expected to smooth out the loop.

In the case of the BOS model, if we increase the number of kinks (and so $\Gamma$), the SGWB in the PTA band now corresponds to the transition between the peak of the spectrum and the high frequency plateau and thus increases the characteristic slope of its PSD while decreasing the amplitude. This increase in spectral index will be disfavored as it will be too steep to fit the data properly. However, to compensate the decrease in amplitude, an increase of the string tension $G\mu$ can both correct the amplitude and the spectral index by placing the high frequency plateau in the PTA frequency band. One can see such interplay between $G\mu$ and $N_k$ in the left panel of Figure 4. The data equally allows two joined solutions: low tension with low number of kinks and very kinky loops ($N_k \gtrsim 120$) with high tension, $\log_{10} G\mu \sim -8.3$.

Contrary to the BOS case, the PSD slope of the LRS increases only slightly at low frequency with increasing $N_k$. Once $F_0 \lesssim 1/T \sim 10^{-9}$, the spectrum is dominated by the extra population of small loops that the LRS model introduces and the associated power is now independent of $\Gamma$ and proportional to $G\mu$ [32]. Therefore, the only shift in amplitude is caused by the domination of kinks (instead of cusp) in the GW power emission leading to an increase by a factor $\sim 1.6$ [32], that can be simply compensated by a slight decrease in the tension leading to $\log_{10} G\mu \sim 10.9$. As a result, adding kinks has practically no effect in the LRS model, and we recover the prior for $N_k$, as seen the right panel of Figure 4.

According to the posterior, the model with a high number of kinks is supported by the observed data, especially for BOS where it suggests a possible very high number of kinks. For example, $N_k = 120$ corresponds to $\Gamma \approx 960$, which is well above the value of 50 expected in the latest simulations of cosmic string loops. Finally, all three methods (depicted with different line styles in Figure 4) show very consistent results, even if the FS method struggles more to recover the double peak solution in the BOS model case.

# B. Model comparison

The fact that the cosmic string parameters are well constrained does not tell us how well it explains the ob-
<figure><img src="image_3.png" /><figcaption>FIG. 2. Comparison of the string tension posteriors (for two string models, BOS and LRS) obtained with the Full method (dashed lines), Resampling (RS) method (dotted lines) and with the free spectrum (FS) method (solid lines). We assume here that the loops are populated by two cusps, leading to $\Gamma = 57$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIG. 3. Left panel: posterior for the two-component SGWB model composed of (i) a signal originating from a population of circular GW-driven SMBHB, parameterised by its PSD amplitude $\log_{10} A$ at $f = 1$/year, and (ii) a SGWB from smooth CS loops background using the BOS loop number density model; different lines styles (dashed, dotted, solid) corresponding to three methods (Full, RS, FS) show good consistency. The 95-quantile for each SGWB parameter posterior is plotted in black, using the line style associated with its respective method. Right panel: the same for the LRS loop number density model.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>FIG. 4. Left panel: posterior for two-dimensional (string tension and the average number of kinks on loops of the network) BOS model; different lines styles corresponding to three methods (Full, RS, FS) show full consistency. Right panel: the same for the LRS loop number density model.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIG. 5. Posterior distributions of the 30 $\rho_k$ coefficients (in log$_{10}$-scale). We over-plotted the best fit (using the FS likelihood of Eq. (20)) for three different PSDs: powerlaw, SGWB using BOS and LRS models (in the case of smoothed loops). We see that all three spectra behave in a similar way at low frequency bins.</figcaption></figure>

served data. The Bayes factor, defined as the ratio of evidence between two models, is often used to quantify the ability of two different models to fit a set of data. Note that the Bayes factor can tell us which model is preferred based on the observations, but still does not tell us how well it describes the data. In what follows, we want to compute the Bayes factors of the considered cosmic string models against the power-law (PL) (SMBHB-inspired) model. Results are reported in Table II. We find that the cosmic string models are slightly preferred over PL model thanks to its simplicity: only one param-

eter. However, Figure 5 shows that all models fit the spectral shape similarly well.

Considering the two-component model of SGWB (SMBHB+CS, see previous subsection), we have found no statistical support for this complex model as compared to a simple SMBHB, $\mathcal{B}_{\text{SMBHB}}^{\text{SMBHB+CS}} \approx 1$ for both BOS and LRS cosmic string networks (using Full, RS and FS methods).

We have demonstrated earlier that the kinks had no strong impact and were not constrained by the data in the models with the variable number of kinks. This is
<figure><img src="image_7.png" /><figcaption>FIG. 6. Marginalized red noise parameters posteriors for the pulsar J1909-3744 obtained with the Full method including HD correlation for the common red noise process. Brown dashed and solid lines were obtained using different prior on the power in the Fourier bins (passing from -10 to -15 for the log₁₀ ρₖ lower bound). The solid line (restricted prior) suggests a truncated correlation of the common red noise with the spin red noise of J1909-3744. We also plotted for comparison those posteriors when using as common red noise, a SGWB from cosmic string following the BOS (blue line)/LRS (orange line) models.</figcaption></figure>

also reflected in the evaluated Bayes Factors.

Both Full and the RS methods give very consistent Bayes factors for all models, as shown in Table II.

The Bayes factors computed with FS method depend on the prior chosen for evaluating the FS amplitude components ρₖ. The main reason is that power at all low frequency bins is not fully constrained. The posteriors have a long tail towards the low values of power (see Figure 5). Choosing a wider prior makes those tails longer, spreading the probability across a long range of ρₖ (draining it from the "bumps" located at high power). Despite the fact that the tails are thin, they go all the way to minus infinity in log-scale (zero in power) and the probability distribution is then strongly affected (through normalization). So we are facing the question of what we should use as a prior range for the log₁₀ ρₖ?

The choice of the prior range has a direct impact on the evidence calculation: some values of the SGWB model parameters have no support (zero contribution to the evidence integral) because of the truncation of the ρₖ imposed by the prior (which is not the case in the Full or RS approach). Even more important is the indirect impact caused by the correlation of pulsar spin noise (red noise intrinsic to each pulsar) and common red noise at the lowest frequencies. In Figure 6, we show the inferred distribution for the parameters of the spin red noise in one of the best pulsars, PSR J1909-3744.

Extending the prior range for ρₖ leads to steeper and stronger spin red noise in that pulsar (dashed brown line compared to the solid brown line). This implies that the (artificially) over-constrained prior on ρₖ truncates the noise estimation in the pulsar J1909-3744. We empirically determined a fiducial value of power in terms of log₁₀ ρₖ (that can be roughly regarded as common to all frequency bins) which represents the threshold below which our data becomes insensitive. We selected a symmetric range centered around this specific value, which we set at log₁₀ ρₖ = -8. For the upper limit, we opted for -4, leading to a lower bound of -12 for the log₁₀ ρₖ prior. This prior answers direct (albeit partially) and indirect (by preserving correlations) impact on the Bayes factor while avoiding an overblow of the prior range.

In addition, it is worth noting <fnref n="5" /> that the posteriors given by the FS method typically have a broader spread compared to other methods. Moreover, the FS method encounters difficulties in accurately capturing the double peak feature observed in the Gμ posterior for the BOS model in Figure 4. This issue is mainly due to the width of the prior for the ρₖ parameters. Indeed, if a ρₖ posterior is not well-constrained, which is the case for several frequency bins as evident in Figure 5, taking a wider prior results in a longer tail and less constrained posterior at its upper bound (due to normalization). This gives a likelihood function (Equation 20) with less pronounced peaks, leading to the spread in the posterior distributions. By running the FS method using a narrower prior range of [-10, -4] for the ρₖ, we recover the posterior distribution which is much closer to the one obtained with the Full method.

In general, we have found the choice for the prior for the free spectrum evaluation (or rather a lack of rigorously defined cut) is a weak point of the FS approach in computing the evidence, despite that it perfectly falls into Bayesian philosophy.

Accumulating high-quality data and using more pulsars in the array, will hopefully result in the well-constrained power at low frequencies. Consequently, this will make the choice of the prior for the ρₖ parameters much less relevant for evaluation of the Bayes factors using the FS method.

## V. CONCLUSION AND DISCUSSION

In this work, we demonstrated that the Free Spectrum and Resampling methods are very powerful tools to determine rapidly the posterior distributions for the parameters describing the SGWB. Moreover, those methods (taking into account caveats for FS method discussed above) could be used to evaluate the Bayes factor between various SGWB models. We have demonstrated these methods using the 6-pulsar early DR2

<footnote n="5">⁵ We thank the referee for pointing us in this direction.</footnote>
<table><caption>TABLE II. For each Bayesian analysis method (columns), we write the cosmic string parameters posteriors for each loop distribution (BOS/LRS) with a varying number of kinks (.kk) or without. The 5% and 95% quantiles are used to set up the credible interval. The second column quotes the linear Bayes Factors comparing each of the cosmic string (CS) models against a power-law (PL) PSD for the SGWB.</caption><thead><tr><th rowspan="2"></th><th colspan="2">Full Method</th><th colspan="2">RS Method</th><th colspan="2">FS Method</th></tr><tr><th>CS Posterior</th><th>BF (PL/CS)</th><th>CS Posterior</th><th>BF (PL/CS)</th><th>CS Posterior</th><th>BF (PL/CS)</th></tr></thead><tbody><tr><td>BOS</td><td>log<sub>10</sub> Gμ = -10.08<sup>+0.32</sup><sub>-0.26</sub></td><td>0.3</td><td>-10.08<sup>+0.30</sup><sub>-0.27</sub></td><td>0.2</td><td>-9.95<sup>+0.56</sup><sub>-0.45</sub></td><td>0.1</td></tr><tr><td>LRS</td><td>log<sub>10</sub> Gμ = -10.60<sup>+0.17</sup><sub>-0.17</sub></td><td>0.2</td><td>-10.60<sup>+0.16</sup><sub>-0.18</sub></td><td>0.3</td><td>-10.59<sup>+0.27</sup><sub>-0.28</sub></td><td>0.2</td></tr><tr><td>BOS.kk</td><td>log<sub>10</sub> Gμ = -8.92<sup>+1.12</sup><sub>-1.42</sub><br>N<sub>k</sub> = 105<sup>+86</sup><sub>-96</sub></td><td>0.2</td><td>-8.98<sup>+1.17</sup><sub>-1.38</sub><br>101<sup>+90</sup><sub>-96</sub></td><td>0.3</td><td>-8.66<sup>+1.24</sup><sub>-1.65</sub><br>100<sup>+88</sup><sub>-87</sub></td><td>0.1</td></tr><tr><td>LRS.kk</td><td>log<sub>10</sub> Gμ = -10.90<sup>+0.29</sup><sub>-0.27</sub><br>N<sub>k</sub> = 102<sup>+89</sup><sub>-92</sub></td><td>0.2</td><td>-10.89<sup>+0.29</sup><sub>-0.27</sub><br>103<sup>+89</sup><sub>-95</sub></td><td>0.3</td><td>-10.86<sup>+0.35</sup><sub>-0.40</sub><br>104<sup>+85</sup><sub>-92</sub></td><td>0.2</td></tr></tbody></table>

EPTA dataset and comparing several models of SGWB produced by a network of cosmic strings.

For this 6-pulsars only dataset, we could perform rigorous inference (Full method) of SGWB parameters (in addition to pulsar noise parameters) and Bayes factors evaluation Table II. Note that the Full method took several days to complete, whereas the execution time for RS/FS methods is reduced to hours/minutes. Extending the data by including more pulsars, more back-ends and new observations will hugely explode the dimensionality of the problem and make the Full method computationally not tractable. The computational power scaling quadratically with the number of pulsars and the number of observations. We hope that this work convinces the reader of the validity of the fast and approximate methods which have to be used in the future PTA data analysis. The attractiveness of FS method is that it is not very sensitive to the number of pulsars in the array. Of course the evaluation of the free spectrum does heavily depend on the volume of the data and dimensionality of parameters space (pulsars noise models), but it can then be cheaply used to infer parameters for multiple SGWB models. The Resampling method is also relatively cheap, the sampling step scales linearly with number of pulsars in the array and can be efficiently parallelized (if needed). Even though we have demonstrated the fast methods using simplistic power-law model for SMBHB and two cosmic string SGWB models, the methods are generic and can be used to infer parameters of any SGWB.

The second result of this paper are the constraints on the cosmic string tension $Gμ$ for the BOS and LRS models. Assuming that the observed red noise process is GWs, we have obtained the preferred value for the string tension $Gμ ≈ 10^{-10.1}$ (resp. $10^{-10.6}$) for BOS (resp. LRS). Moreover, we find that the power-law model with HD correlations is slightly disfavored compared to the simpler one-parameter cosmic strings model. The posterior is not informative about the number of kinks when we consider the two-parameter ($Gμ, N_k$) cosmic string models. However, it also implies that the pos-

sibility of having a large number of kinks is not ruled out, though it might be disfavored on the theoretical grounds. Considering a two-component SGWB model (SMBHB+CS) shows a strong correlation between two components, and we set the upper bounds $log_{10} Gμ ≲ -9.9$ (BOS model) and $log_{10} Gμ ≲ -10.5$ at 95% (LRS model) confidence.

Compared to previous EPTA constraints [63, 64] on $log_{10} Gμ$ which were obtained using the EPTA sensitivity curve, here we use more up to date models and consider directly the PSD of the cosmic string SGWB through a wide range of frequencies to obtain our constraints. Our upper bounds on the string tension are well under the one given by CMB experiment $log_{10} Gμ ≲ -7$ [65] and the LVK collaboration [31] for the BOS model ($log_{10} Gμ ≲ -8$). However, for the LRS model, due to the population of small loops emitting at higher frequencies, the LVK constraint is more stringent, $log_{10} Gμ ≲ -14$.

This paper lays the path for the search and interpretation of a SGWB signal of any origin. The fast methods suggested here will be applied to the extended EPTA dataset containing 25 pulsars [5]. In addition to EPTA data, the new dataset which combines the most sensitive observations from all PTAs, International PTA data combination, is being produced. We expect the future data to be more sensitive to GWs and give better constraints on the string tension in BOS and LRS models.

## ACKNOWLEDGMENTS

H.Q-L, P.A. and D.S. would like to thank the University of Geneva for hospitality while this work was in progress. We are grateful to Chiara Caprini for very useful discussions. H.Q-L thanks Institut Polytechnique de Paris for funding his PhD. The work of P.A. is supported by the Wallonia-Brussels Federation Grant ARC № 19/24-103. S.B. and H.Q-L acknowledge support from ANR-21-CE31-0026, project MBH_waves. D.S. is grateful to CERN for hospitality.
[1] S. Chen et al., Monthly Notices of the Royal Astronomical Society 508, 4970 (2021), https://academic.oup.com/mnras/article-pdf/508/4/4970/40979667/stab2833.pdf.

[2] N. Arzoumanian et al, ApJ 905, L34 (2020), arXiv:2009.04496 [astro-ph.HE].

[3] P. c. B. Goncharov et al, The Astrophysical Journal Letters 917, L19 (2021).

[4] J. Antoniadis et al., Monthly Notices of the Royal Astronomical Society 510, 4873 (2022).

[5] J. Antoniadis et al. (EPTA), Astron. Astrophys. 678, A50 (2023), arXiv:2306.16214 [astro-ph.HE].

[6] G. Agazie et al., The Astrophysical Journal Letters 951, L8 (2023).

[7] D. J. Reardon, A. Zic, R. M. Shannon, G. B. Hobbs, M. Bailes, V. Di Marco, A. Kapur, A. F. Rogers, E. Thrane, J. Askew, N. D. R. Bhat, A. Cameron, M. Curyło, W. A. Coles, S. Dai, B. Goncharov, M. Kerr, A. Kulkarni, Y. Levin, M. E. Lower, R. N. Manchester, R. Mandow, M. T. Miles, R. S. Nathan, S. Osłowski, C. J. Russell, R. Spiewak, S. Zhang, and X.-J. Zhu, ApJ 951, L6 (2023), arXiv:2306.16215 [astro-ph.HE].

[8] H. Xu et al., Res. Astron. Astrophys. 23, 075024 (2023), arXiv:2306.16216 [astro-ph.HE].

[9] G. Agazie et al. (NANOGrav), Astrophys. J. Lett. 952, L37 (2023), arXiv:2306.16220 [astro-ph.HE].

[10] A. Sesana, Monthly Notices of the Royal Astronomical Society: Letters 433, L1 (2013).

[11] E. S. Phinney, "A practical theorem on gravitational wave backgrounds," (2001), arXiv:astro-ph/0108028 [astro-ph].

[12] A. Sesana, Classical and Quantum Gravity 30, 244009 (2013).

[13] J. Ellis and M. Lewicki, Phys. Rev. Lett. 126, 041304 (2021), arXiv:2009.06555 [astro-ph.CO].

[14] S. Blasi, V. Brdar, and K. Schmitz, Phys. Rev. Lett. 126, 041305 (2021), arXiv:2009.06607 [astro-ph.CO].

[15] L. Bian, J. Shu, B. Wang, Q. Yuan, and J. Zong, Phys. Rev. D 106, L101301 (2022), arXiv:2205.07293 [hep-ph].

[16] Z.-C. Chen, Y.-M. Wu, and Q.-G. Huang, Astrophys. J. 936, 20 (2022), arXiv:2205.07194 [astro-ph.CO].

[17] M. Hindmarsh and J. Kume, Journal of Cosmology and Astroparticle Physics 2023, 045 (2023).

[18] H. B. Nielsen and P. Olesen, Nucl. Phys. B 61, 45 (1973).

[19] T. W. B. Kibble, J. Phys. A 9, 1387 (1976).

[20] T. W. B. Kibble, Nucl. Phys. B 252, 227 (1985), [Erratum: Nucl.Phys.B 261, 750 (1985)].

[21] R. Jeannerot, J. Rocher, and M. Sakellariadou, Phys. Rev. D 68, 103514 (2003), arXiv:hep-ph/0308134.

[22] C. Ringeval, M. Sakellariadou, and F. Bouchet, JCAP 02, 023 (2007), arXiv:astro-ph/0511646.

[23] C. J. A. P. Martins and E. P. S. Shellard, Phys. Rev. D 73, 043515 (2006), arXiv:astro-ph/0511792.

[24] K. D. Olum and V. Vanchurin, Phys. Rev. D 75, 063521 (2007), arXiv:astro-ph/0610419.

[25] M. Hindmarsh, S. Stuckey, and N. Bevis, Phys. Rev. D 79, 123504 (2009), arXiv:0812.1929 [hep-th].

[26] D. Matsunami, L. Pogosian, A. Saurabh, and T. Vachaspati, Phys. Rev. Lett. 122, 201301 (2019), arXiv:1903.05102 [hep-ph].

[27] M. Hindmarsh, J. Lizarraga, A. Urio, and J. Urrestilla, Phys. Rev. D 104, 043519 (2021), arXiv:2103.16248 [astro-ph.CO].

[28] P. Auclair, D. A. Steer, and T. Vachaspati, Phys. Rev. D 101, 083511 (2020), arXiv:1911.12066 [hep-ph].

[29] P. Auclair, K. Leyde, and D. A. Steer, JCAP 04, 005 (2023), arXiv:2112.11093 [astro-ph.CO].

[30] B. P. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev. D 97, 102002 (2018), arXiv:1712.01168 [gr-qc].

[31] R. Abbott et al. (LIGO Scientific, Virgo, KAGRA), Phys. Rev. Lett. 126, 241102 (2021), arXiv:2101.12248 [gr-qc].

[32] P. Auclair, J. J. Blanco-Pillado, D. G. Figueroa, A. C. Jenkins, M. Lewicki, M. Sakellariadou, S. Sanidas, L. Sousa, D. A. Steer, J. M. Wachter, and S. Kuroyanagi, Journal of Cosmology and Astroparticle Physics 2020, 034 (2020).

[33] P. Auclair et al. (LISA Cosmology Working Group), (2022), arXiv:2204.05434 [astro-ph.CO].

[34] R. W. Hellings and G. S. Downs, The Astrophysical Journal 265 (1983).

[35] A. Chalumeau et al., Monthly Notices of the Royal Astronomical Society 509, 5538 (2021).

[36] R. van Haasteren and Y. Levin, Monthly Notices of the Royal Astronomical Society 428, 1147 (2012).

[37] R. van Haasteren and M. Vallisneri, Physical Review D 90 (2014), 10.1103/physrevd.90.104012.

[38] S. Hourihane, P. Meyers, A. Johnson, K. Chatziioannou, and M. Vallisneri, Phys. Rev. D 107, 084045 (2023), arXiv:2212.06276 [gr-qc].

[39] S. R. Taylor, J. Simon, L. Schult, N. Pol, and W. G. Lamb, Physical Review D 105 (2022), 10.1103/physrevd.105.084049.

[40] W. G. Lamb, S. R. Taylor, and R. van Haasteren, "The need for speed: Rapid refitting techniques for bayesian spectral characterization of the gravitational wave background using ptas," (2023), arXiv:2303.15442 [astro-ph.HE].

[41] C. J. Moore and A. Vecchio, Nature Astronomy 5, 1268 (2021).

[42] J. J. Blanco-Pillado, K. D. Olum, and B. Shlaer, Physical Review D 89 (2014), 10.1103/physrevd.89.023512.

[43] T. Damour and A. Vilenkin, Physical Review D 64 (2001), 10.1103/physrevd.64.064008.

[44] C. Ringeval and T. Suyama, Journal of Cosmology and Astroparticle Physics 2017, 027 (2017).

[45] P. Binetruy, A. Bohe, T. Hertog, and D. A. Steer, Phys. Rev. D 82, 083524 (2010), arXiv:1005.2426 [hep-th].

[46] P. Binetruy, A. Bohe, T. Hertog, and D. A. Steer, Phys. Rev. D 82, 126007 (2010), arXiv:1009.2484 [hep-th].

[47] A. Vilenkin and E. Shellard, *Cosmic Strings and Other Topological Defects*, Cambridge Monographs on Mathematical Physics (Cambridge University Press, 1994).

[48] N. Aghanim et al. (Planck), Astron. Astrophys. 641, A6 (2020), [Erratum: Astron. Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO].

[49] P. Auclair et al., JCAP 04, 034 (2020), arXiv:1909.00819 [astro-ph.CO].

[50] J. J. Blanco-Pillado and K. D. Olum, Physical Review D 96 (2017), 10.1103/physrevd.96.104046.

[51] L. Lorenz, C. Ringeval, and M. Sakellariadou, Journal of Cosmology and Astroparticle Physics 2010, 003 (2010).

[52] P. G. Auclair, JCAP 11, 050 (2020), arXiv:2009.00334 [astro-ph.CO].
[astro-ph.CO].

[53] J. A. Ellis, M. Vallisneri, S. R. Taylor, and P. T. Baker, "ENTERPRISE: Enhanced Numerical Toolbox Enabling a Robust Pulsar Inference SuitE," (2020).

[54] J. Ellis and R. van Haasteren, "jellis18/ptmcmc sampler: Official release," (2017).

[55] S. Hee, W. J. Handley, M. P. Hobson, and A. N. Lasenby, Monthly Notices of the Royal Astronomical Society **455**, 2461 (2015).

[56] Z. Arzoumanian et al., The Astrophysical Journal **859**, 47 (2018).

[57] D. Scott, *Multivariate Density Estimation: Theory, Practice, and Visualization*, Wiley Series in Probability and Statistics (Wiley, 2015).

[58] J. Skilling, in *Bayesian Inference and Maximum Entropy Methods in Science and Engineering: 24th International Workshop on Bayesian Inference and Maximum Entropy Methods in Science and Engineering*, American Institute of Physics Conference Series, Vol. 735, edited by R. Fischer, R. Preuss, and U. V. Toussaint (2004) pp. 395-405.

[59] J. Skilling, Bayesian Analysis **1**, 833 (2006).

[60] S. Koposov, J. Speagle, K. Barbary, G. Ashton, E. Ben-nett, J. Buchner, C. Scheffler, B. Cook, C. Talbot,

J. Guillochon, P. Cubillos, A. A. Ramos, B. Johnson, D. Lang, Ilya, M. Dartiailh, A. Nitz, A. McCluskey, and A. Archibald, "joshspeagle/dynesty: v2.1.3," (2023).

[61] J. S. Speagle, Monthly Notices of the Royal Astronomical Society **493**, 3132-3158 (2020).

[62] J. J. Blanco-Pillado, K. D. Olum, and B. Shlaer, Physical Review D **92** (2015), 10.1103/physrevd.92.063528.

[63] L. Lentati, S. R. Taylor, C. M. F. Mingarelli, A. Sesana, S. A. Sanidas, A. Vecchio, R. N. Caballero, K. J. Lee, R. van Haasteren, S. Babak, C. G. Bassa, P. Brem, M. Burgay, D. J. Champion, I. Cognard, G. Desvignes, J. R. Gair, L. Guillemot, J. W. T. Hessels, G. H. Janssen, R. Karuppusamy, M. Kramer, A. Lassus, P. Lazarus, K. Liu, S. Oslowski, D. Perrodin, A. Petiteau, A. Posenti, M. B. Purver, P. A. Rosado, R. Smits, B. Stappers, G. Theureau, C. Tiburzi, and J. P. W. Verbiest, Monthly Notices of the Royal Astronomical Society **453**, 2577 (2015).

[64] S. A. Sanidas, R. A. Battye, and B. W. Stappers, Physical Review D **85** (2012), 10.1103/physrevd.85.122003.

[65] T. Charnock, A. Avgoustidis, E. J. Copeland, and A. Moss, Physical Review D **93** (2016), 10.1103/physrevd.93.123503.