Searching for cosmic string induced stochastic gravitational wave background with the Parkes Pulsar Timing Array

We search for stochastic gravitational wave background emitted from cosmic strings using the Parkes Pulsar Timing Array data over 15 years. While we find that the common power-law excess revealed by several pulsar timing array experiments might be accounted for by the gravitational wave background from cosmic strings, the lack of the characteristic Hellings-Downs correlation cannot establish its physical origin yet. The constraints on the cosmic string model parameters are thus derived with conservative assumption that the common power-law excess is due to unknown background. Two representative cosmic string models with different loop distribution functions are considered. We obtain constraints on the dimensionless string tension parameter $G\mu<10^{-11}\sim10^{-10}$, which is more stringent by two orders of magnitude than that obtained by the high-frequency LIGO-Virgo experiment for one model, and less stringent for the other. The results provide the chance to test the Grand unified theories, with the spontaneous symmetry breaking scale of $U(1)$ being two-to-three orders of magnitude below $10^{16}$ GeV. The pulsar timing array experiments are thus quite complementary to the LIGO-Virgo experiment in probing the cosmic strings and the underlying beyond standard model physics in the early Universe.

We search for stochastic gravitational wave background emitted from cosmic strings using the Parkes Pulsar Timing Array data over 15 years. While we find that the common power-law excess revealed by several pulsar timing array experiments might be accounted for by the gravitational wave background from cosmic strings, the lack of the characteristic Hellings-Downs correlation cannot establish its physical origin yet. The constraints on the cosmic string model parameters are thus derived with conservative assumption that the common power-law excess is due to unknown background. Two representative cosmic string models with different loop distribution functions are considered. We obtain constraints on the dimensionless string tension parameter Gµ < 10 −11 ∼ 10 −10 , which is more stringent by two orders of magnitude than that obtained by the high-frequency LIGO-Virgo experiment for one model, and less stringent for the other. The results provide the chance to test the Grand unified theories, with the spontaneous symmetry breaking scale of U(1) being two-to-three orders of magnitude below 10 16 GeV. The pulsar timing array experiments are thus quite complementary to the LIGO-Virgo experiment in probing the cosmic strings and the underlying beyond standard model physics in the early Universe.

I. INTRODUCTION
Cosmic strings (CS) are one-dimensional topological defects supposed to be formed during phase transitions where symmetry gets broken spontaneously [1,2]. The CS dynamics can be described by the Nambu-Goto action for thin and local strings with no internal structures. In this situation, infinite strings can reach the scaling regime [3][4][5] and go to loops through the intercommutation of intersecting string segments [6]. The formed small loops will oscillate and emit gravitational wave (GW) bursts by the structures of Cusps and Kinks [7,8]. The superposition of uncorrelated GW bursts from many CSs form a stochastic gravitational wave background (SGWB). In the Nambu-Goto string scenario, the SGWB is characterized by the loop number density and the string tension (µ) [9]. The dimensionless parameter Gµ (G is the Newtonian constant) that parameterized the gravitational interactions of strings is usually adopted in literature, which is tightly connect with the GUT scale (M GUT ∼ 10 16 GeV) as Gµ ∼ (η/M GUT ) 2 since CSs are generally predicted in the symmetry breaking chain of grand unified theories (GUT). Therefore, the detection of SGWB from CSs provide an intriguing way to access the beyond-standard-model physics close to the GUT scale that are inaccessible by high-energy colliders [10][11][12][13], such as leptogenesis and the type-I seesaw scale [14].
The gravitational waves spectra from CSs span in a wide frequency range characterized by a plateau in the highfrequency region. The SGWB from CSs is one of the most promising target of LIGO-Virgo [15,16], Laser Interferometer Space Antenna (LISA) [17], and pulsar timing arrays (PTA) [18,19]. Recently, both LIGO-Virgo and PTA experiments made great progresses on the search of GWs. The LIGO-Virgo group place stringent constraints on the GW bursts and the SGWB from CSs in the high-frequency window [15]. In the nanohertz range, several PTAs reported the detection of a mysterious common red process [20][21][22][23], which may be explained as SGWB from various sources including CSs [24][25][26][27]. However, the GW interpretation of the common red process is still suspicious due to the lack of the characteristic Hellings-Downs (HD) correlation in the data. In this Letter we search for SGWB signals generated from CSs utilizing the Parkes Pulsar Timing Array (PPTA) data. We consider two representative SGWB models of CSs with two different loop distribution functions, and derive the constraints on the dimensionless Gµ parameter of CS. The results can be translated to constraints on the spontaneous symmetry breaking scale of the GUT.

