The 15 yr pulsar timing data set collected by the North American Nanohertz Observatory for Gravitational Waves NANOGrav shows positive evidence for the presence of a lowfrequency gravitationalwave GW background. In this paper we investigate potential cosmological interpretations of this signal specifically cosmic inflation scalarinduced GWs firstorder phase transitions cosmic strings and domain walls. We find that with the exception of stable cosmic strings of field theory origin all these models can reproduce the observed signal. When compared to the standard interpretation in terms of inspiraling supermassive black hole binaries SMBHBs many cosmological models seem to provide a better fit resulting in Bayes factors in the range from 10 to 100. However these results strongly depend on modeling assumptions about the cosmic SMBHB population and at this stage should not be regarded as evidence for new physics. Furthermore we identify excluded parameter regions where the predicted GW signal from cosmological sources significantly exceeds the NANOGrav signal. These parameter constraints are independent of the origin of the NANOGrav signal and illustrate how pulsar timing data provide a new way to constrain the parameter space of these models. Finally we search for deterministic signals produced by models of ultralight dark matter ULDM and dark matter substructures in the Milky Way. We find no evidence for either of these signals and thus report updated constraints on these models. In the case of ULDM these constraints outperform torsion balance and atomic clock constraints for ULDM coupled to electrons muons or gluons.
['Gravitational waves', 'Cosmology', 'Particle astrophysics', 'Gravitational wave sources', '678', '343', '96', '677', 'Astrophysics - High Energy Astrophysical Phenomena', 'Astrophysics - Cosmology and Nongalactic Astrophysics', 'General Relativity and Quantum Cosmology', 'High Energy Physics - Phenomenology']
The NANOGrav 15 yr Data Set Search for Signals from New Physics
{'The NANOGrav 15-year Data Set: Search for Signals from New Physics': "Adeela Afzal, 1, 2 Gabriella Agazie, 3 Akash Anumarlapudi, 3 Anne M. Archibald, 4 Zaven Arzoumanian, 5 Paul T. Baker, 6 Bence B'ecsy, 7 Jose Juan Blanco-Pillado, 8, 9, 10 Laura Blecha, 11 Kimberly K. Boddy, 12 Adam Brazier, 13, 14 Paul R. Brook, 15 Sarah Burke-Spolaor, 16, 17 Rand Burnette, 7 Robin Case, 7 Maria Charisi, 18 Shami Chatterjee, 13 Katerina Chatziioannou, 19 Belinda D. Cheeseboro, 16, 17 Siyuan Chen, 20 Tyler Cohen, 21 James M. Cordes, 13 Neil J. Cornish, 22 Fronefield Crawford, 23 H. Thankful Cromartie, 13, ∗ Kathryn Crowter, 24 Curt J. Cutler, 25, 19 Megan E. DeCesar, 26 Dallas DeGan, 7 Paul B. Demorest, 27 Heling Deng, 7 Timothy Dolch, 28, 29 Brendan Drachler, 30, 31 Richard von Eckardstein, 32 Elizabeth C. Ferrara, 33, 34, 35 William Fiore, 16, 17 Emmanuel Fonseca, 16, 17 Gabriel E. Freedman, 3 Nate Garver-Daniels, 16, 17 Peter A. Gentile, 16, 17 Kyle A. Gersbach, 18 Joseph Glaser, 16, 17 Deborah C. Good, 36, 37 Lydia Guertin, 38 Kayhan Gultekin, 39 Jeffrey S. Hazboun, 7 Sophie Hourihane, 19 Kristina Islo, 3 Ross J. Jennings, 16, 17, † Aaron D. Johnson, 3, 19 Megan L. Jones, 3 Andrew R. Kaiser, 16, 17 David L. Kaplan, 3 Luke Zoltan Kelley, 40 Matthew Kerr, 41 Joey S. Key, 42 Nima Laal, 7 Michael T. Lam, 30, 31 William G. Lamb, 18 T. Joseph W. Lazio, 25 Vincent S. H. Lee, 19 Natalia Lewandowska, 43 Rafael R. Lino dos Santos, 44, 32 Tyson B. Littenberg, 45 Tingting Liu, 16, 17 Duncan R. Lorimer, 16, 17 Jing Luo, 46, ‡ Ryan S. Lynch, 47 Chung-Pei Ma, 40, 48 Dustin R. Madison, 49 Alexander McEwen, 3 James W. McKee, 50, 51 Maura A. McLaughlin, 16, 17 Natasha McMann, 18 Bradley W. Meyers, 24, 52 Patrick M. Meyers, 19 Chiara M. F. Mingarelli, 37, 36, 53 Andrea Mitridate, 54 Jonathan Nay, 12 Priyamvada Natarajan, 55, 56 Cherry Ng, 57 David J. Nice, 58 Stella Koch Ocker, 13 Ken D. Olum, 59 Timothy T. Pennucci, 60 Benetge B. P. Perera, 61 Polina Petrov, 18 Nihan S. Pol, 18 Henri A. Radovan, 62 Scott M. Ransom, 63 Paul S. Ray, 41 Joseph D. Romano, 64 Shashwat C. Sardesai, 3 Ann Schmiedekamp, 65 Carl Schmiedekamp, 65 Kai Schmitz, 32 Tobias Schroder, 32 Levi Schult, 18 Brent J. Shapiro-Albert, 16, 17, 66 Xavier Siemens, 7, 3 Joseph Simon, 67, § Magdalena S. Siwek, 68 Ingrid H. Stairs, 24 Daniel R. Stinebring, 69 Kevin Stovall, 27 Peter Stratmann, 32 Jerry P. Sun, 7 Abhimanyu Susobhanan, 3 Joseph K. Swiggum, 58, † Jacob Taylor, 7 Stephen R. Taylor, 18 Tanner Trickle, 70 Jacob E. Turner, 16, 17 Caner Unal, 71, 72 Michele Vallisneri, 25, 19 Sonali Verma, 73, 74 Sarah J. Vigeland, 3 Haley M. Wahl, 16, 17 Qiaohong Wang, 18 Caitlin A. Witt, 75, 76 David Wright, 77 Olivia Young, 30, 31 and Kathryn M. Zurek 78  \n1 Institute for Theoretical Physics, University of Munster, D-48149 Munster, Germany  \n2 Department of Physics, Quaid-i-Azam University, Islamabad, 45320, Pakistan  \n3 Center for Gravitation, Cosmology and Astrophysics, Department of Physics, University of Wisconsin-Milwaukee, P.O. Box 413, Milwaukee, WI 53201, USA  \n4 Newcastle University, NE1 7RU, UK  \n- 5 X-Ray Astrophysics Laboratory, NASA Goddard Space Flight Center, Code 662, Greenbelt, MD 20771, USA\n- 6 Department of Physics and Astronomy, Widener University, One University Place, Chester, PA 19013, USA  \n7 Department of Physics, Oregon State University, Corvallis, OR 97331, USA  \n- 8 Department of Physics, University of Basque Country, UPV/EHU, E-48080, Bilbao, Spain  \n9 EHU Quantum Center, University of Basque Country, UPV/EHU, Bilbao, Spain  \n10 IKERBASQUE, Basque Foundation for Science, E-48011, Bilbao, Spain  \n- 11 Physics Department, University of Florida, Gainesville, FL 32611, USA\n- 12 Department of Physics, The University of Texas at Austin, Austin, TX 78712, USA\n- 13 Cornell Center for Astrophysics and Planetary Science and Department of Astronomy, Cornell University, Ithaca, NY 14853, USA 14 Cornell Center for Advanced Computing, Cornell University, Ithaca, NY 14853, USA\n- 15 Institute for Gravitational Wave Astronomy and School of Physics and Astronomy, University of Birmingham, Edgbaston, Birmingham B15 2TT, UK  \n16 Department of Physics and Astronomy, West Virginia University, P.O. Box 6315, Morgantown, WV 26506, USA  \n- 17 Center for Gravitational Waves and Cosmology, West Virginia University, Chestnut Ridge Research Building, Morgantown, WV 26505, USA\n- 18 Department of Physics and Astronomy, Vanderbilt University, 2301 Vanderbilt Place, Nashville, TN 37235, USA\n- 19 Division of Physics, Mathematics, and Astronomy, California Institute of Technology, Pasadena, CA 91125, USA\n- 20 Kavli Institute for Astronomy and Astrophysics, Peking University, Beijing, 100871, People's Republic of China  \nCorresponding author: The NANOGrav Collaboration  \ncomments@nanograv.org", 'the NANOGrav Collaboration': 'Table 4 . Bayesian Estimators, Maximum Posterior Values, and 68% Credible Intervals for the Parameters of the New-physics Models. Values annotated with ∗ are at the boundary of the prior range used in the analysis.', 'NANOGrav 15-year New-Physics Signals': 'Department of Physics, University of Central Florida, Orlando, FL 32816-2385, USA 78 Walter Burke Institute for Theoretical Physics California Institute of Technology, Pasadena, CA 91125 USA', 'ABSTRACT': 'The 15-year pulsar timing data set collected by the North American Nanohertz Observatory for Gravitational Waves (NANOGrav) shows positive evidence for the presence of a low-frequency gravitationalwave (GW) background. In this paper, we investigate potential cosmological interpretations of this signal, specifically cosmic inflation, scalar-induced GWs, first-order phase transitions, cosmic strings, and domain walls. We find that, with the exception of stable cosmic strings of field theory origin, all these models can reproduce the observed signal. When compared to the standard interpretation in terms of inspiraling supermassive black hole binaries (SMBHBs), many cosmological models seem to provide a better fit resulting in Bayes factors in the range from 10 to 100. However, these results strongly depend on modeling assumptions about the cosmic SMBHB population and, at this stage, should not be regarded as evidence for new physics. Furthermore, we identify excluded parameter regions where the predicted GW signal from cosmological sources significantly exceeds the NANOGrav signal. These parameter constraints are independent of the origin of the NANOGrav signal and illustrate how pulsar timing data provide a new way to constrain the parameter space of these models. Finally, we search for deterministic signals produced by models of ultralight dark matter (ULDM) and dark matter substructures in the Milky Way. We find no evidence for either of these signals and thus report updated constraints on these models. In the case of ULDM, these constraints outperform torsion balance and atomic clock constraints for ULDM coupled to electrons, muons, or gluons.  \nKeywords: Gravitational waves - Cosmology: early universe - Methods: data analysis', '1. INTRODUCTION': "The Standard Model (SM) of particle physics currently provides our best description of the laws governing the universe at subatomic scales. However, it fails to explain several observed properties of our universe, such as the origin of the matter-antimatter asymmetry, the nature of dark matter (DM) and dark energy, and the origin of neutrino masses. These shortcomings have motivated the development of several theories for physics beyond the SM, or BSM theories for short, accompanied by a rich experimental program trying to test them. The generation of gravitational waves (GWs) is a ubiquitous feature of many BSM theories (Maggiore 2000; Caprini & Figueroa 2018; Christensen 2019). These GWs form a stochastic background and propagate essentially unimpeded over cosmic distances to be detected today, whereas electromagnetic radiation does not start  \nfree streaming until after recombination. Thus, detecting a stochastic GW background (GWB) of cosmological origin would offer a unique and direct glimpse into the very early universe and herald a new era for using GWs to study fundamental physics.  \nCosmological GWBs can be produced by a number of particle physics models of the early universe. Notably, cosmic inflation generically produces GWs (Guzzetti et al. 2016), which may be observable at nanohertz frequencies if their energy density spectrum is sufficiently blue-tilted. Similarly, an enhanced spectrum of shortwavelength scalar perturbations produced during inflation can source so-called scalar-induced GWs (SIGWs; Dom'enech 2021; Yuan & Huang 2021a). Another potential source of GWs are cosmological first-order phase transitions (Caprini et al. 2016, 2020; Hindmarsh et al. 2021), which proceed through bubble nucleation; bubble collisions and bubble interactions with the primordial plasma giving rise to sound waves contribute to GW production. Finally, topological defects left behind by cosmological phase transitions, such as cosmic strings and domain walls (Vilenkin 1985; Hindmarsh & Kibble 1995; Saikawa 2017), can radiate GWs and hence contribute to the GWB.  \nThe North American Nanohertz Observatory for Gravitational Waves (NANOGrav; McLaughlin 2013) has recently found the first convincing evidence for a stochastic GWB signal, as detailed in Agazie et al. (2023b, hereafter NG15gwb). Analyzing 15-year of pulsar timing observations, NANOGrav has detected a rednoise process whose spectral properties are common among all pulsars and that is spatially correlated among pulsar pairs in a manner consistent with an isotropic GWB. In the following, we will refer to this observation as 'the NANOGrav signal,' 'the GWB signal,' or simply 'the signal,' keeping in mind the level of statistical significance at which the GW nature of the signal has been demonstrated in NG15gwb. While the GWB is primarily expected to arise from a population of inspiraling supermassive black hole binaries (SMBHBs; Rajagopal & Romani 1995; Jaffe & Backer 2003; Wyithe & Loeb 2003; Sesana et al. 2004; Burke-Spolaor et al. 2019), cosmological sources may also contribute to it.  \nThe SMBHB interpretation of the signal is considered in Agazie et al. (2023c, hereafter NG15smbh). In this paper, we analyze the NANOGrav 15-year data set (Agazie et al. 2023a, hereafter NG15) to investigate the possibility that the observed signal is cosmological in nature or that it arises from a combination of SMBHBs and a cosmological source. In particular, we consider phenomenological models of cosmic inflation, SIGWs, first-order phase transitions, cosmic  \nstrings (stable, metastable, and superstrings), and domain walls. We find that all of these models, except for stable cosmic strings of field theory origin, are consistent with the observed GWB signal. Many models provide in fact a better fit of the NANOGrav data than the baseline SMBHB model, which is reflected in the outcome of a comprehensive Bayesian model comparison analysis that we perform: several new-physics models result in Bayes factors between 10 and 100. We also consider composite models where the GWB spectrum receives contributions from new physics and SMBHBs. Comparing these composite models to the SMBHB reference model leads to comparable results, again with many Bayes factors falling into the range from 10 to 100. Cosmic superstrings, as predicted by string theory, are among the models that provide a good fit of the data, while stable cosmic strings of field theory origin only result in Bayes factors in the range from 0.1 to 1.  \nThe reason that some of the Bayes factors reach large values is that the SMBHB signal expected from the theoretical model used in this analysis agrees somewhat poorly (only at the level of 95% regions) with the observed data, leaving room for improvement by adding additional sources or better noise modeling. It is perhaps an intriguing idea that this disagreement may point to the presence of a cosmological source, but the present evidence is quite weak. We stress that Bayes factors for additional models beyond the SMBHB interpretation are highly dependent on the range of priors with which these models are introduced. Thus, one should not assign too much meaning to the exact numerical values of the Bayes factors reported in this work.  \nIn many models, there are ranges of parameter values that would produce signals in conflict with the NG15 data. In those cases, we show the excluded regions and give numerical upper limits for individual parameters. We do so in terms of a new statistical test, introducing what we call the K ratio. These parameter constraints are independent of the origin of the signal in the NG15 data and a testament to the constraining power of PTA data in the search for new physics. In our parameter plots, we label the K -ratio constraints by NG15, and where applicable, we compare them to other existing bounds. In many cases, the NG15 bounds are complementary to existing bounds, highlighting the fact that new-physics searches at the PTA frontier venture into previously unexplored regions of parameter space.  \nAside from cosmological GWBs, signals of new physics can appear in GW detectors in a deterministic manner. Although pulsar timing arrays (PTAs) are primarily used to search for a GWB, we can also leverage their remarkable sensitivity to search for these deterministic  \nsignals. Specifically, DM substructures within the Milky Way can produce a Doppler effect by accelerating the Earth or a pulsar (Seto & Cooray 2007), or a Shapiro delay of the photons' arrival times by perturbing the metric along the photon geodesic (Siegel et al. 2007). PTAs can also probe models of ultralight DM (ULDM), which can cause shifts in the observed pulse timing via metric fluctuations (Khmelnitsky & Rubakov 2014; Porayko & Postnov 2014) or via couplings between ULDM and SM particles (Graham et al. 2016; Kaplan et al. 2022). We search for both of these deterministic signals, and after finding no evidence for either of them, we derive new bounds on both these models.  \nThis paper is organized as follows. We describe the NG15 data set in Section 2 and our general analysis methods in Section 3. In Section 4, we discuss the GWB expected from SMBHBs. We present the analysis and results for new-physics models that generate a cosmological GWB in Section 5 and for models that produce deterministic signals in Section 6. We conclude in Section 7. Additionally, we include a list of parameters for each model, the prior ranges we use in our analysis, and the corresponding recovered posterior ranges in Appendix A. We present median GW spectra for all cosmological models based on our recovered posterior distributions in Appendix B, and we provide supplementary material for specific models in Appendix C .", '2. PTA DATA': "The NANOGrav 15-year (NG15) data set consists of observations of 68 millisecond pulsars made between 2004 July and 2020 August. This updated data set adds 21 pulsars and 3 yr of observations to the previous 12.5 yr data set (Alam et al. 2021). One pulsar, J0614-3329, was observed for less than 3 yr, which is why it is not included in our analysis. The remaining 67 pulsars were all observed for more than 3 yr with an approximate cadence of 1 month (with the exception of six pulsars that were observed weekly as part of a high-cadence campaign, which started in 2013 at the Green Bank Telescope and in 2015 at the Arecibo Observatory).  \nThe pulse times of arrival (TOAs) were generated from the raw data following the procedure discussed in Arzoumanian et al. (2015, 2018b) and Alam et al. (2021). The resulting cleaned TOAs were fit to a timing model that accounts for the pulsar's period and spin period derivative, sky location, proper motion, and parallax. For pulsars in a binary system, we included in the timing model five Keplerian binary parameters and an additional non-Keplerian parameter if they improved the fit as determined by an F -test. Pulse dispersion was  \nmodeled as a piecewise constant with the inclusion of DMX parameters (Arzoumanian et al. 2015; Jones et al. 2017). The timing model fits were performed using the TT(BIPM2019) timescale and the JPL Solar System Ephemeris model DE440 (Park et al. 2021). Additional detail about the data set and the processing of the TOAs can be found in NG15 and Agazie et al. (2023d, hereafter NG15detchar).", '3. DATA ANALYSIS METHODS': 'The statistical tools needed to describe noise sources, GWBs, and deterministic signals in pulsar timing data have already been extensively discussed in the literature (see, e.g., Arzoumanian et al. 2016, 2018a). In the following brief overview, we focus on the implementation of new-physics signals within this framework.', '3.1. Likelihood': "Our search for a new-physics signal utilizes the pulsars' timing residuals, δt . These timing residuals measure the discrepancy between the observed TOAs and the ones predicted by the pulsar timing model described in NG15 and briefly summarized in Section 2. There are three main contributions to these timing residuals: white noise, time-correlated stochastic processes (also known as red noise), and small errors in the fit to the timing-ephemeris parameters (Vallisneri et al. 2020). Specifically, we can model the timing residuals as  \nδt = n + F a + Mϵ . (1)  \nIn the remainder of this section, we will define and discuss each of these three terms and define the PTA likelihood.  \nThe first term on the right-hand side of Eq. (1), n , describes the white noise that is assumed to be left in each of the N TOA timing residuals after subtracting all known systematics. White noise is assumed to be a zero mean normal random variable, fully characterized by its covariance. For the receiver/back-end combination I , the white-noise covariance matrix reads  \n⟨ n i n j ⟩ = F 2 I [ σ 2 i S / N + Q 2 I ] δ ij + J 2 I U ij , (2)  \nwhere i and j index the TOAs, σ i S / N is the TOA uncertainty for the i th observation, F I is the Extra FACtor (EFAC) parameter, Q I is the Extra QUADrature (EQUAD) parameter, and J I is the ECORR parameter. ECORR is modeled using a block diagonal matrix, U , with values of 1 for TOAs from the same observing epoch and zeros for all other entries. Following the approach of previous works (Arzoumanian et al. 2016, 2018a), we fix all white-noise parameters to their values  \nat the maxima in the posterior probability distributions recovered from single pulsar noise studies in order to increase computational efficiency (NG15detchar).  \nTime-correlated stochastic processes, like pulsarintrinsic red noise and GWB signals, are modeled using a Fourier basis of frequencies i/T obs , where i indexes the harmonics of the basis and T obs is the timing baseline, extending from the first to the last recorded TOA in the full PTA data set. Since we are generally interested in processes that exhibit long-timescale correlations, the expansion is truncated after N f frequency bins. In this paper, we use N f = 30 for pulsar-intrinsic red noise and N f = 14 for GWBs. The latter choice stems from the observation that most of the evidence for a GWB comes from the first 14 frequency bins. More specifically, fitting a common-spectrum uncorrelated red-noise process with a broken power-law spectral shape to the NG15 data, the posterior distribution for the break frequency reaches it maximum around the 14th frequency bin (NG15gwb). This set of 2 N f sine-cosine pairs evaluated at the different observation times is contained in the Fourier design matrix, F . The Fourier coefficients of this expansion, a , are assumed to be normally distributed random variables with zero mean and covariance matrix, ⟨ aa T ⟩ = ϕ , given by  \n[ ϕ ] ( ak )( bj ) = δ ij (Γ ab Φ i + δ ab φ a,i ) (3)  \nwhere a and b index the pulsars, i and j index the frequency harmonics, and Γ ab is the GWB overlap reduction function, which describes average correlations between pulsars a and b as a function of their angular separation in the sky. For an isotropic and unpolarized GWB, Γ ab is given by the Hellings & Downs correlation (Hellings & Downs 1983), also known as 'quadrupolar' or 'HD' correlation.  \nThe first term on the right-hand side of Eq. (3) parameterizes the contribution to the timing residuals induced by a GWB in terms of the model-dependent coefficients Φ i . In this work, we consider two kinds of GWB sources: one of astrophysical origin, namely a population of inspiraling SMBHBs (discussed in section 4), and one of cosmological origin, induced by one of the exotic newphysics models under consideration (discussed in section 5). The last term in Eq. (3) models pulsar-intrinsic rednoise in terms of the coefficients φ a,i , where  \nφ a ( f ) = A 2 a 12 π 2 1 T obs ( f 1 yr -1 ) -γ a yr 3 (4)  \nand φ a,i = φ a ( i/T obs ) for all N f frequencies. The priors for the red noise parameters are reported in Table 2.  \nFinally, deviations from the initial best-fit values of the m timing-ephemeris parameters are accounted for by  \nthe term Mϵ . The design matrix , M , is an N TOA × m matrix containing the partial derivatives of the TOAs with respect to each timing-ephemeris parameter (evaluated at the initial best-fit value), and ϵ is a vector containing the linear offset from these best-fit parameters.  \nSince in this analysis we are not interested in the specific realization of the noise but only in its statistical properties, we can analytically marginalize over all the possible noise realizations (i.e., integrate over all the possible values of a and ϵ ). This leaves us with a marginalized likelihood that depends only on the (unknown) parameters describing the red-noise covariance matrix (i.e., A a , γ a , plus any other parameters describing Φ i ; van Haasteren & Levin 2012; Lentati et al. 2013):  \np ( δt | ϕ ) = exp ( -1 2 δt T C -1 δt ) √ det(2 π C ) , (5)  \nwhere C = N + TBT T . Here N is the covariance matrix of white noise, T = [ M , F ], and B = diag( ∞ , ϕ ), where ∞ is a diagonal matrix of infinities, which effectively means that we assume flat priors for the parameters in ϵ . Since in our calculations we always deal with the inverse of B , all these infinities reduce to zeros.  \nEq. (5) can be easily generalized to take into account deterministic signals (like the ones that will be discussed in Sections 6.1 and 6.2). In the presence of a deterministic signal, h ( θ ), which depends on a set of parameters θ , we just need to shift the residuals, δt → δt -h ( θ ).  \nFinally, we relate our characterization of the GWB given in Eq. (3) in terms of Φ i to the commonly adopted spectral representation in terms of the GWB energy density per logarithmic frequency interval, dρ GW /d ln f , as a fraction of the closure density, i.e., the total energy density of our universe, ρ c (Allen & Romano 1999)  \nΩ GW ( f ) ≡ 1 ρ c dρ GW ( f ) d ln f = 8 π 4 f 5 H 2 0 Φ( f ) ∆ f . (6)  \nHere H 0 is the present-day value of the Hubble rate, ∆ f = 1 /T obs is the separation between the N f frequency bins, and Φ( f ) determines the coefficients Φ i in Eq. (3), i.e., Φ i = Φ( i/T obs ). Note that Φ( f ) is identical to the timing residual power spectral density (PSD), S ( f ) = Φ( f ) / ∆ f , up to the constant factor of 1 / ∆ f . In the remainder of this paper, we will often work with h 2 Ω GW instead of Ω GW , where h is the dimensionless Hubble constant, H 0 = h × 100 kms -1 Mpc -1 , such that the explicit value of H 0 cancels in the product h 2 Ω GW .", '3.2. Bayesian analysis': "The goal of this work is to investigate a series of cosmological interpretations of the GWB signal in our  \ndata. Specifically, we would like to answer two questions. First, what is the region in the parameter space of the new-physics models that could produce the observed GWB? And second, is there any preference between the astrophysical and cosmological interpretations of the signal?  \nTo answer these questions, we make use of Bayesian inference. Bayesian inference is a statistical method in which Bayes' rule of conditional probabilities is used to update one's knowledge as observations are acquired. Given a model H , a set of parameters Θ, and data D , we can use Bayes' rule to write  \nP (Θ |D , H ) = P ( D| Θ , H ) P (Θ |H ) P ( D|H ) , (7)  \nwhere P (Θ |D , H ) is the posterior probability distribution for the model parameters, P ( D| Θ , H ) is the likelihood, P (Θ |H ) is the prior probability distribution, and  \nZ ≡ P ( D|H ) = ∫ dΘ P ( D| Θ , H ) P (Θ |H ) (8)  \nis the marginalized likelihood, or evidence. In the context of this work, H is the timing residual model given in Eq. (1), Θ contains the parameter describing the covariance matrix ϕ , and the data are the timing residuals δt . The likelihood function for our analysis is given by Eq. (5) and implemented using the ENTERPRISE (Ellis et al. 2019) and ENTEPRISE EXTENSIONS (Taylor et al. 2021) packages. Our prior choices are summarized in Tables 2 and 3.  \nThe posterior distribution on the left-hand side of Eq. (7) is the central result of the Bayesian analysis and contains all the information needed to answer our two original questions. Indeed, integrating over all the model parameters except one (two) allows us to derive marginalized distributions that can be used to obtain 1D (2D) credible intervals. At the same time, given two models H 0 and H 1 , we can perform model selection by calculating the Bayes factor defined as  \nB 10 ( D ) = Z 1 Z 0 = P ( D|H 1 ) P ( D|H 0 ) . (9)  \nThe numerical value of the Bayes factor for a given model comparison can then be interpreted as evidence against or in favor of model hypothesis H 1 according to the Jeffreys scale (Jeffreys 1961): B 10 < 1 means that H 1 is disfavored, while B 10 values in the ranges [10 0 . 0 , 10 0 . 5 ], [10 0 . 5 , 10 1 . 0 ], [10 1 . 0 , 10 1 . 5 ], [10 1 . 5 , 10 2 . 0 ], [10 2 . 0 , ∞ ) are interpreted as negligibly small, substantial, strong, very strong, and decisive evidence in favor of H 1 , respectively.  \nGiven the large number of parameters, the integration required to derive marginalized distributions and Bayes factors needs to be performed through Monte Carlo sampling. Specifically, we use the Markov Chain Monte Carlo (MCMC) tools implemented in the PTMCMCSampler package (Ellis & van Haasteren 2017) to sample from the posterior distributions. The marginalized posterior densities shown in our plots are then derived by applying kernel density estimates to the MCMC samples via the methods implemented in the GetDist package (Lewis 2019).  \nIn order to compute the Bayes factor between two models, we use product space methods (Carlin & Chib 1995; Godsill 2001; Hee et al. 2015), instead of calculating the evidence Z for each model separately. This procedure recasts model selection as a parameter estimation problem, introducing a model indexing variable that is sampled along with the parameters of the competing models and controls which model likelihood is active at each MCMC iteration. The ratio of samples spent in each bin of the model indexing variable returns the posterior odds ratio between models, which coincides with the Bayes factor for equal model priors, P ( H 1 ) = P ( H 0 ). The Monte Carlo sampling uncertainties associated with this derivation of the Bayes factors can be estimated through statistical bootstrapping (Efron & Tibshirani 1986). Bootstrapping creates new sets of Monte Carlo draws by resampling (with replacement) the original set of draws. These sets of draws act as independent realizations of the sampling procedure and allow us to obtain a distribution for the Bayes factors from which we derive point values and uncertainties on our Bayes factors corresponding to mean and standard deviation. Specifically, the central values and corresponding errors quoted in the following for the Bayes factors were derived by creating 5 × 10 4 realizations of our Monte Carlo draws.  \nFrom Eq. (8), it is evident that models' evidence and, therefore, Bayes factors depend on the prior choice. In our analysis, we will often restrict priors to the region of parameter space for which cosmological models produce an observable signal in the PTA frequency band. However, a more appropriate prior choice would cover the entire allowed region of parameter space. Nonetheless, when working with flat priors, it is easy to rescale the Bayes factors to account for wider prior ranges. Specifically, if the priors are extended to a region of parameter space for which the likelihood P ( D| Θ , H ) is approximately zero, the Bayes factors decrease by a factor proportional to the increase in prior volume.  \nFor each model H considered in our analysis, we use the reconstructed posterior distribution, P (Θ |D , H ), to  \nFigure 1. Comparison of the 68% and 95% probability regions for the amplitude and slope of a power-law fit to the observed GWB signal (green contours) and predicted for purely GW-driven SMBHB populations with circular orbits (blue contours; NG15smbh). The black dashed lines represent a 2D Gaussian fit of the blue contours. The vertical red line indicates γ = 13 / 3, the naive expectation for a GWB produced by a GW-driven SMBHB population (Phinney 2001).  \n<!-- image -->  \nidentify relevant parameter ranges and set upper limits. Specifically, we identify 68% (95%) Bayesian credible intervals (Bernardo & Smith 2000) by integrating the posterior over the regions of highest density until the integral covers 68% (95%) of the posterior probability. Moreover, we give upper limits above which the additional model is 'strongly disfavored' according to the Jeffreys scale (Jeffreys 1961). For instance, to place a bound on a single parameter θ , we first marginalize over all other model parameters and then determine the parameter value at which the likelihood ratio  \nK ( θ ) = P ( D| θ, H ) P ( D| θ 0 , H ) , (10)  \nhas dropped to K = 1 / 10 . Here θ 0 refers to the parameter limit in which the new-physics contribution to the total signal becomes negligible and P ( D| θ, H ) no longer depends on the exact value of θ . Graphing P ( D| θ, H ) as a function of θ , this parameter region appears as a plateau, with P ( D| θ 0 , H ) denoting the height of this plateau. Assuming a flat prior on θ , the ratio in Eq. (10) is identical to the corresponding ratio of marginalized posteriors. Furthermore, multiplying and dividing by the prior on θ ,  \nK ( θ ) = P ( θ |H ) P ( θ 0 |D , H ) P ( θ |D , H ) P ( θ |H ) . (11)  \nFigure 2. Bayes factors for the model comparisons between the new-physics interpretations of the signal considered in this work and the interpretation in terms of SMBHBs alone. Blue points are for the new physics alone, and red points are for the new physics in combination with the SMBHB signal. We also plot the error bars of all Bayes factors, which we obtain following the bootstrapping method outlined in Section 3.2. In most cases, however, these error bars are small and not visible.  \n<!-- image -->  \nThe first factor is the Savage-Dickey density ratio and can hence be identified as the Bayes factor B = P ( D|H ) /P ( D|H 0 ), where H 0 is the model that results from model H when omitting the signal contribution controlled by the parameter θ . The K ratio can thus be written as the product of the global Bayes factor and the local posterior-to-prior ratio for the parameter θ ,  \nK ( θ ) = B P ( θ |D , H ) P ( θ |H ) . (12)  \nOnce B is known, it is straightforward to evaluate Eq. (12) and determine the K -ratio bound on θ . Eq. (12) is useful for numerically evaluating K , as it automatically encodes the height of the plateau in the marginalized posterior, P ( θ 0 |D , H ) = P ( θ |H ) / B , which we would otherwise have to obtain from a fit to our MCMC data. However, we stress that K is defined as a likelihood ratio, which renders it immune to prior effects (prior choice, range, etc.; Azzalini 1996). For more than one parameter dimension, we proceed analogously and derive bounds based on the criterion K (Θ) > 1 / 10 .  \nAll Bayesian inference analyses discussed in this work were implemented into ENTERPRISE via a newly developed wrapper that we call PTArcade (Mitridate et al. 2023a,b). This wrapper is intended to allow easy implementation of new-physics searches in PTA data. We make this wrapper publicly available at https://doi.org/ 10.5281/zenodo.7876429. Similarly, all MCMC chains analyzed in this work can be downloaded at https: //zenodo.org/record/8010909.", '4. GWB SIGNAL FROM SMBHBs': "Most galaxies are expected to host a supermassive black hole (SMBH) at their center (Kormendy & Ho 2013; Akiyama et al. 2019). During the hierarchical merging of galaxies taking place in the course of structure formation (White & Rees 1978), these black holes are expected to sink to the center of the merger remnants, eventually forming binary systems (Begelman et al. 1980). The gravitational radiation emitted by this population of inspiraling SMBHBs forms a GWB in the PTA band (Rajagopal & Romani 1995; Jaffe & Backer 2003; Wyithe & Loeb 2003) and is a natural candidate for the source of the signal observed in our data.  \nThe shape and normalization of this GWB depend on the properties of the SMBHB population and on its dynamical evolution (Enoki & Nagashima 2007; Sesana et al. 2008; Kocsis & Sesana 2011; Kelley et al. 2017). As discussed in NG15smbh, the normalization is primarily controlled by the typical masses and abundance of SMBHBs, while the shape of the spectrum is determined by subparsec-scale binary evolution, which is currently unconstrained by observations. For a population of binaries whose orbital evolution is driven purely by GW emission, the resulting timing residual PSD is a power law with a spectral index (defined below in Eq. (13)) of -γ BHB = -13 / 3 (Phinney 2001), produced by the increasing rate of inspiral and decreasing number of binaries emitting over each frequency interval. However, as GW emission alone is typically insufficient  \nFigure 3. Median GWB spectra produced by a subset of the new-physics models, which we construct by mapping our model parameter posterior distributions to h 2 Ω GW distributions at every frequency f (see Appendix B for more details and Figs. 19 and 20 for the models not included here). We also show the periodogram for an HD-correlated free spectral process (gray violins) and the GWB spectrum produced by an astrophysical population of inspiraling SMBHBs with the parameters A BHB and γ BHB fixed at the central values µ BHB of the 2D Gaussian prior distribution specified in Eq. (A1) (black dashed line).  \n<!-- image -->  \n4  \nto merge SMBHBs within a Hubble time, the number of binaries emitting in the PTA band depends on interactions between binaries and their local galactic environment to extract orbital energy and drive systems toward merger (Begelman et al. 1980). If these environmental effects extend into the PTA band, or if binary orbits are substantially eccentric, then the GWB spectrum can flatten at low frequencies (typically expected at f ≪ 1 yr -1 ; Kocsis & Sesana 2011). At high frequencies, once the expected number of binaries dominating the GWB approaches unity, the spectrum steepens below 13 / 3 (typically expected at f ≫ 1 yr -1 ; Sesana et al. 2008).  \nUnfortunately, current observations and numerical simulations provide only weak constraints on the spectral amplitude or the specific locations and strengths of power-law deviations. Despite these uncertainties, the sensitivity range of PTAs is sufficiently narrowband that it is reasonable, to first approximation, to model the sig-  \nnal by a power law in this frequency range:  \nΦ BHB ( f ) = A 2 BHB 12 π 2 1 T obs ( f yr -1 ) -γ BHB yr 3 , (13)  \nwhere Φ BHB / ∆ f is the timing residual PSD (see Eq. (6)).  \nFollowing Middleton et al. (2021), we can gain some insight into the allowed range of values for the amplitude, A BHB , and slope, γ BHB , of this power law by simulating a large number of SMBHB populations covering the entire range of allowed astrophysical parameters. Specifically, we consider the SMBHB populations contained in the GWOnly-Ext library generated as part of the NG15smbh analysis (and discussed in additional detail there). This library was constructed with the holodeck package (Kelley et al. 2023) using semianalytic models of SMBHB mergers. These models use simple, parameterized forms of galaxy stellar mass functions, pair fractions, merger rates, and SMBH-mass versus galaxy-mass relations to produce binary populations and derived GWB spectra. While some param-  \nFigure 4. Same as Fig. 3 but for a different selection of models and showing a larger frequency range. The solid lines represent median GWB spectra for a subset of new-physics models (see Appendix B for more details); the gray violins correspond to the posteriors of an HD-correlated free spectral reconstruction of the NANOGrav signal; and the shaded regions indicate the power-law-integrated sensitivity (Thrane & Romano 2013) of various existing and planned GW interferometer experiments: LISA (Amaro-Seoane et al. 2017), DECIGO (Kawamura et al. 2011), BBO (Crowder & Cornish 2005), Einstein Telescope (ET; Punturo et al. 2010), Cosmic Explorer (CE; Reitze et al. 2019), the HLVK detector network (consisting of aLIGO in Hanford and Livingston (Aasi et al. 2015), aVirgo (Acernese et al. 2015), and KAGRA (Akutsu et al. 2019)) at design sensitivity, and the HLV detector network during the third observing run (O3). All sensitivity curves are normalized to a signal-to-noise ratio of unity and, for planned experiments, an observing time of one year. For the HLV detector network, we use the O3 observing time. Different signal-to-noise thresholds ρ thr and observing times t obs can be easily implemented by rescaling the sensitivity curves by a factor of ρ thr / √ t obs . More details on the construction of the sensitivity curves can be found in Schmitz (2021). We emphasize that models whose median GWB spectrum exceeds the sensitivity of existing experiments are not automatically ruled out. This applies, e.g., to cosmic superstrings ( super ) and the O3 sensitivity of the HLV detector network. Typically, no single GWB spectrum in a given model will coincide with the median GWB spectrum, which is constructed from distributions of h 2 Ω GW values at any given frequency. Therefore, if the median GWB spectrum is in conflict with existing bounds, typically only some regions in the model parameter space will be ruled out, while others remain viable (see, e.g., Fig. 11 for the super model). Finally, note that any primordial GWB signal is subject to the upper limit on the amount of dark radiation in Eq. (23), which requires the total integrated GW energy density to remain smaller than O (10 -(5 ··· 6) ) (see Section 5.1).  \n<!-- image -->  \nters in these models are fairly well known (e.g., concerning the galaxy stellar mass function), others are almost entirely unconstrained-particularly those governing the dynamical evolution of SMBHBs on subparsec scales (Begelman et al. 1980). The GWOnly-Ext library assumes purely GW-driven binary evolution and uses relatively narrow distributions of model parameters based on literature constraints from galaxy-merger  \nobservations (e.g., Tomczak et al. 2014) in addition to more detailed numerical studies of SMBHB evolution (e.g., Sesana 2013).  \nFor each population contained in the GWOnly-Ext library, we perform a power-law fit of the corresponding GWB spectrum across the first 14 frequency bins that we use in our analysis. The distribution for A BHB and γ BHB obtained in this way is reported in Fig. 1  \n(blue contours) and compared to the results of a simple power-law fit to the GWB signal in the NG15 data set (green contours). The 95% regions of the two distributions barely overlap, signaling a mild tension between the astrophysical prediction and the reconstructed spectral shape of the GWB. In view of this observation, we stress again that while these simulated populations are consistent with systematic investigations of the GWB spectrum (e.g., Sesana 2013), they assume circular orbits and GW-only driven evolution. Adopting models that include either significant coupling between binaries and their local environments or very high eccentricities could serve to flatten the spectral shape and lead to SMBHB signals that better align with the observed data (see NG15smbh for an extended discussion). Neither of these effects, however, is expected to significantly impact the amplitudes of the predicted spectra that, for expected values of astrophysical parameters, remain in mild tension with observed data. As discussed in NG15smbh, in order to reproduce the observed amplitude, SMBHB models require one or more of the astrophysical parameters describing the binaries' population to differ from expected values. For the present analysis, the spectra derived from the GWOnly-Ext library thus represent a convenient benchmark that is simple, well defined, and easy to use. By using theory-motivated priors, our reference model constitutes an important step toward a more realistic modeling of the GWB spectrum from inspiraling SMBHBs that goes beyond a powerlaw parameterization with spectral index γ BHB = 13 / 3, which has been the standard reference model in much of the PTA literature over the past decades.  \nThe black dashed contours in Fig. 1 show the results of a 2D Gaussian fit to the distribution of A BHB and γ BHB values derived from the simulated SMBHB populations (see Eq. (A1) in Appendix A for the parameters of this Gaussian distribution). This fitted distribution is what we adopt as a prior distribution for A BHB and γ BHB in all parts of the analysis described in this paper.", '5. GWB SIGNALS FROM NEW PHYSICS': 'In this section, we discuss the GWB produced by various new-physics models and investigate each model alone and in combination with the SMBHB signal as a possible explanation of the observed GWB signal. For each model, we give a brief review of the mechanism behind the GWB production and discuss the parametrization of its signal prediction. We report the reconstructed posterior distributions of the model parameters and compute the Bayes factors against the baseline SMBHBinterpretation. In Fig. 2, we show a summary of these Bayes factors; in Fig. 3, we present median recon-  \nstructed GWB spectra in the PTA band for a number of select new-physics models; and in Fig. 4, we show similar median reconstructed GWB spectra in the broader landscape of present and future GW experiments.  \nAs discussed in Section 4 and in more detail in NG15smbh, there is a mild tension between the NG15 data and the predictions of SMBHB models. The models generally prefer a weaker and less blue-tilted h 2 Ω GW spectrum than the data. This discrepancy presents an opportunity for new-physics models to fit the data better than the conventional SMBHB signal. Eventually, this tension may grow to the point of giving strong evidence for new physics, or it may be resolved with better modeling and more data. Specifically, models of SMBHB evolution with a significant coupling between binaries and their local environment could lead to a signal that better aligns with the data and reduce the evidence for new physics. For all these reasons, we caution against over-interpreting the observed evidence in favor of some of the new-physics models discussed in the following sections.', 'MODEL DESCRIPTION': "In the ΛCDM model, the structures we observe in the universe are seeded by primordial curvature fluctuations generated during inflation and then imprinted onto the DM density field. CMB observations indicate that these fluctuations have a nearly scale-invariant power spectrum on large scales (i.e., for comoving wavenumbers k ≃ Mpc -1 ). However, on smaller scales, various theories of DM leave unique fingerprints on primordial perturbations or their evolution, resulting in different predictions for the amount of subgalactic DM substruc-  \nFigure 15. The red solid (dashed) lines show our constraints at the 95% credible level for the two vector ULDM models considered in this analysis, for both the correlated (solid line) and uncorrelated (dashed line) regimes. The lower gray shaded regions correspond to regions of parameter space where the signal amplitude is smaller than the purely gravitational signal. The constraints from tests of the equivalence principle ('MICROSCOPE') are shown in teal (Berg'e et al. 2018), and the constraints previously set by the PPTA Collaboration are reported in yellow (Xue et al. 2022).  \n<!-- image -->  \n<!-- image -->  \ng  \ntures. Consequently, measuring local DM substructures  \ncould be crucial in determining the correct model of DM. PBHs are perhaps the simplest example of such smallscale DM substructures. They can be formed in inflationary theories that create large density fluctuations on small scales (like the ones described in Section 5.2). Several studies investigated the possibility of identifying a galactic PBH population by analyzing the Doppler and Shapiro signals they can leave in PTA data (Seto & Cooray 2007; Siegel et al. 2007; Dror et al. 2019; Ramani et al. 2020; Lee et al. 2021a,b). In this analysis, we will closely follow the method outlined by Lee et al. (2021b) to constrain the local PBH abundance. 10  \nThe Doppler signal results from the apparent shift in the pulsar spin frequency, generated by the acceleration induced by the gravitational pull of a passing PBH. According to the timescale of the transit event, τ , the signal can be further classified as dynamic ( static ) if τ is much smaller (larger) than the observation time, T obs . In the static regime, the leading-order term of the Doppler signal that is not degenerate with the timing model is given  \nby (Ramani et al. 2020; Lee et al. 2021a,b)  \nh D, sta ( t ) = A D, sta yr 2 t 3 , (67)  \nwhere A D, sta is a dimensionless parameter that can be related to the kinematic parameters of the transiting event (see Appendix C.6 for more details). In the dynamic limit, and assuming that the signal is dominated by the closest transiting PBH, we get  \nh D, dyn ( t ) = A D, dyn ( t -t 0 ) Θ( t -t 0 ) , (68)  \nwhere Θ is the Heaviside step function, A D, dyn is a dimensionless amplitude that can be related to kinematic parameters of the transiting event, and t 0 is the time of closest approach (see Appendix C.6 for more details).  \nThe Shapiro signal refers to shifts in the TOAs caused by metric perturbations along the photon geodesic produced by PBHs transiting along the observer's line of sight. In the static limit, and after subtracting away degeneracies with timing model parameters, the leadingorder term of this signal can be parameterized as  \nh S, sta ( t ) = A S, sta yr 2 t 3 , (69)  \nwhere, as for the Doppler case, A S, sta is a dimensionless parameter that can be related to the kinematic parameters of the transiting event (see Appendix C.6 for  \nFigure 16. The black solid (dashed) lines show the posterior distributions p (log 10 A sta | δt ) ( p (log 10 A dyn | δt )) for a representative pulsar (J1909-3744). The filled distributions show the conditional probability distributions p (log 10 A | f PBH ) for the same pulsar and different values of f PBH . In this plot, M = 10 -6 (10 -10 ) M ⊙ for the Doppler static (dynamic) signal, and M = 10 M ⊙ for the Shapiro static signal.  \n<!-- image -->  \nmore details). In the dynamic limit, there is no simple  \nparametrization of the Shapiro signal; therefore, we do not search for this signal.  \nAssuming a monochromatic PBH population, our goal is to derive a posterior distribution for the PBH mass fraction, f PBH ≡ Ω PBH / Ω DM , as a function of the PBH mass, M : p ( f PBH | δt , M ). We do this as follows:  \n- · For each given value of f PBH and M , we use the Monte Carlo code developed by Lee et al. (2021a) to generate a PBH population surrounding each of the pulsars in our array. From this distribution, we derive the amplitude of the static Doppler and Shapiro signals generated by the entire PBH population and the amplitude for the dynamic Doppler signal generated by the closest transiting PBH. Finally, we repeat this procedure for 2 . 5 × 10 3 realizations to obtain the conditional probability distributions p (log 10 A I | f PBH ), where I indexes pulsars in the array and A refers to any of the PBH signal amplitudes introduced in Eqs. (67), (68), and (69). In Fig. 16 we report some of the distributions derived in this way. 11\n- · One at a time, we include the PBH signals given in Eqs. (67), (68), and (69) in the timing model, and we analyze our data to derive the posterior distributions for the various PBH signal amplitudes, p (log 10 A | δt ). Since the PBH signal in different pulsars is assumed to be independent, these distributions can be factorized as  \np (log 10 A | δt ) = N P ∏ I =1 p (log 10 A I | δt ) . (70)  \nSome of the p (log 10 A I | δt ) are reported in Fig. 16.  \n- · Finally, we can write  \np ( f PBH | δt ) = N P ∏ I =1 ∫ d log 10 A I p ( f PBH | log 10 A I ) p (log 10 A I | δt ) ∝ N P ∏ I =1 ∫ d log 10 A I p (log 10 A I | f PBH ) p (log 10 A I | δt ) (71)  \nwhere, in the second step, we used Bayes theorem and assumed uniform priors on log 10 A I and f PBH . More details on each of these three steps can be found in Appendix C.6 or in Lee et al. (2021b).  \nDM substructures can also possess macroscopic charges and interact with baryonic matter via long-range Yukawa interactions. These interactions can be modeled  \nby a potential of the form  \nV fifth ( r ) = -˜ α GMM P r e -r/λ , (72)  \nwhere M and M P are the masses of the DM and pulsar, respectively, λ is the range of the interaction, and ˜ α is the effective DM-barion coupling, normalized against the gravitational coupling (also known as the Yukawa parameter). Here the DM can be either a particle or a macroscopic object such as a nugget of asymmetric DM (Detmold et al. (2014); Wise & Zhang (2014); Hardy et al. (2015); Krnjaic & Sigurdson (2015); Gresham et al. (2017, 2018)). These Yukawa interactions can arise from an effective Lagrangian of the form L ⊃  \ng X ϕ ¯ XX + g n ϕ ¯ nn , where X and n are the effective DM and neutron fields, and ϕ is a massive (but potentially light) scalar or vector field. The effective coupling is related to the coupling constants by ˜ α ≈ g n g X 4 πGm X m n , where m n is the neutron mass. These interactions are constrained to be weaker than gravity for the mass range M ≲ 100 GeV by the CMB, Lymanα forest (Xu et al. 2018), and direct detection experiments such as X-ray Quantum Calorimeter (XQC) (Mahdawi & Farrar 2018) and Cryogenic Rare Event Search with Superconducting Thermometers (CRESST) (Angloher et al. 2017) (for a review on these constraints see Xu & Farrar (2020)). However, stronger-than-gravity fifth forces are allowed if M ≫ 100 GeV, even when X accounts for the entirety of the DM population.  \nIf present, these Yukawa interactions will contribute to the pulsar's acceleration induced by a transiting DM substructure and contribute to the Doppler signal discussed before (the expression for the Yukawa contribution to the Doppler signal can be found in Appendix C.6). Therefore, as shown by Gresham et al. (2023), following a procedure similar to the one used to constrain the abundance of PBHs, we can constrain the value of the Yukawa parameter, ˜ α . Specifically, for each given value of ˜ α and M , we use the Monte Carlo code developed by Lee et al. (2021a) to generate a population of DM substructure surrounding each of the pulsars in our array. From this distribution, we derive the amplitude of the static Doppler signal generated by the closest transiting substructure by considering the acceleration induced by both the gravitational and Yukawa interaction. By repeating this procedure for multiple populations of DM substructure, we derive the distribution p (log 10 A I | ˜ α ). By plugging this quantity into an expression similar to the one given in Eq. (71), we can derive p (˜ α | δt ) and use this quantity to constrain ˜ α .", 'RESULTS AND DISCUSSION': "We start by searching for PBH signals on top of a GWBthat we model as a power law with amplitude and spectral index allowed to float within the following prior ranges: log 10 A GWB ∈ [ -11 , -18] and γ GWB ∈ [0 , 7]. We find no statistically significant evidence for any of the PBH signals described in the previous section. Therefore, we proceed to set constraints on the local PBH abundance. The prior distributions used for the PBH signal parameters are reported in Table 3.  \nThe 95% upper limits on f PBH derived from the static Doppler and Shapiro signals are reported in Fig. 17. The dynamic Doppler signal is too weak to produce any detectable signal for any of the f PBH values considered. These are the first constraints on f PBH derived using  \nFigure 17. Constraints at the 95% credible level on the local PBH abundance derived from the search for static Doppler (red shaded region) and static Shapiro signals (blue shaded region). The solid lines interpolate between the PBH masses simulated in this work, while the red dashed line shows an extrapolation of the constraints to higher masses.  \n<!-- image -->  \n/circledot  \nreal PTA data. As expected, our constraints are much weaker compared to the projections that were derived by Lee et al. (2021b) using mock data and including only white noise. Indeed, as already discussed by Lee et al. (2021b), the presence of a common red-noise process significantly reduces the sensitivity to PBH signals.  \nFinally, in Fig. 18 we show the constraints on ˜ α set by NG15 data. These constraints are compared with several other constraints that can be placed on ˜ α . Specifically, in teal we show weak equivalence principle (WEP) constraints (Wagner et al. 2012; Shao et al. 2018; Sun et al. 2019) (properly rescaled to take into account the finite range of the interaction (Gresham et al. 2023)) derived by considering differential acceleration of baryonic test bodies toward the galactic center. In blue we report constraints from neutron star (NS) heating (assuming additional short-range DM-baryon interaction) induced by DM capture (Gresham et al. 2023), derived from the coldest known NS to date - PSR J2144-3933 (Guillot et al. 2019). And in gray we report the indirect constraints that can be derived by combining the fifthforce constraints on baryon-baryon interactions (Berg'e et al. 2018; Fayet 2019), and Bullet Cluster constraints on DM-DM interactions (Spergel & Steinhardt 2000; Kahlhoefer et al. 2014) (see Coskuner et al. (2019); Gresham et al. (2023)). We find that the NG15 constraints  \nFigure 18. The 95% credible level for the fifth-force strength ˜ α derived from the NG15 data (red lines) is compared with constraints from NS kinetic heating (blue lines), equivalence principle constraints (green lines), and Bullet Cluster + equivalence principle constraints (gray line). Solid (dashed) lines are deriving assuming λ = 1pc ( λ = 10 -1 pc), while dashed-dotted (dotted) lines are derived assuming λ = 10 -2 pc ( λ = 10 -3 pc).  \n<!-- image -->  \n/circledot  \nare competitive with WEP and NS kinetic heating, especially in the M ≳ 10 -12 M ⊙ regime. However, the combined constraint from the Bullet Cluster and MICROSCOPE dominates over all other constraints in the entire mass range of interest. We note, however, that the combined constraint is an indirect probe of the Yukawa parameter and depends on the specific form of the light mediator coupling to DM and baryons. Specifically, in deriving the combined constraints, we have assumed that the Yukawa DM-baryon interaction arises from a Lagrangian of the form L ⊃ g X ϕ ¯ XX + g n ϕ ¯ nn . In addition, if only a subcomponent of DM interacts through the long-range fifth force, then the Bullet Cluster constraint will quickly lift, while the other constraints only deteriorate linearly with the DM fraction.", '5.2. Scalar-induced gravitational waves MODEL DESCRIPTION': "The amplitude of the primordial scalar power spectrum is well measured by CMB observations, A s ≃ 2 . 10 × 10 -9 at the CMB pivot scale k CMB = 0 . 05 Mpc -1 (Aghanim et al. 2020). If we naively extrapolate this value down to smaller scales, assuming a fixed and slightly red-tilted h 2 Ω GW spectrum with index n s ∼ 0 . 96, we are led to conclude that there must be increasingly less power in scalar perturbations on shorter scales. This conclusion can, however, be easily avoided in models that deviate from the standard picture of single-field slow-roll inflation giving rise to a nearly scale-invariant spectrum of scalar perturbations. A prominent example, among many other mechanisms, consists in a stage of inflation close to an inflection point in the scalar potential, which readily amplifies the scalar perturbations leaving the horizon (see, e.g., GarciaBellido & Ruiz Morales (2017); Ezquiaga et al. (2018); Ballesteros & Taoso (2018)). An enhanced scalar power spectrum at small scales is, therefore, a viable possibility. Moreover, it promises a rich phenomenology with regard to the production of GWs and potentially the origin of primordial black holes (PBHs) (Carr et al. 2016; Garcia-Bellido et al. 2016; Inomata et al. 2017a; Inomata & Nakama 2019; Wang et al. 2019; Escriv'a et al. 2022b). The possibility of having PBH formation in models of single-field inflation is the subject of ongoing debate (Kristiano & Yokoyama 2022; Riotto 2023a; Choudhury et al. 2023a,b; Kristiano & Yokoyama 2023; Riotto 2023b; Choudhury et al. 2023c; Firouzjahi & Riotto 2023). Below, we comment on the implications of this debate for our PBH-related parameter bounds.  \nIn cosmological perturbation theory, scalar and tensor perturbations evolve independently at linear order. Starting at second order, however, they are coupled, which means that large first-order scalar perturbations can act as a source of second-order tensor perturbations. We refer to these tensor perturbations, which are produced as soon as the enhanced scalar perturbations reenter the Hubble horizon after inflation, as scalar-induced gravitational waves (SIGWs) (Matarrese et al. 1993, 1994, 1998; Mollerach et al. 2004; Ananda et al. 2007; Baumann et al. 2007) (see Dom'enech (2021) for a review and more details on the history of SIGWs). At the same time, large overdensities in the tail of the distribution of scalar perturbations can collapse into PBHs upon horizon reentry. This PBH production mechanism thus results in a PBH population whose properties are strongly  \ncorrelated with the spectral shape of the SIGW signal, as both phenomena are controlled by the scalar spectrum generated during inflation (Yuan & Huang 2021a). For earlier works on the PBH/SIGW interpretation of the signal in recent PTA data sets, see Vaskonen & Veermae (2021); De Luca et al. (2021); Kohri & Terada (2021). For earlier Bayesian searches for an SIGW signal in PTA data sets, see Chen et al. (2020); Bian et al. (2021); Zhao & Wang (2022); Yi & Fei (2023); Dandoy et al. (2023), and for a search in LVK data, see Romero-Rodriguez et al. (2022).  \nIn our analysis, we consider SIGWs in the PTA band and model the associated GW spectrum as follows:  \nΩ ind GW ( f ) = Ω r ( g ∗ ( f ) g 0 ∗ )( g 0 ∗ ,s g ∗ ,s ( f ) ) 4 / 3 ¯ Ω ind GW ( f ) , (26)  \nwhere the first three factors account for the cosmological redshift as in Eq. (17), and where the last factor denotes the GW spectrum at the time of production, which we assume to be during the radiation-dominated era,  \n¯ Ω ind GW ( f ) = ∫ ∞ 0 d v ∫ 1+ v | 1 -v | d u K ( u, v ) P R ( uk ) P R ( vk ) . (27)  \nThe present-day GW frequency is related to the comoving wavenumber by f = k/ (2 πa 0 ), P R denotes the primordial spectrum of the comoving curvature perturbation R , and the integration kernel K is given by (Espinosa et al. 2018; Kohri & Terada 2018; Pi & Sasaki 2020; Gong 2022)  \nK ( u, v ) = 3 ( 4 v 2 -(1 + v 2 -u 2 ) 2 ) 2 ( u 2 + v 2 -3 ) 4 1024 u 8 v 8 [ ( ln ∣ ∣ ∣ ∣ 3 -( u + v ) 2 3 -( u -v ) 2 ∣ ∣ ∣ ∣ -4 uv u 2 + v 2 -3 ) 2 + π 2 Θ( u + v -√ 3) ] . (28)  \nThe expression in Eq. (27) illustrates the dependence of the SIGW signal on the scalar input spectrum; in particular, it shows that Ω ind GW scales as Ω ind GW ∝ P 2 R . We stress that the expression in Eq. (27) assumes Gaussian perturbations, whose statistics are fully described by the power spectrum P R . We do not consider any non-Gaussian contributions to the primordial scalar power spectrum in our analysis. The impact of non-Gaussianities on the SIGW signal was studied in Nakama et al. (2017); Cai et al. (2019); Unal (2019); Atal et al. (2019); Yuan & Huang (2021b); Atal & Dom'enech (2021); Adshead et al. (2021); Ferrante et al. (2023).  \nTo remain as model-independent as possible, we refrain from choosing a particular inflation model capable of generating an enhanced spectrum P R . Instead, we ignore the microphysics of inflation and work with three characteristic templates for P R that reflect the range of possibilities that one may expect in realistic models. Specifically, we consider the following templates:  \nsigw-delta : Sharp spectral feature in P R modeled by a Dirac delta function in logarithmic k space,  \nP R ( k ) = Aδ (ln k -ln k ∗ ) . (29)  \nsigw-gauss : Broad spectral feature in P R modeled by a Gaussian peak in logarithmic k space,  \nP R ( k ) = A √ 2 π ∆ exp [ -1 2 ( ln k -ln k ∗ ∆ ) 2 ] . (30)  \nsigw-box : Flat and continuous spectral feature in P R modeled by a box function in logarithmic k space,  \nP R ( k ) = A Θ(ln k max -ln k ) Θ (ln k -ln k min ) . (31)  \nNote that the Gaussian power spectrum in logarithmic k space that we assume in the sigw-gauss model corresponds to a lognormal power spectrum in linear k space. As evident from the above expressions, sigw-delta represents a two-parameter model, while sigw-gauss and sigw-box are three-parameter models. Our prior choices for the respective parameters are listed in Table 3 in Appendix A, where we use again f = k/ (2 πa 0 ) to convert from wavenumber to frequency. For a given set of parameter values, we are then able to use the scalar power spectrum in Eq. (29), Eq. (30), or Eq. (31) to evaluate the integrals in Eq. (27) and compute the GW spectrum. For sigw-gauss and sigw-box , the integration needs to be carried out numerically; for sigwdelta , we can resort to the exact analytical expression provided in Wang et al. (2019); Yuan & Huang (2021a).", '5.3. Cosmological phase transitions MODEL DESCRIPTION': 'In the cosmological context, first-order phase transitions take place when a potential barrier separates the true and false minima of scalar field potential. 2 The phase transition is triggered by quantum or thermal fluctuations that cause the scalar field to tunnel through or fluctuate over the barrier, which results in the nucleation of bubbles within which the scalar field is in the true vacuum configuration. If sufficiently large, these bubbles expand in the surrounding plasma where the scalar field still resides in the false vacuum. The expansion and collision of these bubbles (Kosowsky et al. 1992a,b; Kosowsky & Turner 1993; Kamionkowski et al. 1994; Caprini et al. 2008; Huber & Konstandin 2008), together with sound waves generated in the plasma (Hindmarsh et al. 2014; Giblin & Mertens 2013, 2014; Hindmarsh et al. 2015), can source a primordial GWB (see Witten (1984); Hogan (1986) for seminal work). 3 For earlier Bayesian searches for a phase transition signal in PTA data, see Arzoumanian et al. (2021); Xue et al. (2021).  \nGenerally, the GWB produced during the phase transition is a superposition of the bubble and sound-wave contributions. However, if the bubble walls interact with the surrounding plasma, most of the energy released in the phase transition is expected to be converted to plasma motion, causing the sound-wave contribution to dominate the resulting GWB. An exception to this scenario is provided by models in which there are no (or only very feeble) interactions between the bubble walls and the plasma, or by models where the energy released in the phase transition is large enough that the friction exerted by the plasma is not enough to stop the walls from accelerating ( runaway scenario ). However, determining whether or not the runaway regime is reached is either model dependent or affected by large theoretical uncertainties (see, e.g., Li et al. (2023); Ai et al. (2023); Krajewski et al. (2023) for recent work on this topic). Therefore, in this work, we perform two separate analyses: a sound-wave-only analysis ( pt-sound ), where we assume that the runaway regime is not reached  \nand sound waves dominate the GW spectrum, and a bubble-collisions-only analysis ( pt-bubble ), where we assume that the runaway regime is reached and bubble collisions dominate the GW spectrum.  \nWe parameterize the GWB produced by both sound waves and bubble collisions in a model-independent way in terms of the following phase transition parameters:  \n- · T ∗ , the percolation temperature, i.e., the temperature of the universe when ∼ 34% of its volume has been converted to the true vacuum (Ellis et al. 2019a). For weak transitions, this temperature coincides with the temperature at the time of bubble nucleation, T n ∼ T ∗ . Conversely, for supercooled transitions, we typically have T n ≪ T ∗ . Barring the case of extremely strong transitions, α ∗ ≫ 1 (see below), which we do not consider in this work, T ∗ also determines the reheating temperature after percolation, T rh ∼ T ∗ (Ellis et al. 2019a).\n- · α ∗ , the strength of the transition, i.e., the ratio of the change in the trace of the energy-momentum tensor across the transition and the radiation energy density at percolation (Caprini et al. 2020; Ellis et al. 2019b).\n- · H ∗ R ∗ = R ∗ /H -1 ∗ , the average bubble separation in units of the Hubble radius at percolation, H -1 ∗ . For relativistic bubble velocities, which is what we consider in the following, R ∗ is related to the bubble nucleation rate, β , by the relation H ∗ R ∗ = (8 π ) 1 / 3 H ∗ /β .  \nIn addition to the parameters T ∗ , α ∗ , and H ∗ R ∗ , the GWB produced by a phase transition also depends on the velocity of the expanding bubble walls, v w . However, deriving the precise value of this quantity is an open theoretical problem, which depends on model-dependent quantities, such as the strength of the interactions between the bubble walls and the SM plasma. Therefore, in our analysis, we fix the bubble velocity to unity (i.e., the speed of light in natural units). This assumption is well justified for strong phase transitions (Bodeker & Moore 2017), which, realistically, are the only ones that could lead to a detectable signal in our current data. In particular, we fix v w = 1 for both phase transition scenarios that we are interested in, pt-sound and ptbubble . In the latter case, v w → 1 is automatically implied by the runaway behavior of the phase transition; in the former case, one actually expects a subluminal terminal velocity, v w < 1. In this sense, our decision to fix v w = 1 amounts to the optimistic assumption that this terminal velocity is numerically close to v w = 1. A similar approach is followed by the authors of the LISA review paper Caprini et al. (2020) who work with v w = 0 . 95 throughout most of their analysis in the absence of more detailed microphysical calculations. Finally, we point out that the parametrization of the GWB  \nsignal in terms of H ∗ R ∗ = (8 π ) 1 / 3 v w H ∗ /β already absorbs a large part of the dependence on the bubble wall velocity. The remaining v w dependence is mostly contained in the efficiency factor κ s (see below). However, in the regime of large α ∗ values, α ∗ ∼ 0 . 3 · · · 10, which turn out to be preferred by the NG15 data (see Fig. 8), this dependence is rather weak (see Fig. 13 in Espinosa et al. (2010)), which justifies again to keeping v w fixed.  \nThe GWB spectrum sourced by bubbles and sound waves can be written in terms of these parameters as  \nΩ b ( f ) = D ˜ Ω b ( α ∗ 1 + α ∗ ) 2 ( H ∗ R ∗ ) 2 S ( f/f b ) (33)  \nΩ s ( f ) = D ˜ Ω s Υ( τ sw ) ( κ s α ∗ 1 + α ∗ ) 2 ( H ∗ R ∗ ) S ( f/f s ) . (34)  \nHere ˜ Ω b = 0 . 0049 (Jinno & Takimoto 2017) and ˜ Ω s = 0 . 036 (Hindmarsh et al. 2017); the efficiency factor κ s = α ∗ / (0 . 73 + 0 . 083 √ α ∗ + α ∗ ) (Espinosa et al. 2010) gives the fraction of the released energy that is transferred to plasma motion in the form of sound waves, and D accounts for the redshift of the GW energy density,  \nD = π 2 90 T 4 0 M 2 Pl H 2 0 g ∗ ( g eq ∗ ,s g ∗ ,s ) 4 / 3 ≃ 1 . 67 × 10 -5 . (35)  \nWe recall that T 0 and H 0 denote the photon temperature and Hubble rate today. The degrees of freedom g ∗ and g ∗ ,s in Eq. (35) are evaluated at T = T ∗ , and g eq ∗ ,s is the number of degrees of freedom contributing to the radiation entropy at the time of matter-radiation equality. The production of GWs from sound waves stops after a period τ sw , when the plasma motion turns turbulent (Ellis et al. 2019a,b, 2020; Guo et al. 2021). In Eq. (34), this effect is taken into account by the suppression factor  \nΥ( τ sw ) = 1 -(1 + 2 τ sw H ∗ ) -1 / 2 , (36)  \nwhere the shock formation time scale, τ sw , can be written in terms of the phase transition parameters as τ sw ≈ R ∗ / ¯ U f , with ¯ U 2 f ≈ 3 κ s α ∗ / [4(1 + α ∗ )] (Weir 2018).  \nThe functions S b,s describe the spectral shape of the GWB and are expected to peak at the frequencies  \nf b,s ≃ 48 . 5 nHz g 1 / 2 ∗ ( g eq ∗ ,s g ∗ ,s ) 1 / 3 ( T ∗ 1 GeV ) f ∗ b,s R ∗ H ∗ R ∗ , (37)  \nwhere the values of the peak frequencies at the time of GW emission are given by f ∗ b = 0 . 58 /R ∗ (Jinno & Takimoto 2017) and f ∗ s = 1 . 58 /R ∗ (Hindmarsh et al. 2017). In passing, we mention that the numerical factors in these estimates may still change in the future, as our understanding of cosmological phase transitions  \nimproves. However, at the level of our Bayesian fit analysis, changes in these prefactors can be absorbed in the temperature scale T ∗ , which in its role as an independent fit parameter only controls the peak frequencies in Eq. (37). A similar argument applies to the numerical factors in Eqs. (33) and (34): changes in these prefactors can always be absorbed in a rescaled version of the fit parameter α ∗ , which only appears in the expressions for the peak amplitudes of the GWB signal.  \nThe shape of the spectral functions can be usually approximated with a broken power law of the form  \nS ( x ) = 1 N ( a + b ) c ( bx -a/c + ax b/c ) c . (38)  \nHere a and b are two real and positive numbers that give the slope of the spectrum at low and high frequencies, respectively; c parametrizes the width of the peak. The normalization constant, N , ensures that the logarithmic integral of S is normalized to unity and is given by  \nN = ( b a ) a/n ( nc b ) c Γ( a/n ) Γ ( b/n ) n Γ( c ) , (39)  \nwhere n = ( a + b ) /c and Γ denotes the gamma function. While the values of the coefficients a , b , and c can in principle be estimated from numerical simulations, we allow them to float within the prior ranges given in Table 3. These prior ranges were chosen to roughly reflect the typical uncertainties of numerical simulations and any possible model dependency of these coefficients (see, e.g., Hindmarsh et al. (2017, 2021); Lewicki & Vaskonen (2020, 2021); Cutting et al. (2018, 2021)). 4', '5.5. Domain walls MODEL DESCRIPTION': 'Domain walls are 2D topological defects that form when a cosmological phase transition results in the spontaneous breaking of a discrete symmetry, filling the universe with patches in different degenerate vacua (Kibble 1976). In the absence of significant interactions with relativistic particles, domain wall networks are expected to reach a scaling regime in which each Hubble volume, H -3 , contains A ∼ O (1) domain walls (see, e.g., Press et al. (1989); Hindmarsh (1996, 2003); Garagounis & Hindmarsh (2003)). In this regime, the energy density stored in domain walls is given by  \nρ DW ∼ A σH , (53)  \nwhere A is nearly constant and σ is the domain wall tension, which gives the domain walls their surface energy density. During the scaling regime, the energy density stored in domain walls dilutes more slowly than that of relativistic radiation and nonrelativistic matter. Indeed, the total energy density of the background always scales like ρ c ∝ H 2 , which means that Ω DW = ρ DW /ρ c grows like the Hubble radius when the domain wall network is in the scaling regime, Ω DW ∝ H -1 . Therefore, domain walls eventually overclose the universe and alter the cosmological evolution in a way that is incompatible with CMB observations (Zeldovich et al. 1974). A possible solution to this problem is to have domain walls decay at some temperature T ∗ by assuming, e.g., that the global symmetry responsible for domain wall formation is explicitly broken. Indeed, an explicit breaking of the symmetry introduces a bias among the possible low-energy vacua, thus lifting their degeneracy, which eventually leads to the collapse of the domain wall network (Vilenkin 1981b; Gelmini et al. 1989; Larsson et al. 1997).  \nDuring the scaling regime, domain walls continuously change their configuration and shrink in order to maintain the scaling relation in Eq. (53). While these pro-  \ncesses take place, a fraction of the domain wall energy is released in the form of GWs, which produce a GWB with a present-day relic abundance given by (Vilenkin 1981b; Preskill et al. 1991; Gleiser & Roberts 1998; Chang et al. 1999; Hiramatsu et al. 2010; Kawasaki & Saikawa 2011)  \nh 2 Ω GW ( f ) = 3 32 π D ˜ ϵ α 2 ∗ S ( f/f p ) , (54)  \nwhere D is the dilution factor defined in Eq. (35), ˜ ϵ = 0 . 7 is an efficiency coefficient derived from numerical simulations (Hiramatsu et al. 2014), and α ∗ is the fraction of the total energy density stored in domain walls at T ∗ ,  \nα ∗ ≡ Ω DW ( T ∗ ) = ρ DW 3 H 2 M 2 Pl ∣ ∣ ∣ ∣ T = T ∗ . (55)  \nThe GWB is dominated by the emission taking place just before the decay of the domain wall network. As the GWs are predominantly sourced by horizon-scale structures in the network, the typical frequency of the GWs emitted around this time is H ∗ , i.e., the Hubble scale at T ∗ . Therefore, after redshifting until today, the GWB is expected to exhibit a peak frequency  \nf p = 1 . 14 nHz ( 10 . 75 g ∗ ,s ) 1 / 3 ( g ∗ 10 . 75 ) 1 / 2 ( T ∗ 10 MeV ) , (56)  \nwith both g ∗ and g ∗ ,s evaluated at T ∗ . For the spectral shape, S ( x ), we use a parameterization similar to the one used for the phase transition spectrum in Section 5.3: 8  \nS ( x ) = ( a + b ) c ( bx -a/c + ax b/c ) c . (57)  \nCausality fixes a = 3, while numerical simulations (Hiramatsu et al. 2014) suggest b ≃ c ≃ 1. In our analysis, we therefore fix the low-frequency slope of the domain wall signal to the value predicted by causality, and, following Ferreira et al. (2022), we allow b and c to float within the uniform prior ranges b ∈ [0 . 5 , 1] and c ∈ [0 . 3 , 3].  \nIn our analysis we consider two possible decay channels for the domain wall network: dark radiation ( dwdr ) and SM particles ( dw-sm ). In the dw-dr model, instead of using α ∗ , we parameterize the strength of the domain wall signal in terms of the dark radiation produced in the decay, ρ DR . This quantity is usually expressed in terms of the effective number of extra neutrino species, ∆ N eff , which we defined in Section 5.1 and which is bounded from above by BBN and CMB observations, ∆ N max eff ∼ few × 0 . 1 (Pisanti et al. 2021; Yeh  \n∗  \n∗  \nFigure 12. Same as in Fig. 5, but for domain walls decaying to SM particles (left panel) and a dark sector (right panel). Fig. 31 in Appendix C.5 shows extended versions of the two plots that include the spectral shape parameters b , c and the SMBHB parameters A BHB and γ BHB .  \n<!-- image -->  \n<!-- image -->  \n∗  \net al. 2021). In our MCMC analysis of the dw-dr model, we follow Ferreira et al. (2022) and use ∆ N max eff = 0 . 39 as the upper bound on our prior for the parameter ∆ N eff . Assuming that all the domain wall energy is converted to dark radiation at T ∗ , ∆ N eff is related to α ∗ by  \n∆ N eff ≃ 0 . 62 ( 10 . 75 g ∗ ) 1 / 3 ( g ∗ ,s g ∗ )( α ∗ 0 . 1 ) , (58)  \nwhere, as before, both g ∗ and g ∗ ,s are evaluated at T ∗ . In the dw-sm model, BBN restricts the possible values of the decay temperature to T ∗ ≳ 2 . 7 MeV (Jedamzik 2006; Bai & Korwar 2022) for any detectable value of α ∗ . Following Ferreira et al. (2022), we also impose α ∗ < 0 . 3 to avoid any possible deviation from radiation domination and to evade bounds from ∆ N eff .', '6. DETERMINISTIC SIGNALS FROM NEW PHYSICS': 'In addition to the GWB signals discussed previously, there are several new-physics theories that can imprint a deterministic signal, described by a time series h , in pulsar timing data. In this section, we consider the deterministic signals induced by ultralight dark matter (ULDM) and DM substructures. After finding no statistically significant evidence for such signals in our data, we report upper limits on the allowed strength of these signals.', '6.1. Ultralight dark matter MODEL DESCRIPTION': "While DM constitutes roughly 27% (Zyla et al. 2020) of the energy density of the universe, very little is known about its fundamental properties. Consequently, a wide range of DM models remain consistent with cosmolog-  \nical and astrophysical observations. The lightest possible DM particles are classified as ultralight, or fuzzy, DM. These particles must be bosonic; otherwise, they could not be packed into galaxies owing to Fermi degeneracy pressure (Tremaine & Gunn 1979; Di Paolo et al. 2018; Savchenko & Rudakovskyi 2019; Alvey et al. 2020; Davoudiasl et al. 2021). These ULDM models also generically suppress structure on small scales, allowing them to be potential solutions to the small-scale structure problems of the standard ΛCDM paradigm (Bullock & Boylan-Kolchin 2017). However, too much suppression on large scales would be in conflict with CMB measurements (Hlozek et al. 2015; Hloˇzek et al. 2017, 2018), which sets a lower bound on the DM mass, 10 -24 eV < m ϕ . PTAs can probe these miniscule masses, since, as we describe below, the frequency of the ULDM signal, f , is generally proportional to the ULDM mass, 2 πf ∼ m ϕ . Therefore, the sensitivity window of NANOGrav is expected to be 10 -23 eV ≲ m ϕ ≲ 10 -20 eV.  \nOther astrophysical constraints, such as measurements of the Lymanα forest (Armengaud et al. 2017; Irˇsiˇc et al. 2017; Kobayashi et al. 2017; Rogers & Peiris 2021), galactic subhalo mass functions (Schutz 2020; Banik et al. 2021; Nadler et al. 2021), and stellar kinematics (Dalal & Kravtsov 2022), can push this bound further up, ranging from 10 -21 eV to 10 -19 eV, and the nonobservation of superradiance at SMBHs (Arvanitaki et al. 2015; Stott 2020; Unal et al. 2021) pushes the bound up to 10 -21 · · · 10 -17 eV. PTA searches for ULDM in the 10 -23 eV ≲ m ϕ ≲ 10 -20 eV window can then be viewed as complementary to these astrophysical searches. Signals in PTAs do not depend on the same astrophysical uncertainties, e.g., modeling the nonlinear small-scale matter power spectrum with analytic or numeric methods (Zhang et al. 2018, 2019) or the evolution of specific density profiles. Moreover, the power of these methods quickly degrades when considering subcomponents of DM, or populations of DM that do not compose the whole DM density. The searches discussed here have a weaker dependence on the subcomponent fraction and are therefore still useful in the hunt for DM.  \nWhile there is a large phenomenology of ULDM signals that can be produced in PTAs, the timing residuals for the I th pulsar, h I , are deterministic and can be written in the form  \nh I ( t ) = ∑ i A i E,I ( x E ) sin( ωt + γ i E ) + A i P,I ( x P,I ) sin( ωt + γ i P,I ) , (59)  \nTable 1. Summary of the different effects induced by ULDM that generate a deterministic signal of the form given in Eq. (59). 'E' and 'P' subscripts denote Earth and pulsar term contributions, respectively. A i is the signal amplitude for the i th DM polarization, ω is the signal frequency, and the 'Spin' column refers to whether the effect occurs for scalar (spin-0) or vector (spin1) ULDM candidates. ρ ϕ is the local DM density, taken to be 0 . 4 GeV / cm 3 , m ϕ is the ULDM mass, g are gauge couplings, Q is the charge under the corresponding gauge group, y parameterizes the sensitivity to different fundamental constant fluctuations, d are defined in Eq. (64), and Λ = M Pl / √ 4 π . ˆ ϕ ( x ) is a random variable representing the fluctuations in the ULDM field in different correlation patches; the correlated limit corresponds to ˆ ϕ ( x ) → ˆ ϕ . n I is a vector pointing from the Earth to the I th pulsar, and ϵ i are the ULDM polarization vectors.  \nwhere i indexes the DM field polarization, 9 ω is the signal frequency, and we have split the signal into two terms: an 'Earth term' with amplitude A i E,I ( x E ) and time-independent phase γ i E and a 'pulsar term' with amplitude A i P,I ( x P,I ) and time-independent phase γ i P,I . The phases can be written as  \nγ i E = α i ( x E ) , (60)  \nγ i P,I = ω | x E -x P,I | + α i ( x P,I ) , (61)  \nwhere α i ( x ) is dependent on the underlying ULDM field phase at point x , and the additional term in γ i P,I is due to the light-travel time between the Earth and the pulsar.  \nAn important scale in understanding these signals is given by the correlation length of the ULDM field, ℓ c ,  \nℓ c ≃ 2 π m ϕ v ϕ ∼ 0 . 4 kpc ( 10 -22 eV m ϕ ) , (62)  \nwhere v ϕ ∼ 10 -3 is the ULDM velocity. If the correlation length is much larger than the distance between Earth and the pulsar, ℓ c ≫ | x E -x P,I | , both experience the same DM field. In this 'correlated' limit of the signals, the amplitudes and phases can be taken to  \nbe position independent:  \nA i E,I ( x E ) → A i E , A i P ( x P,I ) → A i P,I , α i ( x ) → α i . (63)  \nGenerally, a correlated analysis can drastically reduce the number of free amplitude parameters in the MCMC. For example, if ULDM couples to all pulsars identically, then there is only one amplitude, A P , instead of one for each pulsar, to account for fluctuations in the ULDM field. However, since m ϕ | x E -x P,I | ≫ 1, γ P,I and γ E must still be taken to be independent. In the uncorrelated limit, ℓ c ≪| x E -x P,I | , the DM field is no longer correlated, and its amplitude fluctuations in different patches need to be accounted for.  \nIn Table 1, we summarize the ULDM signals searched for in our analysis by providing the corresponding expression for the amplitudes and signal frequency appearing in Eq. (59). In the following, we give a brief review of each of these signals and refer the reader to the original references for more details:  \n- · Metric fluctuations (Khmelnitsky & Rubakov 2014; Porayko & Postnov 2014; Porayko et al. 2018; Kato & Soda 2020; Nomura et al. 2020; Unal et al. 2022; Wu et al. 2022): Oscillations in the ULDM field generate fluctuations in the local stress-energy tensor, T µν , which are independent of any direct couplings the ULDM may have with SM fields. These fluctuations in T µν generate fluctuations in metric perturbations by Einstein's equations. These metric perturbations then affect the photon geodesic on its path from the pulsar and generate a timing residual with pulsar and Earth amplitudes given in Table 1. While the amplitudes shown in Table 1 are  \nspecific to scalar (spin-0) ULDM, vector (spin-1) ULDM also generates this purely gravitational signal. The vector DM scenario is more complicated (Nomura et al. 2020), as off-diagonal components of the metric fluctuations can generate nontrivial angular correlations between signals seen in different pulsars. In our analysis, we ignore these spatial correlations and set approximate bounds by treating the vector as three scalar components, such that the vector signal is three times larger than the scalar one (Nomura et al. 2020).  \n- · DopplerU (1) forces (Graham et al. 2016; Xue et al. 2022): A background of ultralight vector DM can create dark 'electric' fields, which can accelerate pulsars, or the Earth, and Doppler-shift light arrival times. For example, if the DM is the gauge field of a local baryon symmetry, U (1) B , then it couples to the Earth and pulsar baryon current density. This coupling generates a force on both Earth and pulsars proportional to their baryon number, Q B = N B × q B , where N B is the number of baryons and q B is the gauge charge. This scenario is straightforwardly generalized to other U (1) models, e.g., models where the difference between baryon and lepton numbers, B -L , is a local symmetry. These forces cause periodic displacements between the Earth and pulsars and generate timing residuals with amplitudes given in Table 1. These forces also exist for scalar DM coupling to SM operators (e.g., ϕ ¯ nn , where ϕ is the DM and n is the neutron field). However, in the scalar case, these forces originate from the field gradient ∇ ϕ (Graham et al. 2016), which causes them to be velocity suppressed compared to the vector DM scenario.  \nThe next two signals we discuss are specific to scalar ULDM. All the linear couplings of a scalar ULDM field to the SM can be summarized in a single Lagrangian,  \nL ⊃ ϕ Λ [ d γ 4 e 2 F µν F µν + d g β 3 2 g 3 G A µν G µν A -∑ f = e,µ d f m f ¯ ff -∑ q = u,d ( d q + γ q d g ) m q ¯ qq ] , (64)  \nwhere F µν and G µν are the photon and gluon field strengths, respectively, Λ = M Pl / √ 4 π , β 3 is the QCD beta function, γ q are the light quark anomalous dimensions, m f are the fermion masses, and the d values are dimensionless couplings. These couplings induce periodic oscillations in the values of fundamental constants, i.e., particle masses and couplings, which can affect timing residuals in two ways:  \n- · Pulsar spin fluctuations (Kaplan et al. 2022): Particle mass fluctuations change the moment of inertia of the pulsar, which, by conservation of angular momentum,  \nleads to pulsar spin fluctuations. The amplitude of the signal is given in Table 1. The y parameters are pulsar specific and denote the pulsar sensitivity to a specific coupling, e.g., y e parameterizes the pulsar sensitivity to changes in the electron mass. y ˆ m parameterizes the pulsar sensitivity to the mass-weighted combination of quark mass couplings, d ˆ m = ( m u d u + m d d d ) / ( m u + m d ). To be explicit, following Kaplan et al. (2022), we employ in our analysis  \ny g = -5 , y ˆ m = -2 . 4 × 10 -1 , y µ = 2 × 10 -3 , y e = 1 . 7 × 10 -5 , (65)  \nfor each pulsar, which assumes the simplest possible model for the pulsars' moment of inertia.  \n- · Reference clock shifts (Graham et al. 2016; Kaplan et al. 2022): PTA timing residuals are measured with respect to a collection of mostly cesium (McCarthy 2009) atomic clocks. Therefore, changes to these atomic clock frequencies can appear in PTA data as apparent shifts in timing residuals. These shifts are perfectly correlated among pulsars and have the amplitude reported in Table 1. Similar to the pulsar spin fluctuation signal, the y parameters denote the atomic clock sensitivity to specific couplings. To be explicit, following Kaplan et al. (2022), we take these parameters to be  \ny g = 1 y ˆ m = 0 . 158 y e = 2 y γ = 4 . 83 , (66)  \nassuming cesium clocks to be used as a reference.", '7. DISCUSSION': 'The analysis of the NANOGrav 15-year data set has produced the first convincing evidence for a stochastic background of GWs in the nanohertz frequency range. The origin of this background is still unknown. In this work, we considered various cosmological sources and compared them to the commonly studied astrophysical signal produced by a population of inspiraling supermassive black-hole binaries. Specifically, we considered  \nthe signals produced by nonminimal inflation scenarios, scalar-induced GWs, cosmological phase transitions, several cosmic-string models, and domain walls.  \nFor each of these models, we identified regions of the parameter space that are compatible with the observed signal. We find that, with the exception of stable cosmic strings of field theory origin, all new-physics models considered in this analysis are capable of reproducing the GWB signal. Many models allow us in fact to fit the signal better than the SMBHB reference model, which is reflected in Bayes factors ranging from 10 to 100 (see Fig. 2). When the new-physics signals are combined with the astrophysical one, we obtain comparable results. More precisely, in several models, the addition of the SMBHB signal leads to a slight decrease of the Bayes factor, which indicates that the SMBHB contribution does not help to improve the quality of the fit but merely increases the prior volume. In other models, on the other hand, such as pt-sound and dw-dr , adding the SMBHB signal on top of the new-physics signal can lead to a slight increase of the Bayes factor, indicating that the SMBHB signal can in fact play a dominant role in the total GW spectrum. For all four stable-string models, we find Bayes factors between 0 . 1 and 1. Cosmic superstrings, on the other hand, which are also stable but not of field theory origin, can explain the data at a level comparable to other new-physics sources.  \nDespite the fact that some of the Bayes factors derived in this paper might suggest that a purely astrophysical interpretation of the signal is disfavored by the data, we caution against this interpretation. The Bayes factors do not account for the full range of uncertainties in both the cosmological and astrophysical signals and are prior dependent. It is conceivable that, as our understanding of SMBHB signals and our noise models improve, the tension between observations and astrophysical predictions will decrease, potentially weakening the evidence in favor of cosmological signals.  \nFuture data sets will improve the spectral characterization of the signal and improve our ability to discriminate cosmological sources from the SMBHB signal. Unfortunately, similarities in the spectral shape and theoretical uncertainties will make it challenging to definitively determine the origin of the background using power spectrum characterization alone. However, the observation of anisotropies could eventually resolve this debate, as the expected anisotropies generated by black hole binaries are significantly larger than those produced by most cosmological sources. Similarly, the detection of a continuous wave from a single binary would provide convincing evidence in favor of an astrophysical origin of the signal. Ultimately, measurements of the GWB spec-  \ntral shape and angular power spectrum may be complemented by observations of its polarization content and possible deviations from Gaussian statistics, which can again help to discriminate between an astrophysical and a cosmological origin of the signal.  \nIt is worth emphasizing that in all parts of our analysis we described cosmological signals using effective parameters, e.g., T ∗ , α ∗ , and H ∗ R ∗ , for the phase transition models. Moving forward, it will be crucial to identify microscopic models that can reproduce the values of these parameters that we found to best fit the GWB signal. That is, in order to shed more light on the various cosmological interpretations of the signal, we require a better understanding of how the NANOGrav signal could possibly result from the fundamental parameters of a particle physics Lagrangian that describes the generation of GWs in the early universe.  \nAlong with searching for a cosmological GWB, we also analyze our data to search for deterministic signals generated by models of ULDM and DM substructures. We do not find significant evidence for either of these signals. Nonetheless, we are able to place constraints on the parameters space of these models. For a wide range of ULDM models, our constraints compete or outperform laboratory constraints in the 10 -23 eV ≲ m ϕ ≲ 10 -20 eV mass window. The signal from DM substructures is harder to detect; as a consequence, we are able to set very weak constraints on the local abundance of these objects. Future data sets will improve our reach, but a better characterization of the GWB will be needed to probe realistic models of DM substructures.  \nThe discovery of a GWB will lead to significant breakthroughs in our understanding of cosmology and particle physics. As future PTA data sets become available, we will establish the origin of the GWB. Regardless of whether the signal is of cosmological origin, we have shown how PTAs will undoubtedly contribute to exploring new physics, either as a discovery tool or as a new way to constrain the parameter space of BSM models.', 'ACKNOWLEDGMENTS': "Author contributions. An alphabetical-order author list was used for this paper in recognition of the fact that a large, decade timescale project such as NANOGrav is necessarily the result of the work of many people. All authors contributed to the activities of the NANOGrav Collaboration leading to the work presented here, and reviewed the manuscript, text, and figures prior to the paper's submission. Additional specific contributions to this paper are as follows.  \nA.Mi. proposed and initiated the project. A.Mi. and K.Sc. coordinated and led the search. The enterprise  \nwrapper, PTArcade (Mitridate et al. 2023b), used in this analysis was mainly developed by A.Mi., with help from D.W., K.D.O., J.N., R.v.E., T.S., and T.T. R.v.E. and T.S. developed the statistical tools needed to derive the K -ratio bounds. L.Z.K. derived the distributions of power-law fit parameters for characteristic strain spectra from SMBHB simulations. A.Mi. prepared all figures except those in Section 6.1, which were prepared by T.T.; Fig. 18, which was prepared by V.L.; and Fig. 22, which was prepared by R.R.L.d.S. A.Mi. and K.Sc. wrote the paper with help from K.B., K.D.O., S.Ve., and T.T.  \nContributions to specific analyses are as follows. K.Sc. led and discussed the results of the SIGW and IGW analyses, which were performed by D.W. and R.R.L.d.S. A.Af. and R.R.L.d.S. derived the constraints on the parameter space of the SIGW model. D.W. and R.R.L.d.S. derived the LVK and N eff bounds on the IGW model. A.Mi. led the PT analysis, which was performed by both A.Mi. and R.v.E.; A.Mi. and K.Sc. interpreted and discussed the results. The string analyses were led by K.D.O. and K.Sc., with K.D.O., K.Sc., and T.S. conducting the analyses; K.Sc. interpreted and discussed the results with help from K.D.O. and J.J.B.P. A.Mi. led the domain wall analysis, which was performed by both A.Mi. and D.W.; A.Mi. and K.Sc. interpreted and discussed the results. The ULDM analysis was coordinated by A.Mi., K.B., and T.T.; J.N. performed the analysis; and A.Mi., C.U., K.B., J.N., and T.T. interpreted and discussed the results. Finally, the substructure search was led by A.Mi. and T.T.; A.Mi., S.Ve., and V.L. performed the analysis; and A.Mi., S.Ve., V.L., and T.T. interpreted and discussed the results.  \nAcknowledgments. The NANOGrav Collaboration receives support from National Science Foundation (NSF) Physics Frontiers Center award Nos. 1430284 and 2020265, the Gordon and Betty Moore Foundation, NSF AccelNet award No. 2114721, an NSERC Discovery Grant, and CIFAR. The Arecibo Observatory is a facility of the NSF operated under cooperative agreement (AST-1744119) by the University of Central Florida (UCF) in alliance with Universidad Ana G. M'endez (UAGM) and Yang Enterprises (YEI), Inc. The Green Bank Observatory is a facility of the NSF operated under cooperative agreement by Associated Universities, Inc. The National Radio Astronomy Observatory is a facility of the NSF operated under cooperative agreement by Associated Universities, Inc. NANOGrav is part of the International Pulsar Timing Array (IPTA); we would like to thank our IPTA colleagues for their help with this paper.  \nPart of this work was conducted using the High Performance Computing Cluster PALMA II at the Uni-  \nversity of Munster (https://www.uni-muenster.de/IT/ HPC). This work used the Maxwell computational resources operated at Deutsches Elektronen-Synchrotron DESY, Hamburg (Germany). This work was conducted in part using the HPC resources of the Texas Advanced Computing Center (TACC) at the University of Texas at Austin. The Tufts University High Performance Computing Cluster (https://it.tufts.edu/ high-performance-computing) was utilized for some of the research reported in this paper. This research used the computational resources provided by the University of Central Florida's Advanced Research Computing Center.  \nJ.J.B.P. acknowledges the support by the PID2021-123703NB-C21 grant funded by MCIN/AEI /10.13039/501100011033/ and by ERDF; 'A way of making Europe,' the Basque Government grant (IT1628-22); and the Basque Foundation for Science (IKERBASQUE). L.B. acknowledges support from the National Science Foundation under award AST-1909933 and from the Research Corporation for Science Advancement under Cottrell Scholar Award No. 27553. P.R.B. is supported by the Science and Technology Facilities Council, grant No. ST/W000946/1. S.B. gratefully acknowledges the support of a Sloan Fellowship and the support of NSF under award No. #1815664. M.C. and S.R.T. acknowledge support from NSF AST-2007993. M.C. and N.S.P. were supported by the Vanderbilt Initiative in Data Intensive Astrophysics (VIDA) Fellowship. Support for this work was provided by the NSF through the Grote Reber Fellowship Program administered by Associated Universities, Inc./National Radio Astronomy Observatory. Support for H.T.C. is provided by NASA through the NASA Hubble Fellowship Program grant No. #HST-HF2-51453.001 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS526555. K.C. is supported by a UBC Four Year Fellowship (6456). M.E.D. acknowledges support from the Naval Research Laboratory by NASA under contract S-15633Y. T.D. and M.T.L. are supported by an NSF Astronomy and Astrophysics Grant (AAG), award No. 2009468. The work of R.v.E. K.Sc., and T.S. is supported by the Deutsche Forschungsgemeinschaft (DFG) through the Research Training Group, GRK 2149: Strong and Weak Interactions - from Hadrons to Dark Matter. E.C.F. is supported by NASA under award No. 80GSFC21M0002. G.E.F., S.C.S., and S.J.V. are supported by NSF award PHY-2011772. The Flatiron Institute is supported by the Simons Foundation. A.D.J. and M.V. acknowledge support from the  \nCaltech and Jet Propulsion Laboratory President's and Director's Research and Development Fund. A.D.J. acknowledges support from the Sloan Foundation. The work of N.La. and X.S. is partly supported by the George and Hannah Bolinger Memorial Fund in the College of Science at Oregon State University. N.La. acknowledges the support from Larry W. Martin and Joyce B. O'Neill Endowed Fellowship in the College of Science at Oregon State University. Part of this research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (80NM0018D0004). V.S.H.L. is supported by the DoE under contract DE-SC0011632. R.R.L.d.S. is supported by a research grant (29405) from VILLUM FONDEN. D.R.L. and M.A.M. are supported by NSF No. #1458952. M.A.M. is supported by NSF #2009425. C.M.F.M. was supported in part by the National Science Foundation under Grants Nos. NSF PHY-1748958 and AST-2106552. A.Mi. is supported by the Deutsche Forschungsgemeinschaft under Germany's Excellence Strategy - EXC 2121 Quantum Universe - 390833306. The Dunlap Institute is funded by an endowment established by the David Dunlap family and the University of Toronto. K.D.O. was supported in part by NSF grants Nos. 2111738 and 2207267. T.T.P. acknowledges support from the Extragalactic Astrophysics Research Group at Eotvos Lor'and University, funded by the Eotvos Lor'and Research Network (ELKH), which was used during the development of this research. S.M.R. and I.H.S. are CIFAR Fellows. Portions of this work performed at NRL were supported by ONR 6.1 basic research funding. J.D.R. also acknowledges support from start-up funds from Texas Tech University. J.S. is supported by an NSF Astronomy and Astrophysics Postdoctoral Fellowship under award AST-2202388 and acknowledges previous support by the NSF under award 1847938. S.R.T. acknowledges support from an NSF CAREER award No. #2146016. T.T. contribution to this work is supported by the Fermi Research Alliance, LLC, under contract No. DE-AC02-07CH11359 with the U.S. Department of Energy, Office of Science, Office of High Energy Physics. C.U. acknowledges support from BGU (Kreitman fellowship) and the Council for Higher Education and Israel Academy of Sciences and Humanities (Excellence fellowship). C.A.W. acknowledges support from CIERA, the Adler Planetarium, and the Brinson Foundation through a CIERA-Adler postdoctoral fellowship. O.Y. is supported by the National Science Foundation Graduate Research Fellowship under grant No. DGE-2139292. K.Z. is supported by the  \nDoE under contract DE-SC0011632, and by a Simons Investigator award.", 'A. PARAMETER RANGES AND LIMITS': 'In this appendix, we specify the prior assumptions for all model parameters used in our analyses and report characteristic values for these parameters that we extract from the corresponding reconstructed 1D marginalized posterior distributions. In Table 2, we list our prior assumptions for the pulsar-intrinsic red-noise parameters A red and γ red ( A a and γ a in Eq. (4)) as well as for the SMBHB parameters A BHB and γ BHB . In the latter case, we work with a bivariate normal distribution for (log 10 A BHB , γ BHB ) whose mean and covariance matrix are given by  \nµ BHB = ( -15 . 6 4 . 7 ) , σ BHB = 10 -1 × ( 2 . 8 -0 . 026 -0 . 026 1 . 2 ) , (A1)  \nwhich we obtain by fitting the log 10 A BHB and γ BHB distributions extracted from the SMBHB simulations in the GWOnly-Ext library in NG15smbh (see Section 4). In Table 3, we list our prior assumptions for the model parameters of all new-physics models considered in this work, and in Table 4, we summarize various key values extracted from the corresponding 1D marginalized posterior distributions. Specifically, we report the Bayes estimator, including its standard deviation; the maximum posterior estimator, including the 68% credible interval; and (when applicable) the upper bound based on the requirement that the K ratio in Eq. (10) should not be smaller than 1 / 10 .  \nThe Bayes estimator ⟨ θ ⟩ of a parameter θ with marginalized 1D posterior probability distribution P ( θ |D , H ) corresponds to the expectation value with respect to the distribution P ( θ |D , H ), while the standard deviation σ θ of the Bayes estimator corresponds to the positive square root of the associated variance σ 2 θ ,  \n⟨ θ ⟩ = ∫ dθ θ P ( θ |D , H ) , σ 2 θ = ⟨ θ 2 ⟩ - ⟨ θ ⟩ 2 = [∫ dθ θ 2 P ( θ |D , H ) ] -[∫ dθ θ P ( θ |D , H ) ] 2 . (A2)  \nIn practice, in a given analysis and for a given chain of MCMC samples, we compute the Bayes estimator and its standard deviation in terms of the corresponding sample mean and sample variance. The maximum posterior estimator θ max of a parameter θ with marginalized 1D posterior probability distribution P ( θ |D , H ) corresponds to the θ value where P ( θ |D , H ) reaches its global maximum across the predefined prior range,  \nP ( θ max |D , H ) = max θ P ( θ |D , H ) , (A3)  \nand the 68% Bayesian credible interval [ θ min 68 , θ max 68 ] follows from integrating the posterior distribution P ( θ |D , H ) over the regions of highest posterior density such that the integral returns an integrated probability of 68%,  \n∫ θ max 68 θ min 68 dθ P ( θ |D , H ) = 0 . 68 , (A4)  \nwhere P ( θ |D , H ) > P 68 for all θ ∈ [ θ min 68 , θ max 68 ] and some appropriate threshold P 68 . Similarly, we can also construct 95% Bayesian credible intervals. Finally, we mention again that the K -ratio bound in the last column of Table 4 indicates the value θ K of the parameter θ that returns K = 1 / 10 (see Section 3.2),  \nK ( θ K ) = P ( D| θ K , H ) P ( D| θ 0 , H ) = 1 10 . (A5)  \nNote that, unlike all other quantities discussed in this section, the K ratio is not defined in terms of a posterior probability distribution, but rather in terms of a likelihood ratio, which makes it robust against our prior assumptions.  \nTable 2 . Prior distributions for the pulsar-intrinsic red-noise parameters and the parameters of the astrophysical SMBHB signal. The mean and covariance matrix of the Gaussian prior distribution for (log 10 A BHB , γ BHB ) are given in Eq. (A1).  \nTable 3 . Priors distributions for the parameters of the new-physics models considered in this work. The ∗ indicates parameters that are present only in the uncorrelated ULDM analyses.  \n∗', 'Cosmic Superstrings (super)': '-  \n11  \n.  \n67  \n±  \n0  \n.  \n32  \n-  \n11  \n.  \n68  \n±  \n0  \n.  \n35  \n-  \n11  \n.  \n94  \n-  \n11  \n.  \n96  \n[  \n-  \n12  \n.  \n08  \n,  \n-  \n11  \n.  \n50]  \n[  \n-  \n12  \n.  \n09  \n,  \n-  \n11  \n.  \n49]  \n-  \n9  \n.  \n97', 'B. MEDIAN GW POWER SPECTRA': 'In Figs. 3 and 4 in the main text, as well as in Figs. 19 and 20 in this appendix, we show median GW power spectra for all of the new-physics models considered in this work. The purpose of this appendix is to explain how these spectra are constructed. The main idea is to take the output of our Bayesian analysis, i.e., the reconstructed posterior distributions in the parameter space of the respective models, and map these distributions onto distributions of h 2 Ω GW values at each frequency in the NANOGrav frequency band. In practice, this means that we sample model parameter values from our posterior distributions and use these parameter values to evaluate the GW power spectrum at one frequency after another in small steps in the range from f = 1 /T obs to f = 30 /T obs . At each fixed frequency, we thus obtain an induced posterior distribution of h 2 Ω GW values, from which we can derive point values and uncertainty estimates. The median h 2 Ω GW value of a distribution at fixed f defines the value of the median GW spectrum at this frequency. Similarly, the equal-tailed 68% and 95% intervals around the median values provide an uncertainty estimate for the median GW spectrum; these intervals are shown in terms of the blue and red bands in Figs. 19 and 20.  \nWe stress that the median GW spectra in Figs. 3, 4, 19, and 20 represent effective power spectra that encapsulate global properties of the underlying distributions of actual GW spectra. In general, no individual GW spectrum at a given point in parameter space will exactly coincide with the median GW spectrum. For most models, the difference between the individual GW spectra and the median GW spectrum is rather mild. However, for some models, there can be noticeable differences, such as in the case of the SIGW models. Note, e.g., that the median GW spectrum of the sigw-delta model in Fig. 4 does not have the characteristic double-peak structure that each individual GW spectrum in this model has. The reason for this is clear: the median GW spectrum follows from a distribution of many individual GW spectra whose peaks are located at different frequencies. The double-peak structure of the GW signal in the sigw-delta model is therefore washed out owing to the averaging over many individual GW spectra.  \nAnother caveat is that median GW spectra violating an experimental bound (e.g., the LVK bound) do not automatically indicate that the corresponding model is ruled out. Again, as the median GW spectrum is constructed from a distribution of individual GW spectra, the violation of an experimental bound merely signals that some (maybe most)  \nFigure 19. Median GWB spectra (solid lines) for all new-physics models considered in this work except for the cosmic-string models, together with their 68% and 95% posterior envelopes. Median GWB spectra for the cosmic-string models are shown in Fig. 20. In the left column (blue shading), we show the median GWB spectra for the new-physics models alone; in the right column (red shading), we combine the new-physics signals with the signal from SMBHBs. The gray violins are symmetrical representations of the 1D marginalized posterior probability density distributions of the GW energy density at each sampling frequency of the data. The dashed black lines show the GWB spectrum produced by inspiraling SMBHBs (see caption of Fig. 3).  \n<!-- image -->  \n7  \n4  \n4  \nFigure 20. Same as Fig. 19 but for the cosmic-string models considered in this work. See Appendix B for details.  \n<!-- image -->  \nindividual GW spectra are in conflict with the experimental data. It is, however, well possible that some fraction of the underlying distribution of individual GW spectra is not ruled out and is perfectly consistent with all experimental constraints. This is, e.g., true for the GW spectra from cosmic superstrings. The median GW spectrum of the super model violates the LVK bound (see Fig. 4). However, at the level of the model parameter space, this merely means that some parameter regions are experimentally ruled out, while other regions remain viable (see Fig. 11). Another example is the median GWB spectrum of the igw model, which appears to violate the N eff bound if it is extended to high frequencies beyond the NANOGrav band (see Fig. 4). However, the igw model as a whole is not ruled out, as  \nFigure 21. Same as Fig. 5 but including the SMBHB parameters A BHB and γ BHB . The black dashed lines in the three lowerright panels show the prior distributions for A BHB and γ BHB , i.e., one 2D Gaussian and two 1D Gaussian distributions.  \n<!-- image -->  \n∗  \none can still choose the number of e -folds during reheating, N rh , for each individual GW spectrum. In some regions of parameter space, the N eff bound then persists to pose a problem, despite this additional parametric freedom and even in the limit N rh → 0; other regions, however, remain phenomenologically viable (see Appendix C.1 and Fig. 22).', 'C.1. Cosmic inflation': 'In Fig. 22, we present constraints on the parameter space of inflationary GWs, i.e., the igw model discussed in Section 5.1. This parameter space is spanned by four quantities: the tensor-to-scalar ratio r , the spectral index of the primordial tensor power spectrum n t , the reheating temperature T rh , and the number of e -folds during reheating N rh . However, because of the strong covariance between n t and r in Fig. 5, we are able to eliminate r and work on the 3D hypersurface in parameter space where r is given by the linear fit in Eq. (22). The effective parameter space spanned by T rh , n t , and N rh is then subject to two constraints: (i) the upper limit on the allowed amount of extra relativistic radiation, ∆ N eff , at the time of BBN and CMB decoupling (see Eq. (23)) and (ii) the upper limit on the amplitude of the SGWB set by the LVK Collaboration based on their first three observing runs (see Eq. (24)). The N eff bound can be saturated for a critical value of the cutoff frequency f end in the GW spectrum, which, for given T rh , is solely controlled by the Hubble rate at the end of inflation, H end . This means that the N eff bound can be translated to a maximally allowed cutoff frequency f max end , which can then be converted to a maximally allowed Hubble rate H max end and ultimately an upper limit N max rh on the allowed number of e -folds during reheating. Similarly, if the predicted GW signal at f lvk = 25Hz should exceed the LVK bound, we can derive upper limits H max end and N end rh by requiring f end not to exceed the lower end of the LVK band; for definiteness, we use f max end = 20Hz.  \nOur results for N max rh obtained in this analysis are presented in Fig. 22 in terms of black and maroon contour lines, where the black solid lines refer to the N eff bound, while the maroon dashed lines refer to the LVK bound. In the left panel we show our results for N max rh itself, while in the right panel we compare N max rh to the duration of reheating that  \n<!-- image -->  \nFigure 22. Constraints on the parameters of inflationary GWs, i.e., the igw model discussed in Section 5.1. The black solid lines refer to the N eff bound, while the maroon dashed lines refer to the LVK bound. The left panel shows the results for N max rh , while the right panel compares N max rh to the duration of reheating that one would naively expect in models of standard single-field slow-roll inflation. The blue points correspond to the sampled points used to obtain the 2D marginalized posterior probability density for the model parameters T rh and n t in Section 5.1. See Appendix. C.1 for details.  \n<!-- image -->  \none would naively expect in models of standard single-field slow-roll inflation (see Eqs. (15) and (16)),  \nN max rh = 2 3 ln ( H max end H min end ) , N naive rh = 2 3 ln ( H naive end H min end ) , H naive end = ( π 2 2 rA s ) 1 / 2 M Pl , (C6)  \nwhere H min end denotes the minimal Hubble rate that is necessary to realize the desired value of T rh after inflation,  \nH min end = ( π 2 g rh ∗ 90 ) 1 / 2 T 2 rh M Pl . (C7)  \nIn other words, for given T rh , H min end represents the Hubble rate in the limit of instantaneous reheating, N rh = 0. In terms of the maximally allowed number of e -folds N max rh and the naive expectation N naive rh , we can then construct  \n∆ N rh = N max rh -N naive rh = 2 3 ln ( H max end H naive end ) , (C8)  \nwhich is the quantity shown in the right panel of Fig. 22. In addition to the contour lines for N max rh and ∆ N rh , we also display the samples in the T rh -n t plane that we obtain from our MCMC chains. The density of these blue points approximates the 2D marginalized posterior probability density for the model parameters T rh and n t .  \nIn view of the results in Fig. 22, we can draw several conclusions: (i) Both the N eff bound and the LVK bound are only relevant at large values of the spectral index. At T rh ∼ 10 -3 GeV, the number of e -folds during reheating is only constrained as long as n t ≳ 2 . 5, while at T rh ∼ 10 3 GeV, the bounds on N rh only appear at n t ≳ 1. At f ≫ f rh , these n t values translate to a slope of the GW spectrum α in the range from α ∼ -1 to α ∼ 0 . 5, which is in the regime where we expect the approximate expression in Eq. (25) to be valid (see also the discussion in Kuroyanagi et al. (2015)). (ii) The LVK bound only depends on T rh . This follows from the fact that it is derived from the requirement f max end = 20Hz, which is independent of r and n t . Typically, the LVK bound on N rh is much weaker than the N eff bound  \nFigure 23. Same as Fig. 6 but including the SMBHB parameters A BHB and γ BHB .  \n<!-- image -->  \non N rh ; only in the parameter region where the N eff bound begins to disappear, N max rh → ∞ , does the LVK bound becomes competitive. (iii) The N eff bound is particularly strong at large reheating temperatures and large values of the spectral index, where N max rh can become even negative. This region, which we indicate by a dark-gray shading, is phenomenologically not viable and hence ruled out (see also the regions labeled N eff in Fig. 5). In the excluded region, we also find ∆ N rh ≪ 0, which indicates that reheating lasting for the naive number of e -folds N naive rh will definitely result in a violation of the N eff bound. (iv) Farther away from the excluded region, the upper limits on N rh become weak very fast, N max rh ≫ 1. In realistic models of inflation and reheating, where we typically expect N rh ∼ O (1 · · · 10), these bounds should be easy to satisfy. We therefore, conclude that most of the region shaded in light gray as well as the entire white region in Fig. 22 can host viable realizations of the igw model. This includes, in particular, scenarios with a low reheating temperature, T rh ∼ 10 -(3 ··· 0) GeV, and a large spectral index of the primordial tensor power spectrum, n t ∼ 3 · · · 4.', 'C.2. Scalar-induced gravitational waves': "The generation of SIGWs relies on enhanced scalar perturbations on small scales, which can simultaneously lead to the production of PBHs. The parameter space of the SIGW models considered in this work is therefore subject to the constraint that the PBH mass density must not exceed the relic density of DM, f PBH ≤ 1 (see Eq. (32)). This appendix explains how we evaluate this bound on the SIGW model parameters. We use the Press-Schechter formalism for spherical collapse (Press & Schechter 1974; Carr 1975) and mainly follow Inomata et al. (2017b) and Ando et al. (2018). The main quantity of interest is f PBH , which denotes the total fraction of PBH DM,  \nf PBH = ∫ dM M f ( M ) . (C9)  \nwith PBH mass function f given as (Ando et al. 2018)  \nf ( M ) ≃ γ 3 / 2 ( β ( M ) 1 . 6 × 10 -9 )( 10 . 75 g ∗ ( T ) ) 1 / 4 ( 0 . 12 Ω DM h 2 )( M ⊙ M ) 1 / 2 . (C10)  \nFigure 24. Same as Fig. 7 but including the SMBHB parameters A BHB and γ BHB .  \n<!-- image -->  \n∗  \nHere M is the mass of the PBH, T is the temperature at the time of PBH formation, γ is the ratio between the PBH mass and the horizon mass, and β ( M ) is the PBH production rate. In this work, we only focus on PBHs produced during the radiation-dominated era. We choose γ = 0 . 356, following Choptuik (1993),Koike et al. (1995) and Niemeyer & Jedamzik (1998) and assume that the mass of the PBH is proportional to the horizon mass when k = aH ,  \nM = γ M H ≃ γ M ⊙ ( g ∗ ( T ) 10 . 75 ) -1 / 6 ( 4 . 2 × 10 6 Mpc -1 k ) 2 . (C11)  \nWe can then express the PBH mass as a function of the temperature using the standard temperature-wavenumber relation in ΛCDM. Assuming that the density perturbation follows Gaussian statistics, β ( M ) is given by  \nβ ( M ) = ∫ δ c dδ √ 2 π σ ( M ) exp ( -δ 2 2 σ 2 ( M ) ) ≃ σ ( M ) √ 2 π δ c exp ( -δ 2 c 2 σ 2 ( M ) ) , (C12)  \nwhere δ c is the critical density perturbation for the PBH formation, whose exact value depends on the shape of the curvature power spectrum and satisfies 0 . 4 ≲ δ c ≲ 0 . 6 (Musco et al. 2005, 2021), and σ 2 ( M ) is the variance of the coarse-grained density perturbation, which in the radiation domination era is given by  \nσ 2 ( k ) = 16 81 ∫ dq q ( q k ) 4 W 2 ( q k ) T 2 ( q, k -1 ) P R ( q ) . (C13)  \nHere W ( x ) is some window function smearing over k -1 and  \nT ( q, k -1 ) = 3 [ sin ( q √ 3 k ) -( q √ 3 k ) cos ( q √ 3 k )] / ( q √ 3 k ) 3 . (C14)  \nis the transfer function for the radiation-dominated era. In our work, we choose δ c = 0 . 45 as a fiducial value (Ando et al. 2018; Wang et al. 2019) and a Gaussian window function W ( x ) = e -x 2 / 2 .  \nOur analysis is subject to a few caveats. The correct value of the critical density δ c depends on the choice of the power spectrum and the window function (Ando et al. 2018; Musco 2019; Young 2019; Escriv'a 2020; Escriv'a et al. 2020; Gow et al. 2021; Musco et al. 2021; Dandoy et al. 2023), and the nonlinear relation between density perturbations and density contrast (De Luca et al. 2019; Young et al. 2019; Escriv'a 2020). In this work, we use a single value of δ c for three different power spectra for computational reasons. We also disregard any corrections from the QCD equation of state, which are expected to be important at the frequency range probed by PTA observations (Abe et al. 2021; Escriv'a et al. 2022a; Juan et al. 2022; Dandoy et al. 2023; Musco et al. 2023).  \nThe statistics of the primordial scalar perturbations also strongly affect the abundance of PBHs since, in general, In this work, we focus on Gaussian primordial scalar perturbations. With non-Gaussianities, the same PBH abundance can be produced with scalar per-  \nenhanced scalar perturbations come with different levels of non-Gaussianity. turbations whose amplitudes are orders of magnitude smaller than those produced with Gaussian scalar perturbations. Therefore, primordial non-Gaussianity can dramatically affect the PBH bounds (Young & Byrnes 2013; Nakama et al. 2017; Pattison et al. 2017; Garcia-Bellido et al. 2017; Franciolini et al. 2018; Cai et al. 2019; De Luca et al. 2019; Young et al. 2019; Iacconi & Mulryne 2023). In the parameter region where SIGWs manage to explain the NANOGrav signal, this would notably aggravate the PBH overproduction problem, which means that sizable non-Gaussianities could completely rule out the SIGW interpretation of the signal. In this work, we also neglect evolutionary effects on the PBH mass function, namely accretion and mergers, where accretion effects are expected to be small for sub-solar-mass PBHs (Ali-Ha¨ımoud & Kamionkowski 2017; Raidal et al. 2019; Vaskonen & Veermae 2020; De Luca et al. 2020).  \nIn summary, the final PBH DM fraction is quite sensitive to the assumptions discussed above. The bounds presented in the main text should be taken with a grain of salt owing to the uncertainties in the computation of f PBH .", 'C.3. Cosmological phase transitions': 'In the phase transition analysis discussed in the main text, we allow the low-frequency spectral index to float despite that causality predicts a = 3. We do this to capture possible double-peak spectral features with our simple power-law parametrization. However, it is not clear whether or not a strong and fast phase transition like the one needed to explain the observed signal could produce such a double-peak structure (Hindmarsh et al. 2017; Hindmarsh & Hijazi  \n∗  \n∗  \n∗  \n∗  \nFigure 25. Same as in Fig. 5 but for the pt-sound (left panel) and pt-bubble models with a low-frequency slope fixed to the value predicted by causality, i.e., a = 3.  \n<!-- image -->  \n<!-- image -->  \n∗  \n∗  \n∗  \n∗  \n2019). Therefore, in Fig. 25 we report the results of a phase transition analysis where we assume a = 3. Figure 25 shows the reconstructed posterior distributions for the parameters α ∗ , T ∗ , and H ∗ R ∗ of the pt-sound and pt-bubble models, both for the case where the PT is assumed to be the only source of GWs (blue contours) and for the scenario where we consider the superposition of the PT and SMBHBs signal (red contours). For the analyses where the SMBHB signal is included, we also report the posterior distributions for A BHB and γ BHB . For the pt-sound model we notice very minor differences compared to the analysis in the main text. However, for the pt-bubble model we notice how the posterior for T ∗ is peaked to slightly smaller values for the reasons explained in the main text.  \nIn Figs. 26 and 27, we report the posterior distributions for all the parameters of the phase transition models, including the spectral shape parameters a , b , and c that were excluded from Fig. 8. As expected for the pt-bubble model, the low-frequency slope is peaked around a ∼ 2, which is the reconstructed slope of the GWB signal, while the posteriors for b and c are approximately flat. For the pt-sound model, the posterior for a is peaked around the lower limit of the prior range a = 3, and there is also a mild preference for larger values of the width parameter, as this would flatten the spectrum close to the peak.', 'C.4. Cosmic strings': 'In our discussion of stable cosmic strings in the main text, we only present the reconstructed marginalized 1D posterior distributions for the dimensionless cosmic-string tension, Gµ (see Fig. 9). The four strings-plus-SMBHBs models stable-c + smbhb , stable-k + smbhb , stable-m + smbhb , and stable-n + smbhb , however, feature three model parameters in total: Gµ , A BHB , and γ BHB . In this appendix, we complement the discussion in Section 5.4 and show the corner plots for these parameters (see Fig. 28). A notable feature of these corner plots is that the plateau region at small values of Gµ , which we had already observed in Fig. 9, now also appears in the form of flat directions in the 2D posterior distributions for Gµ and A BHB as well as for Gµ and γ BHB . Meanwhile, the 2D posterior distributions for A BHB and γ BHB represent distorted versions of the 2D Gaussian prior distribution that we employ in our analysis, peaking at large A BHB and small γ BHB , where SMBHBs yield the dominant contribution to the signal. All four corner plots in Fig. 28 are qualitatively identical and only display slight numerical differences.  \nNext, let us turn to metastable strings. In the main text, we discuss the number density of closed string loops for the case of stable strings (see Eq. (43)) as well as for metastable strings (see Eq. (46)). In the meta-ls model, we need in addition the number density of string segments that form when long strings and closed loops break apart as a consequence of monopole nucleation. In order to compute the GW signal from segments, we use again Eqs. (40) and  \nFigure 26. Same as Fig. 8 but including the spectral shape parameters a , b , c and SMBHB parameters A BHB and γ BHB .  \n<!-- image -->  \n∗  \n∗  \n∗  \n∗  \n(42), where we replace the loop number density n l by the segment number density n s and the GW power spectrum P k in Eq. (41) by P k = 4 /k , which was derived by Martin and Vilenkin in the approximation of a straight string segment in Martin & Vilenkin (1997). Meanwhile, all relevant expressions for the segment number density n s were computed in Leblond et al. (2009) and Buchmuller et al. (2021), which we shall summarize in this appendix. For more details, we refer to (Leblond et al. 2009) and (Buchmuller et al. 2021).  \nFirst, we consider segments that form when long strings break apart. We denote their number density by ¯ n s , which we decompose into three different contributions that are relevant at different times and in different parameter regimes,  \n¯ n s ( ℓ, t ) = ¯ n rr s ( ℓ, t ) + ¯ n rm s ( ℓ, t ) + ¯ n mm s ( ℓ, t ) . (C15)  \n∗  \n∗  \n∗  \n∗  \nFigure 27. Same as Fig. 26 but for the pt-sound model.  \n<!-- image -->  \nHere ¯ n rr s ( ℓ, t ) is the number density of segments that originate from long strings, assuming the scaling regime to end during radiation domination, evaluated during radiation domination and after the onset of the decay regime,  \n¯ n rr s ( ℓ, t ) = Θ ( t eq -t ) Θ ( t -t s ) Γ 2 d ξ 2 r ( t + t s ) 2 ( t 3 t s ) 1 / 2 exp[ -Γ d ( ℓ ( t + t s ) + 1 / 2 Γ Gµ ( t -t s ) ( t +3 t s ))] , (C16)  \nwhere ¯ n rm s ( ℓ, t ) is the number density of segments that originate from long strings, assuming that the scaling regime ends during radiation domination, evaluated during matter domination and after the onset of the decay regime,  \n¯ n rm s ( ℓ, t ) = Θ ( t -t eq ) Θ ( t eq -t s ) Γ 2 d ξ 2 r ( t eq t ) 2 ( t + t s ) 2 ( t 3 eq t s ) 1 / 2 exp[ -Γ d ( ℓ ( t + t s ) + 1 / 2 Γ Gµ ( t -t s ) ( t +3 t s ))] , (C17)  \nFigure 28. Posterior distributions for the stable+smbhb strings models including the SMBHB parameters A BHB and γ BHB . We also report the marginalized log 10 Gµ posterior distributions for the stable string models (blue lines).  \n<!-- image -->  \nand ¯ n mm s ( ℓ, t ) is the number density of segments that originate from long strings, assuming that the scaling regime ends during matter domination, evaluated during matter domination and after the onset of the decay regime,  \n¯ n mm s ( ℓ, t ) = Θ ( t -t s ) Θ ( t s -t eq ) Γ 2 d ξ 2 m exp[ -Γ d ( ℓ t + 1 / 2 Γ Gµ ( t -t s ) ( t + t s ))] . (C18)  \nIn view of these expressions, several comments are in order: (i) Throughout our analysis, we assume that GW emission by segments is as efficient as GW emission by loops, i.e., we work with Γ = 50 for both loops and segments. (ii) All three expressions depend on the dimensionless correlation length of the long-string network, for which we use the attractor values in the VOS model, ξ r = 0 . 27 and ξ m = 0 . 63 during radiation and matter domination, respectively. (iii) For a given choice of parameter values, ¯ n rr s , ¯ n rm s , and ¯ n mm s never contribute simultaneously to the segment number density. For t eq > t s , only ¯ n rr s and ¯ n rm s yield nonvanishing contributions (first ¯ n rr s at t < t eq and then ¯ n rm s at t > t eq ), whereas for t eq < t s , the number density of segments from long strings is solely determined by ¯ n mm s at all times t > t s .  \n√  \n√  \nFigure 29. Same as Fig. 10 but including the SMBHB parameters A BHB and γ BHB .  \n<!-- image -->  \nFigure 30. Same as Fig. 11 but including the SMBHB parameters A BHB and γ BHB .  \n<!-- image -->  \nIn addition to segments from long strings, there is also a population of segments that form when closed string loops begin to break apart because of monopole nucleation. We denote the number density of this population by ˜ n s, 1 , where the index refers to the fact that ˜ n s, 1 only describes the first generation of segments that form when closed loops break apart for the first time. This first generation then gives rise to a second generation of segments that follow from monopole nucleation on first-generation segments, and so on and so forth. We comment on these higher generations further below. Before that, however, we discuss ˜ n s, 1 , which we decompose again into three contributions,  \n˜ n s, 1 ( ℓ, t ) = ˜ n rr s, 1 ( ℓ, t ) + ˜ n rm s, 1 ( ℓ, t ) + ˜ n mm s, 1 ( ℓ, t ) . (C19)  \nHere ˜ n rr s, 1 ( ℓ, t ) is the number density of first-generation segments that originate from closed loops that formed during radiation domination, evaluated during radiation domination and after the onset of the decay regime,  \n˜ n rr s, 1 ( ℓ, t ) = Θ ( t -t s ) Θ ( t eq -t ) Θ ( t eq -t ∗ ) Γ d [ ℓ ( t -t s ) + 1 / 2 Γ Gµ ( t -t s ) 2 ] n meta l ( ℓ, t ) , (C20)  \nwhere ˜ n rm s, 1 ( ℓ, t ) is the number density of first-generation segments that originate from closed loops that formed during radiation domination, evaluated during matter domination and after the onset of the decay regime,  \n˜ n rm s, 1 ( ℓ, t ) = Θ ( t -t s ) Θ ( t -t eq ) Θ ( t eq -t ∗ ) Γ d [ ℓ ( t -t s ) + 1 / 2 Γ Gµ ( t -t s ) 2 ] n meta l ( ℓ, t ) , (C21)  \nand ˜ n mm s, 1 ( ℓ, t ) is the number density of first-generation segments that originate from closed loops that formed during matter domination, evaluated during matter domination and after the onset of the decay regime,  \n˜ n mm s, 1 ( ℓ, t ) = Θ ( t -t s ) Θ ( t s -t ∗ ) Θ ( t ∗ -t eq ) Γ d e -Γ d [ ℓ ( t -t ∗ )+ 1 / 2 Γ Gµ ( t -t ∗ ) 2 ] t 2 ( ℓ +Γ Gµt ) 2 (C22)  \n× { 0 . 27 [ ℓ ( t -t s ) + 1 / 2 Γ Gµ ( t -t s ) 2 ] +0 . 45 ( ℓ +Γ Gµt ) 1+0 . 31 [ F 2 ( t ) -F 1 ( t ) -F 2 ( t s ) + F 1 ( t s )] } , (C23)  \nwhere F 2 and F 1 are shorthand notations for the following expressions involving the hypergeometric function 2 F 1 ,  \nF n ( x ) = 2 F 1 ( n -β, -β ; n +1 -β ; Γ Gµ ℓ +Γ Gµt x )( Γ Gµ ℓ +Γ Gµt ) n -1 x n -β n -β , n = 1 , 2 , β = 0 . 31 . (C24)  \nThese results were derived in Buchmuller et al. (2021) based on the loop number densities in the so-called BOS model (Blanco-Pillado et al. 2014), which explains the occurrence of factors such as 0 . 27, 0 . 45, and 0 . 31 in Eqs. (C22) and (C23) (see Buchmuller et al. (2021) for more details).  \nFinally, the expressions in Eqs. (C20), (C21), and (C22) enable us to estimate the total number densities of all generations of segments that descend from closed loops breaking apart because of monopole nucleation. In principle, these number densities are described by a partial integro-differential equation where the abundance of the n th segment generation acts as a source term for the ( n +1)th generation. Formally, this partial integro-differential equation can be solved analytically in terms of an infinite recursive series (Buchmuller et al. 2021). The numerical evaluation of this series is, however, technically demanding, which is why we choose to follow a different approach for the purposes of this paper. As shown in Buchmuller et al. (2021), it turns out that the total number densities for segments from closed loops result in predictions for the GW spectrum that are very similar to those obtained from the corresponding first-generation number densities in Eqs. (C20), (C21), and (C22) - up to a numerical fudge factor of O (10). At the level of the number densities, it therefore suffices to rescale all first-generation number densities by a constant factor in order to obtain effective number densities for all generations of segments that descend from closed loops breaking apart,  \n˜ n rr s ( ℓ, t ) → fudge × ˜ n rr s, 1 ( ℓ, t ) , ˜ n rm s ( ℓ, t ) → fudge × ˜ n rm s, 1 ( ℓ, t ) , ˜ n mm s ( ℓ, t ) → fudge × ˜ n mm s, 1 ( ℓ, t ) . (C25)  \nWe stress that the functional dependence on ℓ and t is typically not the same for the first-generation and total number densities; the rescaling in Eq. (C25), however, achieves comparable results at the level of the GW spectrum. In our analysis, we consistently use a fudge factor of 5, which is a characteristic value across large regions of parameter space.', 'C.5. Domain walls': 'In Fig. 31, we report the posterior distribution for all the parameters of the domain wall models, including the spectral shape parameters b and c that were excluded from Fig. 12. For both the dw-sm and dw-dr models we find a flat posterior for the high-frequency slope, b . This result is expected since most of the low frequency bins are fit by the low-frequency tail of the spectrum. For the width parameter, c , instead, both models prefer values near the upper range of our prior range. This preference for wider spectra derives from a mismatch between the reconstructed slope of the GWB, a ∼ 2 and the one predicted by causality, a ∼ 3. Larger width parameters make the spectrum from domain walls flatter near the peak and allow for a better fit to the data.', 'C.6. Dark matter substructures': "In this appendix we provide more details on the procedure that we followed to derive the constraints on f PBH reported in Section 6.2.  \nGiven a PBH with position relative to the pulsar given by r ( t ) = r 0 + vt , where r 0 and v are the initial PBH position and velocity, respectively, we can write the Doppler and Shapiro signals as  \nh D ( t ) = GM v 2 ˆ d · ( √ 1 + x 2 D ˆ b D -sinh -1 ( x )ˆ v ) (C26)  \nh S ( t ) = 2 GM log(1 + x 2 S ) , (C27)  \nwhere x D ≡ ( t -t D,0 ) /τ D and x S ≡ ( t -t S,0 ) /τ S . These expressions only include cubic terms in time t and have been previously derived in (Dror et al. 2019; Lee et al. 2021b). For the static limit in which τ ≫ T obs , these expressions can be expanded in series of t 0 /τ . Since the O ( t 2 ) terms would be degenerate with the timing model, the measurable new-physics signal can then be parametrized as a term ∝ t 3 as  \nh ( t ) D(S) = A D (S), sta yr 2 t 3 , (C28)  \nwhere A D(S), sta for both the Doppler and Shapiro static signal cases are dimensionless amplitudes given by  \nFigure 31. Same as Fig. 12 but including the spectral shape parameters b , c and the SMBHB parameters A BHB and γ BHB .  \n<!-- image -->  \n∗  \nA D, stat = yr 2 GM 2 v 2 ˆ d · [ t D, 0 τ 4 D 1 (1 + t 2 D, 0 /τ 2 D ) 5 / 2 ˆ b D + 1 3 τ 3 D 1 -2 t 2 D, 0 /τ 2 D (1 + t 2 D, 0 /τ 2 D ) 5 / 2 ˆ v ] (C29)  \nA S, stat = -yr 2 4 GM 3 t S, 0 τ 4 S 3 -t 2 S, 0 /τ 2 S (1 + t 2 S, 0 /τ 2 S ) 3 . (C30)  \nFor the Doppler case, in the dynamic limit when τ ≪ T obs , the dominant contribution would come from the first term in Eq. (C26) where √ 1 + x 2 D ∝ | t -t 0 | . Upto linear order in x D , we can write  \nh D, dyn ( t ) = A D, dyn ( t -t 0 )Θ( t -t 0 ) , (C31)  \nwhere A D , dyn is the dimensionless amplitude given by  \nA D, dyn = 2 GM τv 2 ˆ d · ˆ b D . (C32)  \nGiven the expressions in Eqs. (C26) and (C27) for the Doppler and Shapiro signals, we use the MC developed in Lee et al. (2021a) to derive the distributions p (log 10 A I | f PBH ). Specifically, we proceed as follows:  \n- · For each pulsar we generate a population of N PBH PBHs, where N PBH is implicitly defined by the relation  \nN PBH = f PBH ρ DM V M PBH , (C33)  \nwhere the simulation volume, V , is a sphere of radius R = ¯ vT obs centered around the pulsar position for the Doppler signal and a cylinder with the same radius and height given by the Earth-pulsar distance for the Shapiro signal. Here ¯ v ≃ 340 km / s is the average PBH velocity, and T obs ,I is the observation time of the I th pulsar.  \n- · For each of these PBHs we generate a random initial position and velocity. Since PTA searches are only sensitive to DM subhalos in the neighborhood of the solar system, we expect the position distribution to be uniform. Therefore, we use the probability density function f ( r ) = 1 / V to sample initial positions. To sample PBHs' velocity, we use a Maxwell-Boltzmann distribution with v 0 = 325kms -1 , v esc = 600kms -1 , and the angular dependence assumed to be isotropic.\n- · The PBHs' signals are then classified as dynamic if they satisfy the condition T obs ,I -τ > t 0 > τ , and static otherwise.\n- · To evaluate A stat , we sum the static signals of all PBHs computed by using Eqs. (C26) and (C27), and we fit the resulting signal to a cubic polynomial in time to extract the t 3 term. To compute A D , dyn , we take the closest transiting object and compute the signal amplitudes using Eq. (C32).  \nAll the previous points are repeated for 2 . 5 × 10 3 realizations to obtain the distributions p (log 10 A I | f PBH ). Given the conditional distributions p (log 10 A I | f PBH ) and the posterior distribution p (log 10 A I | δt ) derived by analyzing our data, we can write  \np ( f PBH | δt ) = N P ∏ I =1 ∫ p ( f PBH | log 10 A I ) p (log 10 A I | δt ) d log 10 A I . (C34)  \nThen, using Bayes' theorem, we can rewrite  \np ( f PBH | log 10 A I ) = p (log 10 A I | f PBH ) p ( f PBH ) p (log 10 A I ) . (C35)  \nOur priors on p ( f PBH ) and p (log 10 A I ) are uniform, which allows us to eventually write Eq. (C34) as  \np ( f PBH | δt ) ∝ N P ∏ I =1 ∫ p (log 10 A I | f PBH ) p (log 10 A I | δt ) d log 10 A I , (C36)  \nwhere the ∝ implies that the p ( f PBH | δ t ) would be subject to the normalization condition, ∫ ∞ 0 p ( f PBH | δt ) df PBH = 1.  \nIn the presence of a DM-baryon fifth force in the form of a Yukawa potential in Eq. (72), assuming point-mass DM, the pulsar frequency shift due to Doppler effect is given by (Gresham et al. 2023)  \n( δν ν ) fifth = ˜ αGM τ 2 D v ˆ d · ∫ 1 (1 + x 2 D ) 3 / 2 ( 1 + b λ √ 1 + x 2 D ) e -( b/λ ) √ 1+ x 2 D ( ˆ b + x D ˆ v ) dx D . (C37)  \nThe integral in Eq. (C37) has to be computed numerically, and the signal due to the fifth force can be computed by performing an additional integration over time and subtracting away degeneracies with timing model parameters. The total signal is the sum of the gravitational and the fifth-force contribution, h D, total ( t ) = h D, fifth ( t ) + h D ( S ) ( t ). In this analysis of the fifth-foce constraint, we only consider the scenario where the DM substructure makes up the entirety of the DM local density, which is equivalent to taking f PBH = 1 for the gravitational contribution. Parameterizing the signal as h D , total ( t ) = A D, total yr 2 t 3 similar to the PBH case, the amplitude A D, total is a random variable dependent on λ and ˜ α . The probability distribution function P (log 10 A D, total | λ, ˜ α ) can be extracted again by Monte Carlo simulations and Bayes' theorem. Finally, the posterior distribution of ˜ α and λ , P (˜ α, λ | δt ), is given by the analog of Eq. (C36)  \np (˜ α, λ | δt ) ∝ N P ∏ I =1 ∫ p (log 10 A I | ˜ α, λ ) p (log 10 A I | δt ) d log 10 A I . (C38)"}