II. SGWB SPECTRA FOR COSMIC STRING NETWORKS
The SGWB from cosmic string networks comes from the uncorrelated superpositions of GW bursts of three contributions: cusps, kinks, and kink-kink collision. The SGWB spectrum emitted from cosmic strings is given by where ρ c = 3H 2 0 8πG is the critical energy density of the universe, dρ GW d f (t 0 , f ) is the GWs energy density per unit frequency today. Considering different modes of loops oscillation (n), we have where the C n ( f ) is a function of loop distributions depending on the cosmological background. To account for all contributions from cusps, kinks, and kink-kink collisions, we adopt the P n from numerical simulations [28]. Analytically, P n = Γ ζ(q) n −q , Γ ≈ 50 and ζ(q) is the Riemann zeta function with q = 4/3, 5/3, and 2 representing for cusps, kinks, and kink-kink collisions, respectively.
We first consider the SGWB from CSs with the loop production functions for non-self-intersecting loops being obtained directly from CS networks simulation in radiation and matter dominated era 1 by Blanco-Pillado, Olum, Shlaer [30,31] (hereafter denoted as BOS model, which is dubbed as model A by LIGO-Virgo group in Ref. [15]). In the model, one need to take into account contributions from the radiation dominated era and the matter dominated era with the loop number density n(l, t) given in the Supplemental Material. For the SGWB contributions from the radiation era, one has where z eq is the redshift in the radiation-matter equality and z cut is the cutoff redshift. The subscript r here denotes radiation dominated era, similarly, in the following rm is for the case where loops are formed in radiation dominated era but survive past to matter dominated era and m for the case of matter dominated era. We note that here the n r (l, t) is the loop distribution function rather than the oscillation mode n before the integration. The SGWB contributions from matter dominated era constitute two parts: loops that survive to matter dominated era and loops formed in matter dominated era. The C n takes the form of with i being rm and m for these two cases, respectively. Following the same procedure, we can calculate the SGWB spectrum with loop distribution functions being derived analytically by Lorenz, Ringeval, and Sakellariadou [32] (denoted as LRS, which is dubbed as model B by LIGO-Virgo group in Ref. [15]), where the distribution of non-selfintersecting scaling loops are extracted from simulations [33]. In the model, the contribution from the radiation dominated era takes the form of and in the matter dominated era, the GW spectrum is  Tables I and II, respectively). The sensitive frequency regions of PTA and LIGO-Virgo are denoted by green and blue shaded bands. Fig. 1 shows the total SGWB spectra of the CS-induced SGWB for the BOS and LRS models, for a few values of Gµ. It is shown that the BOS model predicts relatively flat spectra (solid lines) and hence the PTA experiment is more suitable to probe it. On the other hand, the LRS model predicts harder spectra (dashed lines) which is more optimal for the highfrequency GW experiment such as LIGO-Virgo. Also shown in Fig. 1 are the upper limits of the free-spectrum SGWB derived from the PPTA data [34]. We find that the CS models with Gµ of 10 −8 ∼ 10 −10 are severely constrained by the PPTA data. The magnitude of the SGWB spectra for the LRS model is much higher than that of the BOS model at high frequencies due to more small loops in the loop distribution function dominate in such frequency range.

III. DATA ANALYSIS
The dataset we used in this work is the second data release (DR2) of the PPTA [35], which is available in the CSIRO pulsar data archive 2 . Searches for SGWBs and other fundamental physics problems using these data (or the subsets) were also carried out [21,34,[36][37][38]. In PTA, to search for such an SGWB signal from CS is to find spatially correlated time residuals among time-of-arrivals (ToA) of different pulsars. The residuals are composed of deterministic timing models, noises, and the hypothetical signal (SGWB from CS here). We use the TEMPO2 tool [39,40] to fit the timing models of pulsars, the ENTERPRISE 3 and the ENTERPRISE EXTENSIONS 4 to model the noises, and PTMCMCSampler [41] to do the Bayesian analysis. The noise model is based on single pulsar analyses of Ref. [42]. Briefly speaking, the stochastic noises mainly include two parts, the white noise and red noise. The white noise may come from the radio frequency interference, pulse profile changes or instrumental artifacts. We use three parameters, EFAC (Error FACtor), EQUAD (Error added in QUADrature) and ECORR (Error of CORRelation between ToAs in a single epoch), to account for white noise which can not be subtracted by fitting the timing model. As usually done [37,38], we fix the white noise parameters as their maximum likelihood values from the single pulsar analyses in the Bayesian analysis. The red noises are mainly caused by the irregularities of the pulsar spin (spin noise) and the dispersion measure of photons when traveling through the interstellar medium (DM noise). For some pulsars there are band noise for ToAs of a certain photon frequency band and chromatic noise which correlates between different photon frequencies.
All the red noises are modeled as power-law forms with amplitude A and slope γ. Several analyses revealed that there is a common power-law (CPL) process in the pulsar ToAs, whose nature is still in debate [20][21][22][23]. In this work we will test different assumptions on the CPL, as either the SGWB signal from CS (see also [24,25]) or an unknown background. The solar system ephemeris uncertainties are modeled with a 11parameter model BayesEphem implemented in ENTERPRISE, including perturbations in masses of major planets, the drift rate of the Earth-Moon barycenter orbit, and the perturbations of the Earth's orbit from Jupiter's average orbital elements described by 6 parameters [43]. We summarize the noise and signal parameters together with their priors adopted in the Bayesian analysis in Table S1 in the Supplemental Material.

IV. RESULT
To address how significant the CS-induced SGWB signals in the data, we test the Bayes factor (BF) of the "signal hy-pothesis" against the null hypothesis, following the Savage-Dickey formula [44] where BF 10 means the BF of hypothesis H1 against H0, D is the observational data, φ is the parameters of the signal model, and φ 0 is the parameters of the null hypothesis which is a subset of φ. P 0 and P 1 are the evidence of the noise and signal hypotheses. The BF 10 is equivalent to the ratio of the prior to the posterior probabilities of the null hypothesis. The null hypothesis (H0) corresponds to the model with only the pulsar timing model and noise. For hypothesis H1, we assume an additional CPL process in the model. The CS model with the HD correlation is included in substitution of the CPL process is labelled as H2. In addition, in H3 and H4 we consider simultaneously the CPL process as an unknown systematics and the CS contribution, in which the autocorrelation of a pulsar's own ToAs is not subtracted (H3) and subtracted (H4), respectively. The results of the fittings for different hypotheses are summarized in Tables I and II. For hypothesis H1, a clear CPL signal with a BF of 10 3.2 is revealed in the data, consistent with previous studies [20][21][22][23]. The posterior distributions of the CPL model parameters are shown in Fig. S1 in the Supplemental Material. Similar signal with comparable BF value is found if we assume a CS-induced SGWB component in the model rather than the CPL (H2). Fitting results of the CS parameter log 10 Gµ are −10.38 ± 0.21 and −10.89 +0.14 −0.17 for the BOS and LRS model, with 1σ error bars presented. The corresponding parameter values are consistent with that employed to account for the NANOGrav data, but with different model assumptions [24]. Since the nature of the CPL process is unclear, and particularly the characteristic HD correlation of GWs is still lack, we also test the hypotheses that treating the CPL as an unknown background. The resulting BF of hypotheses H3 or H4 against H1 is close to 1 which means that the evidence of the CS is not significant in case of a CPL background. We therefore derive the 95% credible level (C.L.) upper limits on log 10 Gµ, as given in Tables I and II. The posterior distributions of the CPL and CS parameters for H4 are shown in Fig. 2. More results of the analysis based on different hypotheses can be found in Figs. S2-S3 in the Supplemental Material. Additional tests assuming that the CPL has an astrophysical origin from the supermassive binary black hole (SMBBH) give similar conclusion (see the Supplemental Material).
Our results can be compared with those obtained by the LIGO-Virgo observations of high-frequency GWs [15], which is displayed in Fig. 3. We also show the bound from the cosmic microwave background (CMB), while the less stronger limit from Big Bang nucleosynthese (BBN) is too weak to be shown in the figure [15,45]. The LIGO-Virgo upper limits of Gµ are 10 −8 − 10 −6 and (4.0 − 6.3) × 10 −15 for the BOS and LRS model (model A and B in [15]), respectively. Compared with LIGO-Virgo results, the PPTA upper limits on log 10 Gµ improve by more than two orders of magnitude for the BOS model. For the LRS model, the LIGO-Virgo result is more stringent because the SGWB spectrum is harder. The PTA ex-   Constraints on Gµ at the 95% C.L. from the PPTA data for different hypotheses, compared with those obtained by LIGO-Virgo experiment [15]. Silver boxes represent the ranges reported in the LIGO-Virgo analysis. The limit from cosmic microwave background (CMB) is also reported by the black dashed line [45]. than 15 years of observations of 26 millisecond pulsars by the PPTA experiment to search for SGWB from CS networks. We show that the CPL process revealed by several PTAs recently can be explained by the CS-induced SGWB signal with log 10 Gµ ∼ −10.38 and −10.89 for the two typical classes of CS models, the BOS and LRS models. While the LRS model explanation would be excluded by the LIGO-Virgo observations at high frequencies, the BOS model explanation is consistent with the LIGO-Virgo data [15]. As a conservative alternative, we also assume the CPL as a background component and turn to set limits on the CS model parameters. The PPTA upper limits on log 10 Gµ reach about −10 ∼ −11, depending on the CS models and the analysis methods. These constraints correspond to the symmetry breaking scale η ∼ 10 13 − 10 14 GeV, that challenge the grand unification theories below the GUT scale already. For the BOS model, the PPTA limits are more stringent than those from LIGO-Virgo, and for the LRS model the LIGO-Virgo experiment is more sensitive. This is due to the fact that, in the LRS model, more small loops are contained in the loop distribution functions that control the behaviours of the SGWB spectra from CS networks in high frequency range. We expect that continuous accumulation of more precise data by PTAs around the world in the near future will critically test the CS models and more generally, a class of GUTs admitting the U(1) symmetry breaking in the early Universe.

A. CS loop distributions
The BOS model. -In this model, the loop production functions are inferred from Nambu-Goto simulations of CS networks in radiation and matter dominated era directly. In the radiation dominated era, we take l/t 0.1 to consider the cutoff of the maximum size of loops. In the matter dominated era, as suggested by simulations, plenty of loops formed in the radiation era will survive till the time of radiation-matter equality and emit GWs in the matter dominated era, which imposes a constrain on loop size, l/t < 0.09t eq /t − ΓGµ. Furthermore, loops can also be produced when the CS networks reach the scaling regime in the matter dominated era, in which case the loop size should satisfy l/t < 0.18. The loop distribution functions of all the circumstances are n r (l, t) = 0.18 t 3/2 (l + ΓGµt) 5/2 , l/t 0.1, (S1) The LRS model. -Different from the BOS model, in the LRS model, the loop distribution of the non-self-intersecting scaling loops rather than the loop production function is inferred from CS network simulations [32]. On scales x = l/t ΓGµ ≡ x d , the simulation gives where the two constants C 0 and p in the radiation dominated era and matter dominated era are Note that, not only loops would emit GWs -which decreases their length l -but the GW emission also back-reacts on the loops. The back-reaction smooths out the loops on the smallest scales (in particular the kinks), hindering the formation of smaller loops. To extend the loop distribution given above down to smaller scales, a further length scale x c x d (the so-called "gravitational back-reaction scale") is introduced, which is estimated as x c = 20(Gµ) 1+2χ through matching the loop distribution on scales x x d (Eq. (S4)), where χ r = 0.2 +0.07 −0.10 for the radiation dominated era and χ m = 0.295 +0.03 −0.04 for the matter dominated era [17]. The specific distributions in all periods are then given by [32] n(x > x d ) where C = C 0 (1 − ν) 2−p , and ν = 1/2 for radiation dominated era and ν = 2/3 for matter dominated era.   Then we search for the signal of SGWB from the CS. Fig.  S2 shows the one-dimensional probability distributions of the CS parameter Gµ without including the CPL process (H2). This test reveals that the CS-induced SGWB can be used to explain the CPL excess in the data, with similar BFs compared with the CPL assumption. Fig. S3 gives the posterior distributions of the CPL parameters and the CS parameter for hypothesis H3. In this case, the CS signal is not evident, and we can derive constraints on the CS parameter. At the 95% C.L., we find log 10 Gµ < −10.02 (−10.64) for the BOS (LRS) model. Similar analysis but re- moving the auto-correlation of the pulsars' ToAs (H4) is pre-sented in the main text. Finally we consider the possibility that the CPL has an astrophysical origin from the supermassive binary black holes (SMBBH). The spectrum of the SMBBH SGWB signal is given by

B. PTA model parameters
where the power-law index γ SMBBH is fixed to be 13/3. The prior of log 10 A SMBBH is assumed to be uniformly distributed in [−18, −14]. We test three model assumptions, with H5 being the SMBBH-only hypothesis and H6 (H7) being the SMBBH + CS for the BOS (LRS) model. The constraints on the model parameters of hypotheses H5 -H7 are shown in Fig. S4, Fig. S5, and Table S2. The CPL process may also be explained by the SMBBH signal, with log 10 A SMBBH = −14.89 +0.10 −0.12 . For hypotheses H6 and H7, we reach similar conclusion with H3 or H4 where the CPL is  treated as a background but with slight changes of the 95% C.L. upper limits of log 10 Gµ.