Corresponding author: The NANOGrav Collaboration
comments@nanograv.org {sup}21Department of Physics, New Mexico Institute of Mining and Technology, 801 Leroy Place, Socorro, NM 87801, USA
{sup}22Department of Physics, Montana State University, Bozeman, MT 59717, USA
{sup}23Department of Physics and Astronomy, Franklin & Marshall College, P.O. Box 3003, Lancaster, PA 17604, USA
{sup}24Department of Physics and Astronomy, University of British Columbia, 6224 Agricultural Road, Vancouver, BC V6T 1Z1, Canada
{sup}25Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, USA
{sup}26George Mason University, resident at the Naval Research Laboratory, Washington, DC 20375, USA
{sup}27National Radio Astronomy Observatory, 1003 Lopezville Rd., Socorro, NM 87801, USA
{sup}28Department of Physics, Hillsdale College, 33 E. College Street, Hillsdale, MI 49242, USA
{sup}29Eureka Scientific, 2452 Delmer Street, Suite 100, Oakland, CA 94602-3017, USA
{sup}30School of Physics and Astronomy, Rochester Institute of Technology, Rochester, NY 14623, USA
{sup}31Laboratory for Multiwavelength Astrophysics, Rochester Institute of Technology, Rochester, NY 14623, USA
{sup}32Institute for Theoretical Physics, University of Münster, 48149 Münster, Germany
{sup}33Department of Astronomy, University of Maryland, College Park, MD 20742, USA
{sup}34Center for Research and Exploration in Space Science and Technology, NASA/GSFC, Greenbelt, MD 20771, USA
{sup}35NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA
{sup}36Department of Physics, University of Connecticut, 196 Auditorium Road, U-3046, Storrs, CT 06269-3046, USA
{sup}37Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue, New York, NY 10010, USA
{sup}38Department of Physics and Astronomy, Haverford College, Haverford, PA 19041, USA
{sup}39Department of Astronomy and Astrophysics, University of Michigan, Ann Arbor, MI 48109, USA
{sup}40Department of Astronomy, University of California, Berkeley, 501 Campbell Hall #3411, Berkeley, CA 94720, USA
{sup}41Space Science Division, Naval Research Laboratory, Washington, DC 20375-5352, USA
{sup}42University of Washington Bothell, 18115 Campus Way NE, Bothell, WA 98011, USA
{sup}43Department of Physics, State University of New York at Oswego, Oswego, NY 13126, USA
{sup}44CP3-Origins, University of Southern Denmark, Campusvej 55, DK-5230 Odense M, Denmark
{sup}45NASA Marshall Space Flight Center, Huntsville, AL 35812, USA
{sup}46Department of Astronomy & Astrophysics, University of Toronto, 50 Saint George Street, Toronto, ON M5S 3H4, Canada
{sup}47Green Bank Observatory, P.O. Box 2, Green Bank, WV 24944, USA
{sup}48Department of Physics, University of California, Berkeley, CA 94720, USA
{sup}49Department of Physics, University of the Pacific, 3601 Pacific Avenue, Stockton, CA 95211, USA
{sup}50E.A. Milne Centre for Astrophysics, University of Hull, Cottingham Road, Kingston-upon-Hull, HU6 7RX, UK
{sup}51Centre of Excellence for Data Science, Artificial Intelligence and Modelling (DAIM), University of Hull, Cottingham Road, Kingston-upon-Hull, HU6 7RX, UK
{sup}52International Centre for Radio Astronomy Research, Curtin University, Bentley, WA 6102, Australia
{sup}53Department of Physics, Yale University, New Haven, CT 06520, USA
{sup}54Deutsches Elektronen-Synchrotron DESY, Notkestr. 85, D-22607 Hamburg, Germany
{sup}55Department of Astronomy, Yale University, 52 Hillhouse Ave., New Haven, CT 06511, USA
{sup}56Black Hole Initiative, Harvard University, 20 Garden Street, Cambridge, MA 02138, USA
{sup}57Dunlap Institute for Astronomy and Astrophysics, University of Toronto, 50 St. George St., Toronto, ON M5S 3H4, Canada
{sup}58Department of Physics, Lafayette College, Easton, PA 18042, USA
{sup}59Institute of Cosmology, Department of Physics and Astronomy, Tufts University, Medford, MA 02155, USA
{sup}60Institute of Physics and Astronomy, Eötvös Loránd University, Pázmány P. s. 1/A, 1117 Budapest, Hungary
{sup}61Arecibo Observatory, HC3 Box 53995, Arecibo, PR 00612, USA
{sup}62Department of Physics, University of Puerto Rico, Mayagüez, PR 00681, USA
{sup}63National Radio Astronomy Observatory, 520 Edgemont Road, Charlottesville, VA 22903, USA
{sup}64Department of Physics, Texas Tech University, Box 41051, Lubbock, TX 79409, USA
{sup}65Department of Physics, Penn State Abington, Abington, PA 19001, USA
{sup}66Giant Army, 915A 17th Ave., Seattle, WA 98122, USA
{sup}67Department of Astrophysical and Planetary Sciences, University of Colorado, Boulder, CO 80309, USA
{sup}68Center for Astrophysics, Harvard University, 60 Garden St., Cambridge, MA 02138, USA
{sup}69Department of Physics and Astronomy, Oberlin College, Oberlin, OH 44074, USA
{sup}70Theoretical Physics Division, Fermi National Accelerator Laboratory, Batavia, IL 60510, USA
{sup}71Department of Physics, Ben-Gurion University of the Negev, Be'er Sheva 84105, Israel
{sup}72Feza Gürsey Institute, Boğaziçi University, Kandilli, 34684, Istanbul, Turkey
{sup}73Scuola Normale Superiore, Piazza dei Cavalieri 7, I-56100 Pisa, Italy
{sup}74Deutsches Elektronen-Synchrotron DESY, Notkestr. 85, 22607 Hamburg, Germany
{sup}75Center for Interdisciplinary Exploration and Research in Astrophysics (CIERA), Northwestern University, Evanston, IL 60208, USA
{sup}76Adler Planetarium, 1300 S. DuSable Lake Shore Dr., Chicago, IL 60605, USA
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 [219 , 70 , 83 ] . These GWs form a stochastic background and propagate essentially unimpeded over cosmic distances to be detected today, whereas electromagnetic radiation does not start free 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.
Cosmological GWBs can be produced by a number of particle physics models of the early universe. Notably, cosmic inflation generically produces GWs [146 ] , which may be observable at nanohertz frequencies if their energy density spectrum is sufficiently blue-tilted. Similarly, an enhanced spectrum of short-wavelength scalar perturbations produced during inflation can source so-called scalar-induced GWs (SIGWs; [100 , 318 ] ). Another potential source of GWs are cosmological first-order phase transitions [71 , 68 , 158 ] , 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 [157 , 267 ] (Vilenkin 1985), can radiate GWs and hence contribute to the GWB.
The North American Nanohertz Observatory for Gravitational Waves (NANOGrav; [228 ] ) has recently found the first convincing evidence for a stochastic GWB signal, as detailed in [13 hereafter NG15gwb] . Analyzing 15-year of pulsar timing observations, NANOGrav has detected a red-noise 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; [258 , 173 , 307 , 274 , 64 ] ), cosmological sources may also contribute to it.
The SMBHB interpretation of the signal is considered in Agazie et al. (2023c, hereafter NG15smbhb). In this paper, we analyze the NANOGrav 15-year data set [12 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 strings (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.
The 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.
In 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 K K K K K
Aside 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 signals. Specifically, DM substructures within the Milky Way can produce a Doppler effect by accelerating the Earth or a pulsar ([276 ] ), or a Shapiro delay of the photons’ arrival times by perturbing the metric along the photon geodesic ([278 ] ). PTAs can also probe models of ultralight DM (ULDM), which can cause shifts in the observed pulse timing via metric fluctuations ([187 , 251 ] ) or via couplings between ULDM and SM particles ([139 , 181 ] ). 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.
This paper is organized as follows. We describe the NG15 data set in § 2 and our general analysis methods in § 3 . In § 4 , we discuss the GWB expected from SMBHBs. We present the analysis and results for new-physics models that generate a cosmological GWB in § 5 and for models that produce deterministic signals in § 11 . We conclude in § 14 . 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 § 16 . We present median GW spectra for all cosmological models based on our recovered posterior distributions in Appendix § 17 , and we provide supplementary material for specific models in Appendix § 18 .
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 ([20 ] ). 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).
The pulse times of arrival (TOAs) were generated from the raw data following the procedure discussed in [36 , 35 , 20 ] . 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 F F [36 , 177 ] . The timing model fits were performed using the TT(BIPM2019) timescale and the JPL Solar System Ephemeris model DE440 [245 ] . Additional detail about the data set and the processing of the TOAs can be found in NG15 and Agazie et al. (2023d), hereafter NG15detchar.
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., [35 ] (Arzoumanian et al. 2016)). In the following brief overview, we focus on the implementation of new-physics signals within this framework.
Likelihood
Our search for a new-physics signal utilizes the pulsars’ timing residuals, δ t \delta t δ t § 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 [292 ] . Specifically, we can model the timing residuals as
δ t = n + F a + M ϵ . (1) \delta t = \mathbf{n} + \mathbf{F}\mathbf{a} + \mathbf{M}\boldsymbol{\epsilon}. \tag*{(1)} δ t = n + Fa + M ϵ . (1) In the remainder of this section, we will define and discuss each of these three terms and define the PTA likelihood.
The first term on the right-hand side of Eq. (1) , n \mathbf{n} n N T O A N_{\mathrm{TOA}} N TOA I I I
⟨ n i n j ⟩ = F I 2 [ σ i S / N 2 + Q I 2 ] δ i j + J I 2 U i j , (2) \langle n_i n_j \rangle= \mathcal{F}_I^2 \left[\sigma_{i\ \mathrm{S/N}}^2 + Q_I^2\right]\delta_{ij} + \mathcal{J}_I^2 U_{ij}, \tag*{(2)} ⟨ n i n j ⟩ = F I 2 [ σ i S/N 2 + Q I 2 ] δ ij + J I 2 U ij , (2) where i i i j j j σ i S / N \sigma_{i\ \mathrm{S/N}} σ i S/N i i i F I \mathcal{F}_I F I Extra FACtor (EFAC) parameter, Q I Q_I Q I Extra QUADrature (EQUAD) parameter, and J I \mathcal{J}_I J I U \mathcal{U} U [35 ] (Arzoumanian et al. 2016), we fix all white-noise parameters to their values at the maxima in the posterior probability distributions recovered from single pulsar noise studies in order to increase computational efficiency (NG15detchar).
Time-correlated stochastic processes, like pulsar-intrinsic red noise and GWB signals, are modeled using a Fourier basis of frequencies i / T o b s i/T_{\mathrm{obs}} i / T obs i i i T o b s T_{\mathrm{obs}} T obs N f N_f N f N f = 30 N_f = 30 N f = 30 N f = 14 N_f = 14 N f = 14 2 N f 2N_f 2 N f F \mathbf{F} F a \mathbf{a} a ⟨ a a T ⟩ = φ \langle\mathbf{a}\mathbf{a}^{\mathsf{T}}\rangle= \boldsymbol{\varphi} ⟨ a a T ⟩ = φ
[ φ ] ( a k ) ( b i ) = δ i j ( Γ a b Φ i + δ a b φ a , i ) (3) \left[\boldsymbol{\varphi}\right]_{(ak)(bi)} = \delta_{ij}\left(\Gamma_{ab}\Phi_i + \delta_{ab}\varphi_{a,i}\right) \tag*{(3)} [ φ ] ( ak ) ( bi ) = δ ij ( Γ ab Φ i + δ ab φ a , i ) (3) where a a a b b b i i i j j j Γ a b \Gamma_{ab} Γ ab a a a b b b Γ a b \Gamma_{ab} Γ ab [150 ] , also known as “quadrupolar” or “HD” correlation.
The 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 \Phi_i Φ i § 4 ), and one of cosmological origin, induced by one of the exotic new-physics models under consideration (discussed in § 5 ). The last term in Eq. (3) models pulsar-intrinsic red-noise in terms of the coefficients φ a , i \varphi_{a,i} φ a , i
φ a ( f ) = A a 2 12 π 2 T o b s 1 1 y r − 1 ( f 1 y r − 1 ) − γ a y r 3 (4) \varphi_a(f) = \frac{A_a^2}{12\pi^2 T_{\mathrm{obs}}}\frac{1}{1\ \mathrm{yr}^{-1}}\left(\frac{f}{1\ \mathrm{yr}^{-1}}\right)^{-\gamma_a}\mathrm{yr}^3 \tag*{(4)} φ a ( f ) = 12 π 2 T obs A a 2 1 yr − 1 1 ( 1 yr − 1 f ) − γ a yr 3 (4) and φ a , i = φ a ( i / T o b s ) \varphi_{a,i} = \varphi_a(i/T_{\mathrm{obs}}) φ a , i = φ a ( i / T obs ) N f N_f N f Table 2 . Finally, deviations from the initial best-fit values of the m m m M ϵ M\epsilon M ϵ design matrix , M M M N T O A × m N_{\mathrm{TOA}} \times m N TOA × m ϵ \epsilon ϵ
Table 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 B H B , γ B H B ) (\log_{10} A_{\mathrm{BHB}}, \gamma_{\mathrm{BHB}}) ( log 10 A BHB , γ BHB ) (73) .
Since 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 a a ϵ \epsilon ϵ A a A_a A a γ a \gamma_a γ a Φ i \Phi_i Φ i [293 , 211 ] ):
p ( δ t ∣ ϕ ) = exp ( − 1 2 δ t T C − 1 δ t ) det ( 2 π C ) (5) p(\delta t|\phi) = \frac{\exp\left(-\frac{1}{2}\delta t^{T}C^{-1}\delta t\right)}{\sqrt{\det(2\pi C)}} \tag*{(5)} p ( δ t ∣ ϕ ) = det ( 2 π C ) exp ( − 2 1 δ t T C − 1 δ t ) (5) where C = N + T B T T C = N + TBT^{T} C = N + T B T T N N N T = [ M , F ] T = [M,F] T = [ M , F ] B = diag ( ∞ , ϕ ) B = \operatorname{diag}(\infty,\phi) B = diag ( ∞ , ϕ ) ∞ \infty ∞ ϵ \epsilon ϵ B B B
Eq. (5) can be easily generalized to take into account deterministic signals (like the ones that will be discussed in § 11.1 and § 12.1 ). In the presence of a deterministic signal, h ( θ ) h(\theta) h ( θ ) θ \theta θ δ t → δ t − h ( θ ) \delta t \to\delta t - h(\theta) δ t → δ t − h ( θ )
Finally, we relate our characterization of the GWB given in Eq. (3) in terms of Φ i \Phi_i Φ i d ρ G W / d ln f d\rho_{\mathrm{GW}}/d\ln f d ρ GW / d ln f ρ c \rho_c ρ c [23 ] )
Ω G W ( f ) ≡ 1 ρ c d ρ G W ( f ) d ln f = 8 π 4 f 5 H 0 2 Φ ( f ) Δ f (6) \Omega_{\mathrm{GW}}(f) \equiv\frac{1}{\rho_c}\frac{d\rho_{\mathrm{GW}}(f)}{d\ln f} = \frac{8\pi^4 f^5}{H_0^2}\frac{\Phi(f)}{\Delta f} \tag*{(6)} Ω GW ( f ) ≡ ρ c 1 d ln f d ρ GW ( f ) = H 0 2 8 π 4 f 5 Δ f Φ ( f ) (6) Here H 0 H_0 H 0 Δ f = 1 / T o b s \Delta f = 1/T_{\mathrm{obs}} Δ f = 1/ T obs N f N_f N f Φ ( f ) \Phi(f) Φ ( f ) Φ i \Phi_i Φ i (3) , i.e., Φ i = Φ ( i / T o b s ) \Phi_i = \Phi(i/T_{\mathrm{obs}}) Φ i = Φ ( i / T obs ) Φ ( f ) \Phi(f) Φ ( f ) S ( f ) = Φ ( f ) / Δ f S(f) = \Phi(f)/\Delta f S ( f ) = Φ ( f ) /Δ f 1 / Δ f 1/\Delta f 1/Δ f h 2 Ω G W h^2\Omega_{\mathrm{GW}} h 2 Ω GW Ω G W \Omega_{\mathrm{GW}} Ω GW h h h H 0 = h × 100 k m s − 1 M p c − 1 H_0 = h \times100\ {\rm km\,s^{-1}Mpc^{-1}} H 0 = h × 100 km s − 1 Mp c − 1 H 0 H_0 H 0 h 2 Ω G W h^2\Omega_{\mathrm{GW}} h 2 Ω GW
Bayesian analysis
The goal of this work is to investigate a series of cosmological interpretations of the GWB signal in our data. 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?
To 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 \mathcal{H} H Θ \Theta Θ D \mathcal{D} D
P ( Θ ∣ D , H ) = P ( D ∣ Θ , H ) P ( Θ ∣ H ) P ( D ∣ H ) (7) P(\Theta|\mathcal{D},\mathcal{H}) = \frac{P(\mathcal{D}|\Theta,\mathcal{H})P(\Theta|\mathcal{H})}{P(\mathcal{D}|\mathcal{H})} \tag*{(7)} P ( Θ∣ D , H ) = P ( D ∣ H ) P ( D ∣Θ , H ) P ( Θ∣ H ) (7) where P ( Θ ∣ D , H ) P(\Theta|\mathcal{D},\mathcal{H}) P ( Θ∣ D , H ) P ( D ∣ Θ , H ) P(\mathcal{D}|\Theta,\mathcal{H}) P ( D ∣Θ , H ) P ( Θ ∣ H ) P(\Theta|\mathcal{H}) P ( Θ∣ H )
Z ≡ P ( D ∣ H ) = ∫ d Θ P ( D ∣ Θ , H ) P ( Θ ∣ H ) (8) \mathcal{Z} \equiv P(\mathcal{D}|\mathcal{H}) = \int d\Theta\,P(\mathcal{D}|\Theta,\mathcal{H})P(\Theta|\mathcal{H}) \tag*{(8)} Z ≡ P ( D ∣ H ) = ∫ d Θ P ( D ∣Θ , H ) P ( Θ∣ H ) (8) is the marginalized likelihood, or evidence. In the context of this work, H \mathcal{H} H (1) , Θ \Theta Θ Φ \Phi Φ δ t \delta t δ t (5) and implemented using the ENTERPRISE (Ellis et al. 2019) and ENTERPRISE_EXTENSIONS ([283 ] ) packages. Our prior choices are summarized in Table 2 and Table 3 .
Table 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.
The 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 \mathcal{H}_0 H 0 H 1 \mathcal{H}_1 H 1
B 10 ( D ) = Z 1 Z 0 = P ( D ∣ H 1 ) P ( D ∣ H 0 ) (9) \mathcal{B}_{10}(\mathcal{D}) = \frac{\mathcal{Z}_1}{\mathcal{Z}_0} = \frac{P(\mathcal{D}|\mathcal{H}_1)}{P(\mathcal{D}|\mathcal{H}_0)} \tag*{(9)} B 10 ( D ) = Z 0 Z 1 = P ( D ∣ H 0 ) P ( D ∣ H 1 ) (9) The 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 \mathcal{H}_1 H 1 [175 ] ): B 10 < 1 \mathcal{B}_{10} < 1 B 10 < 1 H 1 \mathcal{H}_1 H 1 B 10 \mathcal{B}_{10} B 10 [ 10 0.0 , 10 0.5 ] [10^{0.0},10^{0.5}] [ 1 0 0.0 , 1 0 0.5 ] [ 10 0.5 , 10 1.0 ] [10^{0.5},10^{1.0}] [ 1 0 0.5 , 1 0 1.0 ] [ 10 1.0 , 10 1.5 ] [10^{1.0},10^{1.5}] [ 1 0 1.0 , 1 0 1.5 ] [ 10 1.5 , 10 2.0 ] [10^{1.5},10^{2.0}] [ 1 0 1.5 , 1 0 2.0 ] [ 10 2.0 , ∞ ) [10^{2.0},\infty) [ 1 0 2.0 , ∞ ) H 1 \mathcal{H}_1 H 1
Given 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 [108 ] ) 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 [214 ] ).
In order to compute the Bayes factor between two models, we use product space methods (Carlin & Chib [72 ] ; Godsill [136 ] ; Hee et al. [148 ] ), instead of calculating the evidence Z Z Z P ( H 1 ) = P ( H 0 ) P(\mathcal{H}_1)=P(\mathcal{H}_0) P ( H 1 ) = P ( H 0 ) [103 ] ). 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 5 \times10^4 5 × 1 0 4
From 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 ) P(D|\Theta,\mathcal{H}) P ( D ∣Θ , H )
For each model H \mathcal{H} H P ( Θ ∣ D , H ) P(\Theta|\mathcal{D},\mathcal{H}) P ( Θ∣ D , H )
identify relevant parameter ranges and set upper limits. Specifically, we identify 68% (95%) Bayesian credible intervals (Bernardo & Smith [52 ] ) 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 [175 ] ). For instance, to place a bound on a single parameter θ \theta θ
K ( θ ) = P ( D ∣ θ , H ) P ( D ∣ θ 0 , H ) , (10) K(\theta)=\frac{P(D|\theta,\mathcal{H})}{P(D|\theta_0,\mathcal{H})}, \tag*{(10)} K ( θ ) = P ( D ∣ θ 0 , H ) P ( D ∣ θ , H ) , (10) has dropped to K = 1 / 10 K=1/10 K = 1/10 θ 0 \theta_0 θ 0 P ( D ∣ θ , H ) P(D|\theta,\mathcal{H}) P ( D ∣ θ , H ) θ \theta θ P ( D ∣ θ , H ) P(D|\theta,\mathcal{H}) P ( D ∣ θ , H ) θ \theta θ P ( D ∣ θ 0 , H ) P(D|\theta_0,\mathcal{H}) P ( D ∣ θ 0 , H ) θ \theta θ (10) is identical to the corresponding ratio of marginalized posteriors. Furthermore, multiplying and dividing by the prior on θ \theta θ
K ( θ ) = P ( θ ∣ H ) P ( θ 0 ∣ D , H ) P ( θ ∣ D , H ) P ( θ ∣ H ) . (11) K(\theta)=\frac{P(\theta|\mathcal{H})}{P(\theta_0|\mathcal{D},\mathcal{H})}\frac{P(\theta|\mathcal{D},\mathcal{H})}{P(\theta|\mathcal{H})}. \tag*{(11)} K ( θ ) = P ( θ 0 ∣ D , H ) P ( θ ∣ H ) P ( θ ∣ H ) P ( θ ∣ D , H ) . (11) The 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 ) \mathcal{B}=P(D|\mathcal{H})/P(D|\mathcal{H}_0) B = P ( D ∣ H ) / P ( D ∣ H 0 ) H 0 \mathcal{H}_0 H 0 H \mathcal{H} H θ \theta θ K K K θ \theta θ
K ( θ ) = B P ( θ ∣ D , H ) P ( θ ∣ H ) . (12) K(\theta)=\mathcal{B}\frac{P(\theta|D,\mathcal{H})}{P(\theta|\mathcal{H})}. \tag*{(12)} K ( θ ) = B P ( θ ∣ H ) P ( θ ∣ D , H ) . (12) Once B \mathcal{B} B (12) and determine the K K K θ \theta θ (12) is useful for numerically evaluating K K K P ( θ 0 ∣ D , H ) = P ( θ ∣ H ) / B P(\theta_0|D,\mathcal{H})=P(\theta|\mathcal{H})/\mathcal{B} P ( θ 0 ∣ D , H ) = P ( θ ∣ H ) / B K K K [42 ] ). For more than one parameter dimension, we proceed analogously and derive bounds based on the criterion K ( Θ ) > 1 / 10 K(\Theta)>1/10 K ( Θ ) > 1/10
All Bayesian inference analyses discussed in this work were implemented into ENTERPRISE via a newly developed wrapper that we call PTArcade ([230 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 .
GWB signal from SMBHBs
Most galaxies are expected to host a supermassive black hole (SMBH) at their center ([196 , 18 ] ). During the hierarchical merging of galaxies taking place in the course of structure formation ([303 ] ), these black holes are expected to sink to the center of the merger remnants, eventually forming binary systems ([47 ] ). The gravitational radiation emitted by this population of inspiraling SMBHBs forms a GWB in the PTA band ([258 , 173 , 307 ] ) and is a natural candidate for the source of the signal observed in our data.
The shape and normalization of this GWB depend on the properties of the SMBHB population and on its dynamical evolution ([111 , 275 , 192 , 186 ] ). 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 γ B H B = − 13 / 3 \gamma_{\mathrm{BHB}}=-13/3 γ BHB = − 13/3 [247 ] ), 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. to 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 ([47 ] ). 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 y r − 1 f \ll1\ \mathrm{yr}^{-1} f ≪ 1 yr − 1 [192 ] ). At high frequencies, once the expected number of binaries dominating the GWB approaches unity, the spectrum steepens below 13 / 3 13/3 13/3 f ≫ 1 y r − 1 f \gg1\ \mathrm{yr}^{-1} f ≫ 1 yr − 1 [275 ] ).
Unfortunately, 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 signal by a power law in this frequency range:
Φ B H B ( f ) = A B H B 2 12 π 2 1 T o b s ( f y r − 1 ) − γ B H B y r 3 (13) \Phi_{\rm BHB}(f) = \frac{A_{\rm BHB}^{2}}{12\pi^{2}}\frac{1}{T_{\rm obs}}\left(\frac{f}{\mathrm{yr}^{-1}}\right)^{-\gamma_{\rm BHB}}\mathrm{yr}^{3} \tag*{(13)} Φ BHB ( f ) = 12 π 2 A BHB 2 T obs 1 ( yr − 1 f ) − γ BHB yr 3 (13) where Φ B H B / Δ f \Phi_{\rm BHB}/\Delta f Φ BHB /Δ f (6) ).
Following Middleton et al. [229 ] , we can gain some insight into the allowed range of values for the amplitude, A B H B A_{\rm BHB} A BHB γ B H B \gamma_{\rm BHB} γ BHB [185 ] ) 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 parameters 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 ([47 ] ). 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 observations (e.g., Tomczak et al. [285 ] ) in addition to more detailed numerical studies of SMBHB evolution (e.g., Sesana [273 ] ).
For 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 B H B \mathcal{A}_{\mathrm{BHB}} A BHB γ B H B \gamma_{\mathrm{BHB}} γ BHB Figure 1 . (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 [273 ] ), 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 power-law parameterization with spectral index γ B H B = 13 / 3 \gamma_{\rm BHB} = 13/3 γ BHB = 13/3
Figure 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 \gamma=13/3 γ = 13/3 [247 ] ).
The black dashed contours in Figure 1 show the results of a 2D Gaussian fit to the distribution of A B H B A_{\rm BHB} A BHB γ B H B \gamma_{\rm BHB} γ BHB (73) in Appendix § 16 for the parameters of this Gaussian distribution). This fitted distribution is what we adopt as a prior distribution for A B H B A_{\rm BHB} A BHB γ B H B \gamma_{\rm BHB} γ BHB
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 SMBHB interpretation. In Figure 2 , we show a summary of these Bayes factors; in Figure 3 , we present median reconstructed GWB spectra in the PTA band for a number of select new-physics models; and in Figure 4 , we show similar median reconstructed GWB spectra in the broader landscape of present and future GW experiments.
Figure 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 § 3.2 . In most cases, however, these error bars are small and not visible.
Figure 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 Ω G W h^2\Omega_{\rm GW} h 2 Ω GW f f f § 17 for more details and Figure 19 and Figure 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 B H B A_{\rm BHB} A BHB γ B H B \gamma_{\rm BHB} γ BHB μ B H B \mu_{\rm BHB} μ BHB (73) (black dashed line).
Figure 4. Same as Figure 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 § 17 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 ([284 ] ) of various existing and planned GW interferometer experiments: LISA ([25 ] ), DECIGO ([183 ] ), BBO ([86 ] ), Einstein Telescope (ET; [256 ] ), Cosmic Explorer (CE; [261 ] ), the HLVK detector network (consisting of aLIGO in Hanford and Livingston ([1 ] ), aVirgo ([8 ] ), and KAGRA ([19 ] ) 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 ρ t h r \rho_{\mathrm{thr}} ρ thr t o b s t_{\mathrm{obs}} t obs ρ t h r / t o b s \rho_{\mathrm{thr}/\sqrt{t_{\mathrm{obs}}}} ρ thr / t obs [271 ] . 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 Ω G W h^2\Omega_{\mathrm{GW}} h 2 Ω GW Figure 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 ) ) \mathcal{O}(10^{-(5\cdots6)}) O ( 1 0 − ( 5 ⋯ 6 ) ) § 5.1 ).
As discussed in § 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 Ω G W h^2\Omega_{\rm GW} h 2 Ω GW
Cosmic inflation
MODEL DESCRIPTION
Cosmic inflation denotes a stage of exponential expansion in the early universe that provides an explanation for the initial conditions of Big Bang cosmology (Liddle & Lyth [217 ] ). At the level of the background expansion, inflation accounts for the size, homogeneity, isotropy, and flatness of the observable universe on cosmological scales; at the level of perturbations, it provides the seeds for structure formation in the form of primordial density fluctuations. In the standard scenario of inflation, these primordial perturbations are sourced by scalar quantum vacuum fluctuations of the spacetime metric and inflaton field, which are first stretched to superhorizon scales during inflation and then reenter the horizon in the form of classical density perturbations after inflation. In addition to scalar perturbations, inflation also leads to the amplification of tensor perturbations of the metric, which reenter the horizon in the form of stochastic GWs after inflation. These primordial or inflationary gravitational waves (IGWs) (Grishchuk [143 ] ; Starobinsky [280 ] ; Rubakov et al. [266 ] ; Fabbri & Pollock [119 ] ; Abbott & Wise [3 ] ) represent a prime GW signal from the early universe. For earlier work on the IGW interpretation of the signal in recent PTA data sets, see Vagnozzi [291 ] , Kuroyanagi et al. (2021), and Benetti et al. [48 ] .
IGWs leave an imprint in the temperature and polarization anisotropies of the cosmic microwave background (CMB) whose relative strength compared to the contributions from scalar perturbations is quantified in terms of the tensor-to-scalar ratio, r r r h 2 Ω G W h^2\Omega_{\rm GW} h 2 Ω GW n t n_t n t n t = − r / 8 < 0 n_t=-r/8<0 n t = − r /8 < 0 r r r r ≤ 0.036 r\leq0.036 r ≤ 0.036 [10 ] . A vanishing tensor spectral index, n t ≈ 0 n_t\approx0 n t ≈ 0 Ω G W h 2 ∼ 10 − 16 \Omega_{\rm GW}h^2\sim10^{-16} Ω GW h 2 ∼ 1 0 − 16
We remain agnostic about the microphysics of inflation and restrict ourselves to a model-independent analysis, in which we parameterize the IGW signal in terms of four parameters: the tensor-to-scalar ratio r r r n t n_t n t f C M B = 0.05 M p c − 1 / ( 2 π a 0 ) ≃ 7.73 × 10 − 17 H z f_{\rm CMB}=0.05\,\mathrm{Mpc}^{-1}/(2\pi a_0)\simeq7.73\times10^{-17}\,\mathrm{Hz} f CMB = 0.05 Mpc − 1 / ( 2 π a 0 ) ≃ 7.73 × 1 0 − 17 Hz T r h T_{\rm rh} T rh N r h N_{\rm rh} N rh a 0 a_0 a 0 f C M B f_{\rm CMB} f CMB a ( t ) a(t) a ( t ) a 0 = 1 a_0=1 a 0 = 1
We do not impose the standard consistency relation between r r r n t n_t n t log 10 r ∈ [ − 40 , 0 ] \log_{10}r\in[-40,0] log 10 r ∈ [ − 40 , 0 ] n t ∈ [ 0 , 6 ] n_t\in[0,6] n t ∈ [ 0 , 6 ] [27 , 84 , 240 , 99 , 67 ] or in other non-minimal inflation models (see [249 , 269 , 191 , 110 , 127 ] for an incomplete list).
Similarly, we allow for more flexibility for T r h T_{\rm rh} T rh N r h N_{\rm rh} N rh N r h N_{\rm rh} N rh [217 ]
N r h = 1 3 ( 1 + w r h ) ln ( 3 H e n d 2 M P l 2 π 2 / 30 g ∗ r h T r h 4 ) , (14) N_{\rm rh}=\frac{1}{3(1+w_{\rm rh})}\ln\left(\frac{3H_{\rm end}^{2}M_{\rm Pl}^{2}}{\pi^{2}/30\,g_{*}^{\rm rh}T_{\rm rh}^{4}}\right), \tag*{(14)} N rh = 3 ( 1 + w rh ) 1 ln ( π 2 /30 g ∗ rh T rh 4 3 H end 2 M Pl 2 ) , (14) where H e n d H_{\rm end} H end w r h w_{\rm rh} w rh M P l ≃ 2.44 × 10 18 G e V M_{\rm Pl}\simeq2.44\times10^{18}\,\mathrm{GeV} M Pl ≃ 2.44 × 1 0 18 GeV g ∗ r h g_{*}^{\rm rh} g ∗ rh w r h = 0 w_{\rm rh}=0 w rh = 0 H e n d H_{\rm end} H end
N r h n a i v e ≃ 1 3 ln ( 3 H n a i v e 2 M P l 2 π 2 / 30 g ∗ r h T r h 4 ) , (15) N_{\rm rh}^{\rm naive}\simeq\frac{1}{3}\ln\left(\frac{3H_{\rm naive}^{2}M_{\rm Pl}^{2}}{\pi^{2}/30\,g_{*}^{\rm rh}T_{\rm rh}^{4}}\right), \tag*{(15)} N rh naive ≃ 3 1 ln ( π 2 /30 g ∗ rh T rh 4 3 H naive 2 M Pl 2 ) , (15) where H n a i v e H_{\rm naive} H naive r r r A s ≃ 2.10 × 10 − 9 A_s\simeq2.10\times10^{-9} A s ≃ 2.10 × 1 0 − 9 [16 ] ,
H e n d ≈ H n a i v e = ( π 2 2 r A s ) 1 / 2 M P l . (16) H_{\rm end}\approx H_{\rm naive}=\left(\frac{\pi^{2}}{2}rA_s\right)^{1/2}M_{\rm Pl}. \tag*{(16)} H end ≈ H naive = ( 2 π 2 r A s ) 1/2 M Pl . (16) However, we already assume nonminimal dynamics in order to motivate a strongly blue-tilted h 2 Ω G W h^2\Omega_{\rm GW} h 2 Ω GW H e n d H_{\rm end} H end N r h N_{\rm rh} N rh r r r T r h T_{\rm rh} T rh (15) and (16) . This flexibility provides us with more parametric freedom that we can use in order to ensure that the IGW signal does not violate constraints on the amplitude of the stochastic GWB set by the LIGO–Virgo–KAGRA (LVK) Collaboration ([5 ] ) and on the amount of dark radiation, i.e., the effective number of neutrino species, N e f f N_{\rm eff} N eff [250 , 313 ] . As for the latter constraint, we specifically work with Δ N e f f = ρ D R / ρ ν \Delta N_{\rm eff}=\rho_{\rm DR}/\rho_\nu Δ N eff = ρ DR / ρ ν ρ D R \rho_{\rm DR} ρ DR ρ ν \rho_\nu ρ ν Δ N e f f \Delta N_{\rm eff} Δ N eff e + e − e^+e^- e + e − N e f f = N e f f S M + Δ N e f f N_{\rm eff}=N_{\rm eff}^{\rm SM}+\Delta N_{\rm eff} N eff = N eff SM + Δ N eff N e f f S M ≃ 3.0440 N_{\rm eff}^{\rm SM}\simeq3.0440 N eff SM ≃ 3.0440 [50 ] .
Under the assumptions outlined above, we are able to model the IGW spectrum at PTA frequencies as
Ω G W i n f ( f ) = Ω r 24 ( g ∗ ( f ) g ∗ 0 ) ( g ∗ , s 0 g ∗ , s ( f ) ) 4 / 3 P t ( f ) T ( f ) . (17) \Omega_{\rm GW}^{\rm inf}(f)=\frac{\Omega_r}{24}\left(\frac{g_*(f)}{g_*^0}\right)\left(\frac{g_{*,s}^0}{g_{*,s}(f)}\right)^{4/3}\mathcal{P}_t(f)\mathcal{T}(f). \tag*{(17)} Ω GW inf ( f ) = 24 Ω r ( g ∗ 0 g ∗ ( f ) ) ( g ∗ , s ( f ) g ∗ , s 0 ) 4/3 P t ( f ) T ( f ) . (17) Here Ω r / g ∗ 0 ≃ 2.72 × 10 − 5 \Omega_r/g_*^0\simeq2.72\times10^{-5} Ω r / g ∗ 0 ≃ 2.72 × 1 0 − 5 g ∗ , s 0 ≃ 3.93 g_{*,s}^0\simeq3.93 g ∗ , s 0 ≃ 3.93 g ∗ ( f ) g_*(f) g ∗ ( f ) g ∗ , s ( f ) g_{*,s}(f) g ∗ , s ( f ) k = 2 π a 0 f k = 2\pi a_0 f k = 2 π a 0 f g ∗ ( f ) g_*(f) g ∗ ( f ) g ∗ , s ( f ) g_{*,s}(f) g ∗ , s ( f ) [268 ] , g ∗ ( T ) g_*(T) g ∗ ( T ) g ∗ , s ( T ) g_{*,s}(T) g ∗ , s ( T ) k = a ( T ) H ( T ) k = a(T)H(T) k = a ( T ) H ( T ) g ∗ g_* g ∗ g ∗ , s g_{*,s} g ∗ , s g ∗ ( f ) g_*(f) g ∗ ( f ) g ∗ , s ( f ) g_{*,s}(f) g ∗ , s ( f ) g ∗ ( T ) g_*(T) g ∗ ( T ) g ∗ , s ( T ) g_{*,s}(T) g ∗ , s ( T )
The inflationary dynamics give rise to the primordial tensor power spectrum P t \mathcal{P}_t P t T \mathcal{T} T n t n_t n t
P t ( f ) = r A s ( f f C M B ) n t . (18) \mathcal{P}_t(f) = r A_s \left(\frac{f}{f_{\mathrm{CMB}}}\right)^{n_t}. \tag*{(18)} P t ( f ) = r A s ( f CMB f ) n t . (18) Meanwhile, the only relevant contribution to T \mathcal{T} T
T ( f ) ≈ Θ ( f e n d − f ) 1 − 0.22 ( f / f r h ) 1.5 + 0.65 ( f / f r h ) 2 . (19) \mathcal{T}(f) \approx\frac{\Theta(f_{\mathrm{end}}-f)}{1-0.22(f/f_{\mathrm{rh}})^{1.5}+0.65(f/f_{\mathrm{rh}})^2}. \tag*{(19)} T ( f ) ≈ 1 − 0.22 ( f / f rh ) 1.5 + 0.65 ( f / f rh ) 2 Θ ( f end − f ) . (19) Here the fit function in the denominator of this expression is taken from Kuroyanagi et al. [205 ] (2021) and describes the spectral turnover, f n t → f n t − 2 f^{n_t} \to f^{n_t-2} f n t → f n t − 2 f ∼ f r h f \sim f_{\mathrm{rh}} f ∼ f rh
f r h = 1 2 π ( g ∗ , s 0 g ∗ , s r h ) 1 / 3 ( π 2 g ∗ r h 90 ) 1 / 2 T r h T 0 M P l , (20) f_{\mathrm{rh}} = \frac{1}{2\pi}\left(\frac{g_{*,s}^0}{g_{*,s}^{\mathrm{rh}}}\right)^{1/3}\left(\frac{\pi^2 g_*^{\mathrm{rh}}}{90}\right)^{1/2}\frac{T_{\mathrm{rh}}T_0}{M_{\mathrm{Pl}}}, \tag*{(20)} f rh = 2 π 1 ( g ∗ , s rh g ∗ , s 0 ) 1/3 ( 90 π 2 g ∗ rh ) 1/2 M Pl T rh T 0 , (20) with g ∗ r h = g ∗ ( T r h ) g_*^{\mathrm{rh}} = g_*(T_{\mathrm{rh}}) g ∗ rh = g ∗ ( T rh ) g ∗ , s r h = g ∗ , s ( T r h ) g_{*,s}^{\mathrm{rh}} = g_{*,s}(T_{\mathrm{rh}}) g ∗ , s rh = g ∗ , s ( T rh ) T 0 ≃ 2.73 K T_0 \simeq2.73\ \mathrm{K} T 0 ≃ 2.73 K [124 ] ). The Heaviside theta function in Eq. (19) denotes the end-point of the IGW spectrum at f = f e n d f = f_{\mathrm{end}} f = f end
f e n d = 1 2 π ( g ∗ , s 0 g ∗ , s r h ) 1 / 3 ( π 2 g ∗ r h 90 ) 1 / 3 T r h 1 / 3 H e n d 1 / 3 T 0 M P l 2 / 3 . (21) f_{\mathrm{end}} = \frac{1}{2\pi}\left(\frac{g_{*,s}^0}{g_{*,s}^{\mathrm{rh}}}\right)^{1/3}\left(\frac{\pi^2 g_*^{\mathrm{rh}}}{90}\right)^{1/3}\frac{T_{\mathrm{rh}}^{1/3}H_{\mathrm{end}}^{1/3}T_0}{M_{\mathrm{Pl}}^{2/3}}. \tag*{(21)} f end = 2 π 1 ( g ∗ , s rh g ∗ , s 0 ) 1/3 ( 90 π 2 g ∗ rh ) 1/3 M Pl 2/3 T rh 1/3 H end 1/3 T 0 . (21) For fixed T r h T_{\mathrm{rh}} T rh f e n d f_{\mathrm{end}} f end H e n d H_{\mathrm{end}} H end N r h N_{\mathrm{rh}} N rh (14) (recall that we set w r h = 0 w_{\mathrm{rh}} = 0 w rh = 0 f e n d f_{\mathrm{end}} f end
Instead, we constrain its maximally allowed value (i.e., N r h N_{\mathrm{rh}} N rh r r r n t n_t n t T r h T_{\mathrm{rh}} T rh N e f f N_{\mathrm{eff}} N eff f r h f_{\mathrm{rh}} f rh (20) , we have f r h ∼ 30 n H z ( T r h / 1 G e V ) f_{\mathrm{rh}} \sim30\ \mathrm{nHz}\ (T_{\mathrm{rh}}/1\ \mathrm{GeV}) f rh ∼ 30 nHz ( T rh /1 GeV ) T r h = 1 G e V T_{\mathrm{rh}} = 1\ \mathrm{GeV} T rh = 1 GeV T B B N ∼ 1 M e V T_{\mathrm{BBN}} \sim1\ \mathrm{MeV} T BBN ∼ 1 MeV log 10 ( T r h / 1 G e V ) ∈ [ − 3 , + 3 ] \log_{10}(T_{\mathrm{rh}}/1\ \mathrm{GeV}) \in[-3,+3] log 10 ( T rh /1 GeV ) ∈ [ − 3 , + 3 ]
Results and Discussion
We now discuss the outcome of our Bayesian fit analysis. First, we note that the IGW model fits the NG15 data slightly better than the baseline SMBHB model. This is evident from the Bayes factor that we find for the IGW versus SMBHB model comparison, B = 8.8 ± 0.3 \mathcal{B} = 8.8 \pm0.3 B = 8.8 ± 0.3 B = 7.9 ± 0.3 \mathcal{B} = 7.9 \pm0.3 B = 7.9 ± 0.3
The reconstructed posterior distributions for the parameters of the IGW model and its IGW+SMBHB extension are shown in Figure 5 .1 n t n_t n t r r r
Figure 5. Reconstructed posterior distributions for the parameters of the IGW (blue) and IGW+SMBHB (red) models. On the diagonal of the corner plot, we report the 1D marginalized distributions for the parameters with the 68% Bayesian credible intervals (vertical lines), while the off-diagonal panels show the 68% (darker) and 95% (lighter) Bayesian credible regions in the 2D posterior distributions. We construct all credible intervals and regions by integrating over the regions of highest posterior density. The gray-shaded regions are strongly disfavored by the NG15 data as they result in a K K K 1 / 10 1/10 1/10 (10) ]. The black-shaded region results in a violation of the N e f f N_{\mathrm{eff}} N eff (23) (see Appendix § 18 .§ 1 and Figure 22 ), assuming N r h = 0 N_{\mathrm{rh}} = 0 N rh = 0 N r h = 5 N_{\mathrm{rh}} = 5 N rh = 5 N r h = 10 N_{\mathrm{rh}} = 10 N rh = 10 Figure 21 in Appendix § 18 .§ 1 shows an extended version of this plot that includes the SMBHB parameters A B H B A_{\mathrm{BHB}} A BHB γ B H B \gamma_{\mathrm{BHB}} γ BHB
Figure 21. Same as Figure 5 but including the SMBHB parameters A B H B A_{\mathrm{BHB}} A BHB γ B H B \gamma_{\mathrm{BHB}} γ BHB A B H B A_{\mathrm{BHB}} A BHB γ B H B \gamma_{\mathrm{BHB}} γ BHB
n t = − 0.14 log 10 r + 0.58 , (22) n_t = -0.14\log_{10} r + 0.58, \tag*{(22)} n t = − 0.14 log 10 r + 0.58 , (22) and which can be explained as follows: the IGW interpretation of the PTA signal requires the primordial tensor power spectrum P t \mathcal{P}_t P t O ( 10 − 4 ) \mathcal{O}(10^{-4}) O ( 1 0 − 4 ) r ( f P T A / f C M B ) n t r(f_{\mathrm{PTA}}/f_{\mathrm{CMB}})^{n_t} r ( f PTA / f CMB ) n t (18) and thus allows us to estimate the coefficients in Eq. (22) as 1 / log 10 ( f P T A / f C M B ) ≈ 0.14 1/\log_{10}(f_{\mathrm{PTA}}/f_{\mathrm{CMB}}) \approx0.14 1/ log 10 ( f PTA / f CMB ) ≈ 0.14 log 10 ( P t / A s ) / log 10 ( f P T A / f C M B ) ≈ 0.58 \log_{10}(\mathcal{P}_{t}/A_{s})/\log_{10}(f_{\mathrm{PTA}}/f_{\mathrm{CMB}}) \approx0.58 log 10 ( P t / A s ) / log 10 ( f PTA / f CMB ) ≈ 0.58 f P T A = 1 n H z f_{\mathrm{PTA}} = 1\,\mathrm{nHz} f PTA = 1 nHz P t = 0.3 × 10 − 4 \mathcal{P}_{t} = 0.3 \times10^{-4} P t = 0.3 × 1 0 − 4
In addition to the strong covariance between n t n_{t} n t r r r ( T r h , n t ) (T_{\mathrm{rh}}, n_{t}) ( T rh , n t ) ( T r h , r ) (T_{\mathrm{rh}}, r) ( T rh , r ) ( log 10 T r h / G e V , n t ) ∼ ( − 0.5 , 2.75 ) (\log_{10} T_{\mathrm{rh}}/\mathrm{GeV}, n_{t}) \sim(-0.5, 2.75) ( log 10 T rh / GeV , n t ) ∼ ( − 0.5 , 2.75 ) ( log 10 T r h / G e V , log 10 r ) ∼ ( − 0.5 , − 15 ) (\log_{10} T_{\mathrm{rh}}/\mathrm{GeV}, \log_{10} r) \sim(-0.5, -15) ( log 10 T rh / GeV , log 10 r ) ∼ ( − 0.5 , − 15 ) Figure 5 indicate that the IGW model can operate in two regimes: for T r h ≫ 1 G e V T_{\mathrm{rh}} \gg1\,\mathrm{GeV} T rh ≫ 1 GeV f r h f_{\mathrm{rh}} f rh T r h ≪ 1 G e V T_{\mathrm{rh}} \ll1\,\mathrm{GeV} T rh ≪ 1 GeV f r h f_{\mathrm{rh}} f rh n t n_{t} n t n t − 2 n_{t}-2 n t − 2 ( T r h , n t ) (T_{\mathrm{rh}}, n_{t}) ( T rh , n t ) ( T r h , r ) (T_{\mathrm{rh}}, r) ( T rh , r ) g ∗ g_{*} g ∗ g ∗ , s g_{*,s} g ∗ , s (19) .
At small or large values of the reheating temperature, the posterior distributions develop approximately flat directions along the T r h T_{\mathrm{rh}} T rh ( n t , log 10 r ) ∼ ( 2 , − 10 ) (n_{t}, \log_{10} r) \sim(2, -10) ( n t , log 10 r ) ∼ ( 2 , − 10 ) T r h T_{\mathrm{rh}} T rh ( n t , log 10 r ) ∼ ( 4 , − 20 ) (n_{t}, \log_{10} r) \sim(4, -20) ( n t , log 10 r ) ∼ ( 4 , − 20 ) T r h T_{\mathrm{rh}} T rh (22) , the reflection symmetry discussed above, and the fact that, past a certain point, raising or lowering T r h T_{\mathrm{rh}} T rh n t = 2 n_{t} = 2 n t = 2 T r h T_{\mathrm{rh}} T rh γ = 3 \gamma= 3 γ = 3 A A A γ \gamma γ Figure 1 . The same conclusion holds at low T r h T_{\mathrm{rh}} T rh γ = 3 \gamma= 3 γ = 3 n t − 2 = 2 n_{t} - 2 = 2 n t − 2 = 2
Finally, we derive the N e f f N_{\mathrm{eff}} N eff N e f f N_{\mathrm{eff}} N eff f e n d f_{\mathrm{end}} f end
Here, f B B N ∼ 10 − 12 H z f_{\mathrm{BBN}} \sim10^{-12}\,\mathrm{Hz} f BBN ∼ 1 0 − 12 Hz T ∼ 10 − 4 G e V T \sim10^{-4}\,\mathrm{GeV} T ∼ 1 0 − 4 GeV [70 ] ), and Δ N e f f max ∼ f e w × 0.1 \Delta N_{\mathrm{eff}}^{\max} \sim\mathrm{few} \times0.1 Δ N eff m a x ∼ few × 0.1 f B B N = 10 − 12 H z f_{\mathrm{BBN}} = 10^{-12}\,\mathrm{Hz} f BBN = 1 0 − 12 Hz Δ N e f f max = 0.5 \Delta N_{\mathrm{eff}}^{\max} = 0.5 Δ N eff m a x = 0.5 Δ N e f f max \Delta N_{\mathrm{eff}}^{\max} Δ N eff m a x [250 , 313 ] ). For given values of the parameters T r h T_{\mathrm{rh}} T rh r r r n t n_{t} n t f e n d max f_{\mathrm{end}}^{\max} f end m a x (23) . If this is the case, f e n d max f_{\mathrm{end}}^{\max} f end m a x e e e N r h max N_{\mathrm{rh}}^{\max} N rh m a x (14) and (21) . A second constraint on N r h N_{\mathrm{rh}} N rh [5 ] ),
∫ f B B N f e n d d f f h 2 Ω G W i n f ( f ) ≲ 5.6 × 10 − 6 Δ N e f f max . (23) \int_{f_{\mathrm{BBN}}}^{f_{\mathrm{end}}} \frac{\mathrm{d}f}{f} h^{2}\Omega_{\mathrm{GW}}^{\mathrm{inf}}(f) \lesssim5.6 \times10^{-6}\Delta N_{\mathrm{eff}}^{\max}. \tag*{(23)} ∫ f BBN f end f d f h 2 Ω GW inf ( f ) ≲ 5.6 × 1 0 − 6 Δ N eff m a x . (23) Ω G W ≤ 1.7 × 10 − 8 at f L V K = 25 H z , (24) \Omega_{\mathrm{GW}} \leq1.7 \times10^{-8} \quad\text{at} \quad f_{\mathrm{LVK}} = 25\,\mathrm{Hz}, \tag*{(24)} Ω GW ≤ 1.7 × 1 0 − 8 at f LVK = 25 Hz , (24) assuming a flat GW spectrum. For a power-law spectrum, Ω G W = Ω G W , α ( f / f l v k ) α \Omega_{\mathrm{GW}} = \Omega_{\mathrm{GW},\alpha} (f/f_{\mathrm{lvk}})^\alpha Ω GW = Ω GW , α ( f / f lvk ) α α = n t − 2 \alpha= n_t - 2 α = n t − 2 f ≫ f r h f \gg f_{\mathrm{rh}} f ≫ f rh [205 ] (2021),
Ω G W ≤ 1.7 × 10 − 8 ( 5 − 2 α 5 ) 1 / 2 ( 20 H z f l v k ) − α . (25) \Omega_{\mathrm{GW}} \leq1.7 \times10^{-8} \left(\frac{5-2\alpha}{5}\right)^{1/2} \left(\frac{20\,\mathrm{Hz}}{f_{\mathrm{lvk}}}\right)^{-\alpha}. \tag*{(25)} Ω GW ≤ 1.7 × 1 0 − 8 ( 5 5 − 2 α ) 1/2 ( f lvk 20 Hz ) − α . (25) which approximately holds if α ≪ 5 / 2 \alpha\ll5/2 α ≪ 5/2 T r h T_{\mathrm{rh}} T rh r r r n t n_t n t Ω G W i n f \Omega_{\mathrm{GW}}^{\mathrm{inf}} Ω GW inf f l v k f_{\mathrm{lvk}} f lvk N r h N_{\mathrm{rh}} N rh f e n d max = 20 H z f_{\mathrm{end}}^{\max} = 20\,\mathrm{Hz} f end m a x = 20 Hz
The outcome of our analysis is shown in Figure 22 in Appendix C.1. In both plots of this figure, the dependence on the tensor-to-scalar ratio r r r (22) . In the left panel of Figure 22 , we present contour lines of N r h max N_{\mathrm{rh}}^{\max} N rh m a x N e f f N_{\mathrm{eff}} N eff T r h − n t T_{\mathrm{rh}}-n_t T rh − n t N r h N_{\mathrm{rh}} N rh N r h max < N r h N_{\mathrm{rh}}^{\max} < N_{\mathrm{rh}} N rh m a x < N rh N e f f N_{\mathrm{eff}} N eff n t n_t n t T r h T_{\mathrm{rh}} T rh N r h < N r h max N_{\mathrm{rh}} < N_{\mathrm{rh}}^{\max} N rh < N rh m a x N r h max N_{\mathrm{rh}}^{\max} N rh m a x N r h max ∼ O ( 10 ⋯ 100 ) N_{\mathrm{rh}}^{\max} \sim\mathcal{O}(10\cdots100) N rh m a x ∼ O ( 10 ⋯ 100 ) N r h ∼ O ( 1 ⋯ 10 ) N_{\mathrm{rh}} \sim\mathcal{O}(1\cdots10) N rh ∼ O ( 1 ⋯ 10 ) Figure 22 , we compare our result to the naive expectation N r h n a i v e N_{\mathrm{rh}}^{\mathrm{naive}} N rh naive (15) ]. Across large parts of parameter space, we find that N r h n a i v e N_{\mathrm{rh}}^{\mathrm{naive}} N rh naive N r h n a i v e ≫ 10 N_{\mathrm{rh}}^{\mathrm{naive}} \gg10 N rh naive ≫ 10
Figure 22. Constraints on the parameters of inflationary GWs, i.e., the IGW model discussed in § 5.1 . The black solid lines refer to the N e f f N_{\rm eff} N eff N r h max N_{\rm rh}^{\max} N rh m a x N r h max N_{\rm rh}^{\max} N rh m a x T r h T_{\rm rh} T rh n t n_{\rm t} n t § 5.1 . See Appendix § 18 .§ 1 for details.
In Figure 5 , we highlight the constraints on T r h T_{\mathrm{rh}} T rh n t n_t n t T r h T_{\mathrm{rh}} T rh r r r N e f f N_{\mathrm{eff}} N eff N r h N_{\mathrm{rh}} N rh N r h = 0 N_{\mathrm{rh}} = 0 N rh = 0 N r h = 0 N_{\mathrm{rh}} = 0 N rh = 0 T r h − n t T_{\mathrm{rh}}-n_t T rh − n t T r h T_{\mathrm{rh}} T rh n t n_t n t N r h = 5 N_{\mathrm{rh}} = 5 N rh = 5 e e e Figure 22 in Appendix C.1.
In view of Figure 5 and Figure 22 , we conclude that the IGW model is indeed capable of fitting the NANOGrav signal across large regions in parameter space. An interesting viable scenario consists, e.g., of a large tensor spectral index, n t ∼ 3 ⋯ 4 n_t \sim3\cdots4 n t ∼ 3 ⋯ 4 r ∼ 10 − ( 23 ⋯ 16 ) r \sim10^{-(23\cdots16)} r ∼ 1 0 − ( 23 ⋯ 16 ) T r h ∼ 10 − ( 3 ⋯ 0 ) G e V T_{\mathrm{rh}} \sim10^{-(3\cdots0)}\,\mathrm{GeV} T rh ∼ 1 0 − ( 3 ⋯ 0 ) GeV e e e N r h ≲ 20 N_{\mathrm{rh}} \lesssim20 N rh ≲ 20
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 A_s \simeq2.10 \times10^{-9} A s ≃ 2.10 × 1 0 − 9 k C M B = 0.05 M p c − 1 k_{\mathrm{CMB}} = 0.05\,\mathrm{Mpc}^{-1} k CMB = 0.05 Mpc − 1 [16 ] ). If we naively extrapolate this value down to smaller scales, assuming a fixed and slightly red-tilted h 2 Ω G W h^2\Omega_{\mathrm{GW}} h 2 Ω GW n s ∼ 0.96 n_s \sim0.96 n s ∼ 0.96 [131 , 118 , 44 ] ). 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) ([73 , 129 , 168 , 170 , 301 , 115 ] ). 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.
In 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) ([225 , 224 , 232 , 26 , 46 ] (Matarrese et al. 1994)) (see Domènech [100 ] 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 correlated with the spectral shape of the SIGW signal, as both phenomena are controlled by the scalar spectrum generated during inflation [318 ] . For earlier works on the PBH/SIGW interpretation of the signal in recent PTA data sets, see Vaskonen & Veermäe (2021); [95 ] ; Kohri & Terada (2021). For earlier Bayesian searches for an SIGW signal in PTA data sets, see [78 ] ; [53 ] ; [323 ] ; Yi & Fei (2023); [91 ] , and for a search in LVK data, see [265 ] .
In our analysis, we consider SIGWs in the PTA band and model the associated GW spectrum as follows:
Ω G W i n d ( f ) = Ω r ( g ∗ ( f ) g ∗ 0 ) ( g ∗ , s 0 g ∗ , s ( f ) ) 4 / 3 Ω ‾ G W i n d ( f ) , (26) \Omega_{\mathrm{GW}}^{\mathrm{ind}}(f)=\Omega_r\left(\frac{g_*(f)}{g_*^0}\right)\left(\frac{g_{*,s}^0}{g_{*,s}(f)}\right)^{4/3}\overline{\Omega}_{\mathrm{GW}}^{\mathrm{ind}}(f), \tag*{(26)} Ω GW ind ( f ) = Ω r ( g ∗ 0 g ∗ ( f ) ) ( g ∗ , s ( f ) g ∗ , s 0 ) 4/3 Ω GW ind ( f ) , (26) K ( 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 u v u 2 + v 2 − 3 ) 2 + π 2 Θ ( u + v − 3 ) ] . (27) \mathcal{K}(u,v)=\frac{3\left(4v^2-\left(1+v^2-u^2\right)^2\right)^2\left(u^2+v^2-3\right)^4}{1024u^8v^8}\left[\left(\ln\left|\frac{3-(u+v)^2}{3-(u-v)^2}\right|-\frac{4uv}{u^2+v^2-3}\right)^2+\pi^2\Theta(u+v-\sqrt{3})\right]. \tag*{(27)} K ( u , v ) = 1024 u 8 v 8 3 ( 4 v 2 − ( 1 + v 2 − u 2 ) 2 ) 2 ( u 2 + v 2 − 3 ) 4 [ ( ln 3 − ( u − v ) 2 3 − ( u + v ) 2 − u 2 + v 2 − 3 4 uv ) 2 + π 2 Θ ( u + v − 3 ) ] . (27) The expression in Eq. (30) illustrates the dependence of the SIGW signal on the scalar input spectrum; in particular, it shows that Ω G W i n d \Omega_{\mathrm{GW}}^{\mathrm{ind}} Ω GW ind Ω G W i n d ∝ P R 2 \Omega_{\mathrm{GW}}^{\mathrm{ind}}\propto\mathcal{P}_{\mathcal{R}}^2 Ω GW ind ∝ P R 2 (30) assumes Gaussian perturbations, whose statistics are fully described by the power spectrum P R \mathcal{P}_{\mathcal{R}} P R [239 ] ; [65 ] ; [287 ] ; [40 ] ; [318 ] ; [39 ] ; [11 ] ; [121 ] .
To remain as model-independent as possible, we refrain from choosing a particular inflation model capable of generating an enhanced spectrum P R \mathcal{P}_{\mathcal{R}} P R P R \mathcal{P}_{\mathcal{R}} P R
SIGW-DELTA: Sharp spectral feature in P R \mathcal{P}_{\mathcal{R}} P R delta function in logarithmic k k k
P R ( k ) = A δ ( ln k − ln k ∗ ) . (28) \mathcal{P}_{\mathcal{R}}(k)=A\delta(\ln k-\ln k_*). \tag*{(28)} P R ( k ) = A δ ( ln k − ln k ∗ ) . (28) SIGW-GAUSS: Broad spectral feature in P R \mathcal{P}_{\mathcal{R}} P R Gaussian peak in logarithmic k k k
P R ( k ) = A 2 π Δ exp [ − 1 2 ( ln k − ln k ∗ Δ ) 2 ] . (29) \mathcal{P}_{\mathcal{R}}(k)=\frac{A}{\sqrt{2\pi}\Delta}\exp\left[-\frac{1}{2}\left(\frac{\ln k-\ln k_*}{\Delta}\right)^2\right]. \tag*{(29)} P R ( k ) = 2 π Δ A exp [ − 2 1 ( Δ ln k − ln k ∗ ) 2 ] . (29) where 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,
Ω ‾ G W i n d ( f ) = ∫ 0 ∞ d v ∫ ∣ 1 − v ∣ 1 + v d u K ( u , v ) P R ( u k ) P R ( v k ) . (30) \overline{\Omega}_{\mathrm{GW}}^{\mathrm{ind}}(f)=\int_0^\infty\mathrm{d}v\int_{\lvert1-v\rvert}^{1+v}\mathrm{d}u\,\mathcal{K}(u,v)\mathcal{P}_{\mathcal{R}}(uk)\mathcal{P}_{\mathcal{R}}(vk). \tag*{(30)} Ω GW ind ( f ) = ∫ 0 ∞ d v ∫ ∣ 1 − v ∣ 1 + v d u K ( u , v ) P R ( u k ) P R ( v k ) . (30) The present-day GW frequency is related to the comoving wavenumber by f = k / ( 2 π a 0 ) f=k/(2\pi a_0) f = k / ( 2 π a 0 ) P R \mathcal{P}_{\mathcal{R}} P R R \mathcal{R} R K \mathcal{K} K [117 , 193 , 137 ] (Pi & Sasaki 2020).
SIGW-BOX: Flat and continuous spectral feature in P R \mathcal{P}_{\mathcal{R}} P R box function in logarithmic k k k
P R ( k ) = A Θ ( ln k max − ln k ) Θ ( ln k − ln k min ) . (31) \mathcal{P}_{\mathcal{R}}(k)=A\Theta(\ln k_{\max}-\ln k)\Theta(\ln k-\ln k_{\min}). \tag*{(31)} P R ( k ) = A Θ ( ln k m a x − ln k ) Θ ( ln k − ln k m i n ) . (31) Note that the Gaussian power spectrum in logarithmic k k k k k k Table 3 in Appendix § 16 , where we use again f = k / ( 2 π a 0 ) f=k/(2\pi a_0) f = k / ( 2 π a 0 ) (28) , Eq. (29) , or Eq. (31) to evaluate the integrals in Eq. (30) and compute the GW spectrum. For SIGW-GAUSS and SIGW-BOX, the integration needs to be carried out numerically; for SIGW-DELTA, we can resort to the exact analytical expression provided in [301 ] ; [318 ] .
Results and Discussion
We now turn to the outcome of our Bayesian fit analysis. Compared to the IGW model discussed in § 5.1 , we obtain even larger Bayes factors, indicating that SIGWs tend to provide an even better fit of the NG15 data than IGWs. Specifically, we obtain B = 44 ± 2 \mathcal{B}=44\pm2 B = 44 ± 2 B = 57 ± 3 \mathcal{B}=57\pm3 B = 57 ± 3 B = 21 ± 1 \mathcal{B}=21\pm1 B = 21 ± 1 B = 38 ± 2 \mathcal{B}=38\pm2 B = 38 ± 2 B = 59 ± 3 \mathcal{B}=59\pm3 B = 59 ± 3 B = 23 ± 1 \mathcal{B}=23\pm1 B = 23 ± 1 Figure 2 ). This improvement in the quality of the fit reflects the fact that the SIGW spectra deviate from a pure power law and thus manage to provide a better fit across the whole frequency range probed by NANOGrav (see Figure 3 ). For SIGW-DELTA, we observe again that adding the SMBHB signal does not improve the quality of the fit. The larger prior volume of SIGW-DELTA+SMBHB compared to SIGW-DELTA therefore results in a slight decrease of the Bayes factor. For the other two SIGW models, the Bayes factors remain more or less the same, within the statistical uncertainty of our bootstrapping analysis.
The reconstructed posterior distributions for all SIGW models under consideration are shown in Figure 6 and Figure 7 . Our first conclusion in view of these results is that a successful explanation of the NANOGrav signal in terms of SIGWs requires the primordial scalar power spectrum to have a large amplitude A A A A A A log 10 A ≳ − 1.55 \log_{10} A \gtrsim-1.55 log 10 A ≳ − 1.55 log 10 A ≳ − 1.43 \log_{10} A \gtrsim-1.43 log 10 A ≳ − 1.43 log 10 A ≳ − 1.90 \log_{10} A \gtrsim-1.90 log 10 A ≳ − 1.90 f min f_{\min} f m i n f max f_{\max} f m a x f ∗ f_* f ∗ Figure 6 and Figure 7 and that are summarized in Table 4 .
Figure 6. Same as in Fig. 5 but for the SIGW-BOX model. The regions above the teal contour lines labeled f P B H = 1 f_{\mathrm{PBH}}=1 f PBH = 1 Figure 23 in Appendix C.2 shows an extended version of this plot that includes the SMBHB parameters A B H B A_{\mathrm{BHB}} A BHB γ B H B \gamma_{\mathrm{BHB}} γ BHB
Figure 23. Same as Figure 6 but including the SMBHB parameters A B H B A_{\mathrm{BHB}} A BHB γ B H B \gamma_{\mathrm{BHB}} γ BHB
Figure 7. Same as in Figure 6 but for the SIGW-DELTA (left panel) and SIGW-GAUSS (right panel) models. The solid and dashed teal contour lines labeled f P B H = 1 f_{\rm PBH}=1 f PBH = 1 Δ = 1 \Delta=1 Δ = 1 Δ = 2 \Delta=2 Δ = 2 Figure 24 in Appendix § 18 .§ 2 shows extended versions of the two plots that include the SMBHB parameters A B H B A_{\rm BHB} A BHB γ B H B \gamma_{\rm BHB} γ BHB
Figure 24. Same as Figure 7 but including the SMBHB parameters A B H B A_{\mathrm{BHB}} A BHB γ B H B \gamma_{\mathrm{BHB}} γ BHB
Table 4. 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.
A notable feature in this context is that the posterior distributions for f min f_{\min} f m i n f max f_{\max} f m a x f ∗ f_* f ∗ Figure 19 ). Beyond the NANOGrav band, the SIGW spectrum keeps increasing until it reaches its maximal value at f ≫ 1 n H z f \gg1\ \mathrm{nHz} f ≫ 1 nHz A A A f min f_{\min} f m i n Figure 6 and the 2D posterior distributions for A A A f ∗ f_* f ∗ Figure 7 . A simultaneous increase in A A A f min f_{\min} f m i n f ∗ f_* f ∗
Figure 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 Figure 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 Figure 3 ).
We also observe that the 2D posterior distributions for A A A f ∗ f_* f ∗
SIGW-GAUSS in the limit Δ → 0 \Delta\to0 Δ → 0 A A A f ∗ f_* f ∗ Figure 7 thus reflects the extra dimension in the parameter space of the SIGW-GAUSS model and the additional parametric freedom that comes with it.
Finally, we comment on the bounds on the parameter space of the three SIGW models. As in the case of the IGW model, we identify regions of the parameter space where the K K K (10) falls below 1 / 10 1/10 1/10 Figure 7 and labeled NG15, adding the SIGW contribution to the GW signal leads to much worse agreement with the data than in the case of no SIGW contribution at all. In fact, parameter points in these regions lead to an SIGW signal that exceeds the observed signal—in other words, the gray shaded regions are ruled out because they result in too strong of a SIGW signal. For SIGW-BOX, we are not able to place a K K K A A A f min f_{\min} f m i n f max f_{\max} f m a x A A A f min f_{\min} f m i n f max f_{\max} f m a x f min f_{\min} f m i n
These bounds are an important result of our analysis that is independent of the origin of the NANOGrav signal. They provide valid constraints on the parameter space for both the SIGW-DELTA and SIGW-GAUSS models, regardless of whether these models contribute to the explanation of the observed GWB. In addition, they are complementary to other existing bounds, in particular, the requirement that the mass density of PBHs produced alongside SIGWs must not exceed the energy density of DM,
f P B H ≤ 1 , (32) f_{\rm PBH} \le1, \tag*{(32)} f PBH ≤ 1 , (32) where f P B H = Ω P B H / Ω D M f_{\rm PBH}=\Omega_{\rm PBH}/\Omega_{\rm DM} f PBH = Ω PBH / Ω DM f P B H f_{\rm PBH} f PBH § 18 .§ 2 ; here we simply report our final results in terms of the teal contour lines labeled f P B H = 1 f_{\rm PBH}=1 f PBH = 1 Figure 6 and Figure 7 . For SIGW-BOX, the PBH bound in the A A A f max f_{\max} f m a x f min = 10 − 11 H z f_{\min}=10^{-11}\,\mathrm{Hz} f m i n = 1 0 − 11 Hz A A A f min f_{\min} f m i n f max = 10 − 5 H z f_{\max}=10^{-5}\,\mathrm{Hz} f m a x = 1 0 − 5 Hz A A A f ∗ f_\ast f ∗ Δ = 1 \Delta=1 Δ = 1 Δ = 2 \Delta=2 Δ = 2
In all three cases, we find that the PBH bound is very restrictive, ruling out most of the parameter space favored by the NG15 data. If taken at face value, the PBH bound therefore renders the SIGW interpretation of the NANOGrav signal less likely. However, we stress that the computation of f P B H f_{\rm PBH} f PBH f P B H f_{\rm PBH} f PBH Figure 6 and Figure 7 need to be treated with caution.
In view of the large conceptional uncertainty in the computation of f P B H f_{\rm PBH} f PBH Figure 6 and Figure 7 illustrate that small regions of parameter space do remain viable. In fact, for SIGW-DELTA and SIGW-GAUSS, it is even possible to realize f P B H = 1 f_{\rm PBH}=1 f PBH = 1 f P B H = 1 f_{\rm PBH}=1 f PBH = 1
In closing, we remark that the above conclusions are endangered by the recent claim of a no-go theorem for PBH formation from single-field inflation (Kristiano & Yokoyama [201 ] (2023)). Kristiano & Yokoyama [201 ] (2023) argue that enhanced scalar perturbations at small scales lead to unacceptably large one-loop corrections to the scalar power spectrum at large scales. In terms of the model parameters discussed in this section, this means that the amplitude A A A log 10 A ≪ − 2 \log_{10} A \ll-2 log 10 A ≪ − 2 A s ≃ 2.10 × 10 − 9 A_s \simeq2.10 \times10^{-9} A s ≃ 2.10 × 1 0 − 9 [262 ] ; Firouzjahi & Riotto [123 ] . However, if it should prove to be valid, the requirement of log 10 A ≪ − 2 \log_{10} A \ll-2 log 10 A ≪ − 2 A A A [201 ] (2023) and still lead to a large SIGW signal.
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 [154 , 134 , 133 ] (Hindmarsh et al. 2015)), can source a primordial GWB (see Witten [305 ] ; Hogan [164 ] for seminal work).3 [34 ] ; Xue et al. [311 ] .
Generally, 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. [200 ] 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
and 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.
We parameterize the GWB produced by both sound waves and bubble collisions in a model-independent way in terms of the following phase transition parameters:
T ⋆ T_{\star} T ⋆ ∼ 34 % \sim34\% ∼ 34% [105 ] ). For weak transitions, this temperature coincides with the temperature at the time of bubble nucleation, T n ∼ T ⋆ T_n \sim T_{\star} T n ∼ T ⋆ T n ≪ T ⋆ T_n \ll T_{\star} T n ≪ T ⋆ α ⋆ ≫ 1 \alpha_{\star} \gg1 α ⋆ ≫ 1 T ⋆ T_{\star} T ⋆ T r h ∼ T ⋆ T_{\rm rh} \sim T_{\star} T rh ∼ T ⋆ [105 ] ).
α ⋆ \alpha_{\star} α ⋆ [68 ] ; Ellis et al. [107 ] ).
H ⋆ R ⋆ = R ⋆ / H ⋆ − 1 H_{\star}R_{\star} = R_{\star}/H_{\star}^{-1} H ⋆ R ⋆ = R ⋆ / H ⋆ − 1 H ⋆ − 1 H_{\star}^{-1} H ⋆ − 1 R ⋆ R_{\star} R ⋆ β \beta β H ⋆ R ⋆ = ( 8 π ) 1 / 3 H ⋆ / β H_{\star}R_{\star} = (8\pi)^{1/3}H_{\star}/\beta H ⋆ R ⋆ = ( 8 π ) 1/3 H ⋆ / β
In addition to the parameters T ⋆ T_{\star} T ⋆ α ⋆ \alpha_{\star} α ⋆ H ⋆ R ⋆ H_{\star}R_{\star} H ⋆ R ⋆ v w v_w v w [61 ] ), which, realistically, are the only ones that could lead to a detectable signal in our current data. In particular, we fix v w = 1 v_w = 1 v w = 1 v w = 1 v_w = 1 v w = 1 v w < 1 v_w < 1 v w < 1 v w = 1 v_w = 1 v w = 1 v w = 1 v_w = 1 v w = 1 [68 ] who work with v w = 0.95 v_w = 0.95 v w = 0.95 H ∗ R ∗ = ( 8 π ) 1 / 3 v w H ∗ / β H_*R_*=(8\pi)^{1/3}v_w H_*/\beta H ∗ R ∗ = ( 8 π ) 1/3 v w H ∗ / β v w v_w v w κ s \kappa_s κ s α ∗ \alpha_* α ∗ α ∗ ∼ 0.3 ⋯ 10 \alpha_* \sim0.3\cdots10 α ∗ ∼ 0.3 ⋯ 10 Figure 8 ), this dependence is rather weak (see Figure 13 in Espinosa et al. [116 ] ), which justifies again to keeping v w v_w v w
Figure 8. Same as in Figure 5 but for the PT-BUBBLE (left panel) and PT-SOUND model (right panel). Figure 25 in Appendix § 18 .§ 3 shows the same plots but with the parameter a a a a = 3 a=3 a = 3 Figure 26 and Figure 27 in Appendix § 18 .§ 3 show extended versions of the two plots that include the spectral shape parameters a a a b b b c c c A B H B A_{\mathrm{BHB}} A BHB γ B H B \gamma_{\mathrm{BHB}} γ BHB
Figure 13. Constraints on the local ULDM density for the correlated (solid line) and uncorrelated (dashed line) signals at the 95% credible level (see the discussion in the text). The gray dashed line indicates the predicted DM abundance.
Ω b ( f ) = D Ω ~ b ( α ∗ 1 + α ∗ ) 2 ( H ∗ R ∗ ) 2 S ( f / f b ) (33) \Omega_b(f)=\mathcal{D}\widetilde{\Omega}_b\left(\frac{\alpha_*}{1+\alpha_*}\right)^2(H_*R_*)^2S(f/f_b) \tag*{(33)} Ω b ( f ) = D Ω b ( 1 + α ∗ α ∗ ) 2 ( H ∗ R ∗ ) 2 S ( f / f b ) (33) Ω s ( f ) = D Ω ~ s Υ ( τ s w ) ( κ s α ∗ 1 + α ∗ ) 2 ( H ∗ R ∗ ) S ( f / f s ) . (34) \Omega_s(f)=\mathcal{D}\widetilde{\Omega}_s\Upsilon(\tau_{\rm sw})\left(\frac{\kappa_s\alpha_*}{1+\alpha_*}\right)^2(H_*R_*)S(f/f_s). \tag*{(34)} Ω s ( f ) = D Ω s Υ ( τ sw ) ( 1 + α ∗ κ s α ∗ ) 2 ( H ∗ R ∗ ) S ( f / f s ) . (34) Here Ω ~ b = 0.0049 \widetilde{\Omega}_b=0.0049 Ω b = 0.0049 [176 ] ) and Ω ~ s = 0.036 \widetilde{\Omega}_s=0.036 Ω s = 0.036 κ s = α ∗ / ( 0.73 + 0.083 α ∗ + α ∗ ) \kappa_s=\alpha_* /(0.73+0.083\sqrt{\alpha_*}+\alpha_*) κ s = α ∗ / ( 0.73 + 0.083 α ∗ + α ∗ ) [116 ] ) gives the fraction of the released energy that is transferred to plasma motion in the form of sound waves, and D \mathcal{D} D
D = π 2 90 T 0 4 M P l 2 H 0 2 g ∗ ( g ∗ , s e q g ∗ , s ) 4 / 3 ≃ 1.67 × 10 − 5 . (35) \mathcal{D}=\frac{\pi^2}{90}\frac{T_0^4}{M_{\mathrm{Pl}}^2H_0^2}g_*\left(\frac{g_{*,s}^{\mathrm{eq}}}{g_{*,s}}\right)^{4/3}\simeq1.67\times10^{-5}. \tag*{(35)} D = 90 π 2 M Pl 2 H 0 2 T 0 4 g ∗ ( g ∗ , s g ∗ , s eq ) 4/3 ≃ 1.67 × 1 0 − 5 . (35) We recall that T 0 T_0 T 0 H 0 H_0 H 0 g ∗ g_* g ∗ g ∗ , s g_{*,s} g ∗ , s (35) are evaluated at T = T ∗ T=T_* T = T ∗ g ∗ , s e q g_{*,s}^{\mathrm{eq}} g ∗ , s eq τ s w \tau_{\rm sw} τ sw [105 ] (2020); Guo et al. [145 ] ). In Eq. (34) , this effect is taken into account by the suppression factor
Υ ( τ s w ) = 1 − ( 1 + 2 τ s w H ∗ ) − 1 / 2 , (36) \Upsilon(\tau_{\rm sw})=1-(1+2\tau_{\rm sw}H_*)^{-1/2}, \tag*{(36)} Υ ( τ sw ) = 1 − ( 1 + 2 τ sw H ∗ ) − 1/2 , (36) where the shock formation time scale, τ s w \tau_{\rm sw} τ sw τ s w ≈ R ∗ / U ‾ f \tau_{\rm sw}\approx R_*/\overline{U}_f τ sw ≈ R ∗ / U f U ‾ f 2 ≈ 3 κ s α ∗ / [ 4 ( 1 + α ∗ ) ] \overline{U}_f^2\approx3\kappa_s\alpha_* /[4(1+\alpha_*)] U f 2 ≈ 3 κ s α ∗ / [ 4 ( 1 + α ∗ )] [302 ] ). The functions S b , s S_{b,s} S b , s
f b , s ≃ 48.5 n H z g ∗ 1 / 2 ( g ∗ , s e q g ∗ , s ) 1 / 3 ( T ∗ 1 G e V ) f b , s ∗ R ∗ H ∗ R ∗ , (37) f_{b,s}\simeq48.5\,\mathrm{nHz}\,g_*^{1/2}\left(\frac{g_{*,s}^{\mathrm{eq}}}{g_{*,s}}\right)^{1/3}\left(\frac{T_*}{1\,\mathrm{GeV}}\right)\frac{f_{b,s}^*R_*}{H_*R_*}, \tag*{(37)} f b , s ≃ 48.5 nHz g ∗ 1/2 ( g ∗ , s g ∗ , s eq ) 1/3 ( 1 GeV T ∗ ) H ∗ R ∗ f b , s ∗ R ∗ , (37) where the values of the peak frequencies at the time of GW emission are given by f b ∗ = 0.58 / R ∗ f_b^*=0.58/R_* f b ∗ = 0.58/ R ∗ [176 ] ) and f s ∗ = 1.58 / R ∗ f_s^*=1.58/R_* f s ∗ = 1.58/ R ∗
However, at the level of our Bayesian fit analysis, changes in these prefactors can be absorbed in the temperature scale T ∗ T_* T ∗ (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 α ∗ \alpha_* α ∗
S ( x ) = 1 N ( a + b ) c ( b x − a / c + a x b x / c ) c . (38) S(x)=\frac{1}{\mathcal{N}}\frac{(a+b)^c}{(bx^{-a/c}+axb^{x/c})^c}. \tag*{(38)} S ( x ) = N 1 ( b x − a / c + a x b x / c ) c ( a + b ) c . (38) Here a a a b b b c c c N \mathcal{N} N S S S
N = ( b a ) a / n ( n c b ) c Γ ( a / n ) Γ ( b / n ) n Γ ( c ) , (39) \mathcal{N}=\left(\frac{b}{a}\right)^{a/n}\left(\frac{nc}{b}\right)^c\frac{\Gamma(a/n)\Gamma(b/n)}{n\Gamma(c)}, \tag*{(39)} N = ( a b ) a / n ( b n c ) c n Γ ( c ) Γ ( a / n ) Γ ( b / n ) , (39) where n = ( a + b ) / c n=(a+b)/c n = ( a + b ) / c Γ \Gamma Γ a a a b b b c c c 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. [158 ] (2017); ); Lewicki & Vaskonen [212 ] (2021); ); Cutting et al. [88 , 87 ] ).4
Results and Discussion
The reconstructed posterior distributions for the parameters α ∗ \alpha_* α ∗ T ∗ T_* T ∗ H ∗ R ∗ H_*R_* H ∗ R ∗ Figure 8 , both for the case where the phase transition is assumed to be the only source of GWs (blue contours) and for the scenario where we consider the superposition of the phase transition and SMBHB signals (red contours).5
plots including the posterior distributions for the spectral shape parameters a a a b b b c c c A B H B A_{\mathrm{BHB}} A BHB γ B H B \gamma_{\mathrm{BHB}} γ BHB Figure 26 and Figure 27 in Appendix § 18 .§ 3 .
Figure 26. Same as Figure 8 but including the spectral shape parameters a a a b b b c c c A B H B A_{\mathrm{BHB}} A BHB γ B H B \gamma_{\mathrm{BHB}} γ BHB
Figure 27. Same as Figure 26 but for the PT-SOUND model.
In all analyses, the data prefer a relatively strong and slow phase transition. Specifically, for PT-BUBBLE, we find α ∗ > 1.1 \alpha_* > 1.1 α ∗ > 1.1 H ∗ R ∗ > 0.28 H_*R_* > 0.28 H ∗ R ∗ > 0.28 α ∗ > 1.0 \alpha_* > 1.0 α ∗ > 1.0 H ∗ R ∗ > 0.26 H_*R_* > 0.26 H ∗ R ∗ > 0.26 α ∗ > 0.42 \alpha_* > 0.42 α ∗ > 0.42 H ∗ R ∗ ∈ [ 0.053 , 0.27 ] H_*R_* \in[0.053, 0.27] H ∗ R ∗ ∈ [ 0.053 , 0.27 ] α ∗ ∈ [ 0.46 , 5.4 ] \alpha_* \in[0.46, 5.4] α ∗ ∈ [ 0.46 , 5.4 ] > 0.16 > 0.16 > 0.16 H ∗ R ∗ ∈ [ 0.054 , 0.35 ] H_*R_* \in[0.054, 0.35] H ∗ R ∗ ∈ [ 0.054 , 0.35 ] > 0.0015 > 0.0015 > 0.0015
As can be seen in Figure 3 , for both phase transition models, the reconstructed GWB spectrum tends to peak around the higher frequency bins and fit the signal in the lower frequency bins with the left tail of the spectrum. Specifically, for the PT-BUBBLE model we find T ∗ ∈ [ 0.047 , 0.41 ] T_* \in[0.047, 0.41] T ∗ ∈ [ 0.047 , 0.41 ] T ∗ ∈ [ 4.7 , 33 ] T_* \in[4.7, 33] T ∗ ∈ [ 4.7 , 33 ] T ∗ T_* T ∗ f s ∗ f_s^* f s ∗ f b ∗ f_b^* f b ∗ (37) ). As explained below Eq. (37) , any change in these numerical factors can be reabsorbed in a redefinition of the fit parameter T ∗ T_* T ∗
The inclusion of the SMBHB signal, by adding power to the lowest frequency bins, allows the T ∗ T_* T ∗ T ∗ ∈ [ 4.9 , 50 ] T_* \in[4.9, 50] T ∗ ∈ [ 4.9 , 50 ] 2 × 10 6 2 \times10^6 2 × 1 0 6 T ∗ T_* T ∗ Figure 8 . The drastic increase in the 95% upper limit, on the other hand, is related to the fact that adding the SMBHB signal to the GWB results in a flat plateau region in the posterior distribution of the PT-SOUND model parameters where the NANOGrav signal is mostly explained by the SMBHB contribution to the GWB. The 95% credible regions for the PT-SOUND+SMBHB model cover much of this plateau, which explains their large extent and noisy appearance in Figure 8 . For the PT-BUBBLE model, the inclusion of the SMBHB signal is less significant, and we find T ∗ ∈ [ 0.046 , 0.46 ] T_* \in[0.046, 0.46] T ∗ ∈ [ 0.046 , 0.46 ]
The larger phase transition temperatures observed for the PT-BUBBLE model are a consequence of the smaller value of the peak frequency at the time of emission, f b ∗ f_b^* f b ∗ § 18 .§ 3 , we report the results of the analysis in which the low-frequency slope is set to the value predicted by causality (a = 3 a = 3 a = 3
The corner plots in Figure 8 also illustrate that, as expected from the expression for the peak frequency in Eq. (37) , there is an approximately linear correlation between log 10 T ∗ \log_{10} T_\ast log 10 T ∗ log 10 H ∗ R ∗ \log_{10} H_\ast R_\ast log 10 H ∗ R ∗ α ∗ ≲ 1 \alpha_\ast\lesssim1 α ∗ ≲ 1 α ∗ ∼ 1 / ( H ∗ R ∗ ) \alpha_\ast\sim1/(H_\ast R_\ast) α ∗ ∼ 1/ ( H ∗ R ∗ ) α ∗ 2 ∼ 1 / ( H ∗ R ∗ ) \alpha_\ast^2 \sim1/(H_\ast R_\ast) α ∗ 2 ∼ 1/ ( H ∗ R ∗ ) (33) and Eq. (34) .
We also notice that, for both models, the posterior distribution for T ∗ T_\ast T ∗ [34 ] ). This shift results from the reconstructed h 2 Ω G W h^2\Omega_{\rm GW} h 2 Ω GW [43 , 96 ] ) for both the PT-BUBBLE and PT-SOUND models, which constrain the phase transition parameter space at temperatures around T ∗ ∼ 1 T_\ast\sim1 T ∗ ∼ 1
Instead, we conclude that the reconstructed posterior distribution of T ∗ T_\ast T ∗ [241 ] ; Li et al. (2021)) and (ii) strong first-order phase transitions in a dark sector composed of new BSM degrees of freedom (see, e.g., Nakai et al. [238 ] ; Ratzinger & Schwaller [260 ] ). In view of the NG15 data, both of these options for the particle physics origin of the phase transition signal remain viable. A third option may consist in a strongly supercooled first-order electroweak phase transition (Kobakhidze et al. [189 ] ).
Finally, we report that that, like the models studied in § 5.1 and § 6.1 , the phase transition models provide a better fit of the NG15 data than the base SMBHB model. The Bayes factors for PT-BUBBLE and PT-SOUND versus SMBHB are B = 18.1 ± 0.6 \mathcal{B} = 18.1 \pm0.6 B = 18.1 ± 0.6 B = 3.7 ± 0.1 \mathcal{B} = 3.7 \pm0.1 B = 3.7 ± 0.1 B = 12.6 ± 0.5 \mathcal{B} = 12.6 \pm0.5 B = 12.6 ± 0.5 B = 6.5 ± 0.3 \mathcal{B} = 6.5 \pm0.3 B = 6.5 ± 0.3
tively. An interesting observation in view of these results is that adding the SMBHB contribution to the GWB signal does not help to improve the quality of the fit for PT-BUBBLE—in this case, we find again a decrease in the Bayes factor going from PT-BUBBLE to PT-BUBBLE+SMBHB because of the larger prior volume—but it does lead to a better fit for PT-SOUND. This model benefits from the additional SMBHB contribution because it can add power to the low frequency bins in the GW spectrum that the PT-SOUND model alone struggles to fit well on its own (see Figure 3 ). The reason for this, in turn, is the prior range for the spectral index at low frequencies, a a a a = 1 a = 1 a = 1 a = 3 a = 3 a = 3 Table 3 ). Another consequence of this interplay between the phase transition and SMBHB signals is that the NANOGrav signal may in fact be dominated by SMBHBs. This possibility is realized when the PT-SOUND model parameters fall into the plateau region in Figure 8 (i.e., the red 95% credible regions in the right panel) and the SMBHB parameters are close to log 10 A B H B ∼ − ( 15 ÷ 14 ) \log_{10} A_{\rm BHB} \sim-(15\div14) log 10 A BHB ∼ − ( 15 ÷ 14 ) γ B H B ∼ 3 ÷ 4 \gamma_{\rm BHB} \sim3\div4 γ BHB ∼ 3 ÷ 4 Figure 27 ).
Cosmic strings
MODEL DESCRIPTION
Cosmic strings are effectively 1D topological defects that can form in the early universe as a consequence of a cosmological phase transition (Kibble [188 ] ). Rigorously speaking, the criterion for the formation of a cosmic-string network is that the underlying phase transition must entail the spontaneous breaking of a local or global symmetry and end on a vacuum manifold with a nontrivial first homotopy group. In practice, the most relevant scenario satisfying this criterion is the cosmological breaking of a U ( 1 ) U(1) U ( 1 ) U ( 1 ) U(1) U ( 1 ) U ( 1 ) U(1) U ( 1 ) [297 , 165 , 290 , 7 Bennett &] Bouchet 1991; Caldwell & Allen 1992). For earlier work on the cosmic strings interpretation of the signal in recent PTA data sets, see, e.g., Blasi et al. [60 ] , Ellis & Lewicki [104 ] , and Blanco-Pillado et al. [59 ] .
In our analysis, we study ordinary stable cosmic strings as well as metastable strings, which are unstable against the nucleation of GUT monopoles. In passing, we also comment on cosmic superstrings, which do not have a particle-physics origin but are present in some string-theoretic models of the early universe. Our starting point for describing all these scenarios is the Nambu–Goto action, which treats cosmic strings as featureless 1D objects that can be characterized by a single parameter: their tension, i.e., energy per unit length, μ \mu μ κ \kappa κ μ \mu μ P P P [54 ] ), and we refer the reader to Section 3.4 of Auclair et al. [41 ] for an extended discussion. Furthermore, we only consider the GWB produced by cosmic-string loops in our analysis and disregard the subdominant contribution from long (infinite) strings (see a discussion of this point in Section 4.4 of Auclair et al. [41 ] ).
Shortly after their formation, cosmic strings enter a scaling regime where the total energy stored in the network, ρ c s \rho_{\rm cs} ρ cs Ω c s = ρ c s / ρ c ≈ c o n s t \Omega_{\rm cs} = \rho_{\rm cs}/\rho_{\rm c} \approx\mathrm{const} Ω cs = ρ cs / ρ c ≈ const [157 ] ). This behavior is possible because long strings, stretching over superhorizon distances, frequently intercommute with each other, thereby producing an abundance of closed string loops on subhorizon scales that radiate energy in the form of GWs. These GWs are sourced by the oscillations of the loops under their own tension, as well as by localized features (“cusps” and “kinks”) propagating along the loops. The superposition of the GWs emitted by the individual loops in the network thus results in a stochastic GWB,
Ω G W c s ( f ) = 8 π 3 H 0 2 ( G μ ) 2 ∑ k = 1 k max P k I k ( f ) . (40) \Omega_{\rm GW}^{\rm cs}(f) = \frac{8\pi}{3H_0^2}(G\mu)^2 \sum_{k=1}^{k_{\max}} P_k \mathcal{I}_k(f). \tag*{(40)} Ω GW cs ( f ) = 3 H 0 2 8 π ( G μ ) 2 k = 1 ∑ k m a x P k I k ( f ) . (40) Here the dimensionless factor G μ G\mu G μ log 10 G μ ∈ [ − 14 , − 6 ] \log_{10} G\mu\in[-14,-6] log 10 G μ ∈ [ − 14 , − 6 ] log 10 G μ ∈ [ − 14 , − 1.5 ] \log_{10} G\mu\in[-14,-1.5] log 10 G μ ∈ [ − 14 , − 1.5 ] (40) runs over the harmonic excitations of the closed string loops that, given a loop of length ℓ \ell ℓ f = 2 k / ℓ f = 2k/\ell f = 2 k / ℓ k = 1 k = 1 k = 1 k max = 10 5 k_{\max} = 10^5 k m a x = 1 0 5
The dimensionless factor P k P_k P k G μ 2 G\mu^2 G μ 2 k k k k max k_{\max} k m a x
P k = Γ ζ ( q ) k q , (41) P_k = \frac{\Gamma}{\zeta(q) k^q}, \tag*{(41)} P k = ζ ( q ) k q Γ , (41) where the prefactor Γ / ζ ( q ) \Gamma/\zeta(q) Γ/ ζ ( q ) ∑ k P k = Γ \sum_k P_k = \Gamma ∑ k P k = Γ q q q Γ = 50 \Gamma= 50 Γ = 50 [55 ] ), and discuss different possibilities for q q q
STABLE-C: Stable strings emitting GWs only in the form of GW bursts from cusps on closed loops, q = 4 / 3 q = 4/3 q = 4/3 [290 ] ).
STABLE-K: Stable strings emitting GWs only in the form of GW bursts from kinks on closed loops, q = 5 / 3 q = 5/3 q = 5/3 [90 ] ).
STABLE-M: Stable strings emitting monochromatic GWs at the fundamental frequency f = 2 / ℓ f = 2/\ell f = 2/ ℓ
STABLE-N: Stable strings described by numerical simulations including GWs from cusps and kinks (Blanco-Pillado et al. [56 ] (2015)).
For STABLE-N, P k P_k P k k k k q ≈ 4 / 3 q \approx4/3 q ≈ 4/3 k ≳ 100 k \gtrsim100 k ≳ 100 k k k q q q q ∼ 1.5 q \sim1.5 q ∼ 1.5 k ∼ 10 k \sim10 k ∼ 10 (41) is irrelevant for STABLE-M; for this model, we simply set P k = Γ P_k = \Gamma P k = Γ k = 1 k = 1 k = 1 P k = 0 P_k = 0 P k = 0 [59 ] .
Finally, the frequency-dependent factor I k \mathcal{I}_k I k (40) represents an integral of the number density of closed string loops, n l ( ℓ , t ) n_l(\ell,t) n l ( ℓ , t )
I k ( f ) = 2 k f ∫ t i n i t 0 d t ( a ( t ) a ( t 0 ) ) 5 n l ( 2 k f a ( t ) a ( t 0 ) , t ) . (42) \mathcal{I}_k(f) = \frac{2k}{f} \int_{t_{\rm ini}}^{t_0} \mathrm{d}t \left(\frac{a(t)}{a(t_0)}\right)^5 n_l\left(\frac{2k}{f}\frac{a(t)}{a(t_0)},t\right). \tag*{(42)} I k ( f ) = f 2 k ∫ t ini t 0 d t ( a ( t 0 ) a ( t ) ) 5 n l ( f 2 k a ( t 0 ) a ( t ) , t ) . (42) Here t 0 t_0 t 0 t i n i t_{\rm ini} t ini t i n i t_{\mathrm{ini}} t ini n l n_l n l
n l ( ℓ , t ) = F C ∗ Θ ( t − t ∗ ) Θ ( t ∗ − t i n i ) α ∗ ( α ∗ + Γ G μ + α ˙ ∗ t ∗ ) t ∗ 4 ( a ∗ a ( t ) ) 3 (43) n_l(\ell,t)=\mathcal{F}\frac{C_* \Theta(t-t_*)\Theta(t_*-t_{\mathrm{ini}})}{\alpha_*\left(\alpha_*+\Gamma G\mu+\dot{\alpha}_*t_*\right)t_*^4}\left(\frac{a_*}{a(t)}\right)^3 \tag*{(43)} n l ( ℓ , t ) = F α ∗ ( α ∗ + Γ G μ + α ˙ ∗ t ∗ ) t ∗ 4 C ∗ Θ ( t − t ∗ ) Θ ( t ∗ − t ini ) ( a ( t ) a ∗ ) 3 (43) where F = 0.1 \mathcal{F}=0.1 F = 0.1 [41 ] (Blanco-Pillado et al. 2014)) and an asterisk indicates that the corresponding quantity needs to be evaluated at the time of loop formation, which follows from solving the following relation for t ∗ t_* t ∗
t ∗ = ℓ + Γ G μ t α ∗ + Γ G μ , α ∗ = α ( t ∗ ) . (44) t_*=\frac{\ell+\Gamma G\mu t}{\alpha_*+\Gamma G\mu},\qquad\alpha_*=\alpha(t_*) . \tag*{(44)} t ∗ = α ∗ + Γ G μ ℓ + Γ G μ t , α ∗ = α ( t ∗ ) . (44) The time-dependent functions C C C α \alpha α
C ( t ) = c ~ 2 v ˉ ( t ) ξ 3 ( t ) , α ( t ) = α L ξ ( t ) . (45) C(t)=\frac{\tilde{c}}{\sqrt{2}}\frac{\bar{v}(t)}{\xi^3(t)},\qquad\alpha(t)=\alpha_L\xi(t) . \tag*{(45)} C ( t ) = 2 c ~ ξ 3 ( t ) v ˉ ( t ) , α ( t ) = α L ξ ( t ) . (45) Here c ~ \tilde{c} c ~ α L \alpha_L α L L = ξ t L=\xi t L = ξ t c ~ ≃ 0.23 \tilde{c}\simeq0.23 c ~ ≃ 0.23 α L ≃ 0.37 \alpha_L\simeq0.37 α L ≃ 0.37 [222 , 56 ] (Blanco-Pillado et al. 2014)). With c ~ \tilde{c} c ~ ξ \xi ξ v ˉ \bar{v} v ˉ C C C α \alpha α a a a Λ \Lambda Λ ξ r ≃ 0.27 \xi_r\simeq0.27 ξ r ≃ 0.27 v ˉ r ≃ 0.66 \bar{v}_r\simeq0.66 v ˉ r ≃ 0.66 α ≃ 0.10 \alpha\simeq0.10 α ≃ 0.10 n l n_l n l g ∗ g_* g ∗ g ∗ , s g_{*,s} g ∗ , s
In the case of metastable strings, the loop number density in Eq. (43) receives two correction factors,
n l m e t a ( ℓ , t ) = Θ ( t s − t ∗ ) E ( ℓ , t ) n l ( ℓ , t ) . (46) n_l^{\mathrm{meta}}(\ell,t)=\Theta(t_s-t_*)E(\ell,t)n_l(\ell,t) . \tag*{(46)} n l meta ( ℓ , t ) = Θ ( t s − t ∗ ) E ( ℓ , t ) n l ( ℓ , t ) . (46) The Heaviside function accounts for the fact that we expect loop formation to cease when monopole nucleation becomes efficient and the network transitions from the scaling regime to the decay regime at times around
t s = 1 Γ d 1 / 2 , (47) t_s=\frac{1}{\Gamma_d^{1/2}} , \tag*{(47)} t s = Γ d 1/2 1 , (47) with the decay rate, Γ d \Gamma_d Γ d
Γ d = μ 2 π e − π κ . (48) \Gamma_d=\frac{\mu}{2\pi}e^{-\pi\kappa}. \tag*{(48)} Γ d = 2 π μ e − π κ . (48) Here κ \kappa κ U ( 1 ) U(1) U ( 1 ) κ \sqrt{\kappa} κ U ( 1 ) U(1) U ( 1 )
κ = m G U T μ 1 / 2 ∼ Λ G U T Λ U ( 1 ) . (49) \sqrt{\kappa}=\frac{m_{\mathrm{GUT}}}{\mu^{1/2}}\sim\frac{\Lambda_{\mathrm{GUT}}}{\Lambda_{U(1)}} . \tag*{(49)} κ = μ 1/2 m GUT ∼ Λ U ( 1 ) Λ GUT . (49) Meanwhile, the second new factor in Eq. (46) represents an exponential suppression term, reflecting the depletion of the loop number density in the decaying network,
E ( ℓ , t ) = e − Γ d [ ℓ ( t − t ∗ ) 1 / 2 Γ G μ ( t − t ∗ ) 2 ] . (50) E(\ell,t)=e^{-\Gamma_d[\ell(t-t_*)^{1/2}\Gamma G\mu(t-t_*)^2]} . \tag*{(50)} E ( ℓ , t ) = e − Γ d [ ℓ ( t − t ∗ ) 1/2 Γ G μ ( t − t ∗ ) 2 ] . (50) As evident from Eqs. (48) and (50) , the loop number density of a metastable network is exponentially sensitive to the decay parameter κ \kappa κ κ ≫ 10 \sqrt{\kappa}\gg10 κ ≫ 10 κ ∼ 1 \sqrt{\kappa}\sim1 κ ∼ 1 κ \sqrt{\kappa} κ κ ∈ [ 7.0 , 9.5 ] \sqrt{\kappa}\in[7.0,9.5] κ ∈ [ 7.0 , 9.5 ]
Monopole nucleation in a metastable network results in string segments, dumbbell-like shaped pieces of string with monopoles attached on either end. In many GUT models, these monopoles still carry unconfined flux, such that we actually need to distinguish between monopoles and antimonopoles. Monopoles and antimonopoles with unconfined flux occur, e.g., in GUT models involving, schematically, the symmetry breaking pattern
S U ( 2 ) 1 × U ( 1 ) 2 → U ( 1 ) 1 × U ( 1 ) 2 → U ( 1 ) . (51) SU(2)_1\times U(1)_2\to U(1)_1\times U(1)_2\to U(1). \tag*{(51)} S U ( 2 ) 1 × U ( 1 ) 2 → U ( 1 ) 1 × U ( 1 ) 2 → U ( 1 ) . (51) In this case, we expect the string segments to dissipate most of their energy via the emission of massless vector bosons soon after their formation. Alternatively, monopoles may carry no unconfined flux, in which case they are able to lose energy only via the emission of GWs. In our analysis, we cover both possibilities:
META-L: Metastable strings; monopoles carry unconfined flux; GW emission from loops only.
META-LS: Metastable strings; monopoles carry no unconfined flux; GW emission from loops and segments .
For META-LS, we also require the number density of segments that result from the decaying network; the relevant expressions are summarized in Appendix § 18 .§ 4 .
Moreover, we need to specify the power spectrum of GWs emitted by string segments. Following Buchmuller et al. [62 ] , we use again Eq. (41) and set q = 1 q=1 q = 1 [221 ] ). At the same time, we assume a pure cusp spectrum (q = 4 / 3 q=4/3 q = 4/3
Finally, we also consider cosmic superstrings:
SUPER: Cosmic superstrings; suppressed intercommutation probability P P P q = 4 / 3 q=4/3 q = 4/3
Cosmic superstrings are characterized by a reduced intercommutation probability P P P 10 3 10^{3} 1 0 3 (40) ,
Ω G W c s ( f ) ⟶ 1 P Ω G W c s ( f ) , (52) \Omega_{\mathrm{GW}}^{\mathrm{cs}}(f) \quad\longrightarrow\quad\frac{1}{P}\Omega_{\mathrm{GW}}^{\mathrm{cs}}(f), \tag*{(52)} Ω GW cs ( f ) ⟶ P 1 Ω GW cs ( f ) , (52) where we allow for P P P log 10 P ∈ [ − 3 , 0 ] \log_{10} P \in[-3,0] log 10 P ∈ [ − 3 , 0 ] F F F D D D [172 ] ).
Results and Discussion
We now turn to the outcome of our Bayesian fit analysis, discussing first our results for stable strings. As evident from Figure 2 and Figure 3 , we find that stable strings are not able to provide a better fit of the PTA data than the benchmark SMBHB model. In fact, among all exotic GWB sources considered in this work, stable cosmic strings represent the only case that consistently yields a Bayes factor in favor of the SMBHB interpretation. Comparing the STABLE-C, STABLE-K, STABLE-M, and STABLE-N models to the SMBHB model, we obtain B = 0.277 ± 0.006 \mathcal{B}=0.277\pm0.006 B = 0.277 ± 0.006 B = 0.364 ± 0.008 \mathcal{B}=0.364\pm0.008 B = 0.364 ± 0.008 B = 0.379 ± 0.008 \mathcal{B}=0.379\pm0.008 B = 0.379 ± 0.008 B = 0.307 ± 0.006 \mathcal{B}=0.307\pm0.006 B = 0.307 ± 0.006 B = 0.76 ± 0.01 \mathcal{B}=0.76\pm0.01 B = 0.76 ± 0.01 B = 0.89 ± 0.02 \mathcal{B}=0.89\pm0.02 B = 0.89 ± 0.02 B = 0.84 ± 0.02 \mathcal{B}=0.84\pm0.02 B = 0.84 ± 0.02 B = 0.83 ± 0.01 \mathcal{B}=0.83\pm0.01 B = 0.83 ± 0.01
The reason for the poor performance of the stable-string models is straightforward. In order to explain the relatively large amplitude of the signal, a comparatively large value of the cosmic-string tension G μ G\mu G μ G μ G\mu G μ Figure 20 (see also, e.g., Figure 1 in Blanco-Pillado et al. [56 ] or Figure 1 in Blasi et al. [60 ] ), which clashes with the fact that the data seem to prefer a blue-tilted h 2 Ω G W h^{2}\Omega_{\mathrm{GW}} h 2 Ω GW γ − A \gamma-A γ − A Figure 1 . In terms of γ \gamma γ A A A γ ≲ 4 \gamma\lesssim4 γ ≲ 4 log 10 A ≲ − 15.0 \log_{10} A\lesssim-15.0 log 10 A ≲ − 15.0 log 10 A ∼ − ( 14.5 ⋯ 14.0 ) \log_{10} A\sim-(14.5\cdots14.0) log 10 A ∼ − ( 14.5 ⋯ 14.0 ) γ ∼ 5 \gamma\sim5 γ ∼ 5
Figure 20. Same as Figure 19 but for the cosmic-string models considered in this work. See Appendix § 17 for details.
In Figure 9 , we present the reconstructed posterior distributions for the cosmic-string tension in our four stable-string models. The left panel of Figure 9 shows our results in the absence of an additional SMBHB contribution to the signal, while the right panel displays the posterior distributions in the combined strings-plus-SMBHB models. In the former case, we find peaked distributions centered around values of the order of log 10 G μ ∼ − ( 10.5 ⋯ 10.0 ) \log_{10}G\mu\sim-(10.5\cdots10.0) log 10 G μ ∼ − ( 10.5 ⋯ 10.0 ) log 10 G μ \log_{10}G\mu log 10 G μ log 10 G μ \log_{10}G\mu log 10 G μ [59 ] at γ ≲ 4.6 \gamma\lesssim4.6 γ ≲ 4.6 G μ G\mu G μ
Figure 9. Reconstructed 1D marginalized posterior distributions for the dimensionless cosmic-string tension parameter G μ G\mu G μ G μ G\mu G μ K K K (10) reaches a value of 1 / 10 1/10 1/10 Figure 28 in Appendix § 18 .§ 4 shows an extended version of this plot that includes the SMBHB parameters A B H B A_{\mathrm{BHB}} A BHB γ B H B \gamma_{\mathrm{BHB}} γ BHB
If we extend the stable-string models by including an SMBHB contribution to the GW signal, the posterior distributions feature not only a peak toward the upper end of the log 10 G μ \log_{10}G\mu log 10 G μ log 10 G μ \log_{10}G\mu log 10 G μ K K K (10) and hence allows us to derive the following upper limits on the cosmic-string tension in the STABLE-C+SMBHB, STABLE-K+SMBHB, STABLE-M+SMBHB, and STABLE-N+SMBHB models, respectively: log 10 G μ < − 9.67 \log_{10}G\mu<-9.67 log 10 G μ < − 9.67 − 9.87 -9.87 − 9.87 − 10.10 -10.10 − 10.10 − 9.71 -9.71 − 9.71 1 / 10 1/10 1/10 G μ G\mu G μ 6 log 10 G μ < − 9.88 \log_{10} G\mu< -9.88 log 10 G μ < − 9.88 − 10.04 -10.04 − 10.04 − 10.26 -10.26 − 10.26 − 9.90 -9.90 − 9.90 Figure 2 , provide a better fit to the data. In the absence of an SMBHB contribution to the GW signal, the Bayes factors for the META-L and META-
LS models are B = 13.4 ± 0.4 \mathcal{B}=13.4\pm0.4 B = 13.4 ± 0.4 B = 21.3 ± 0.8 \mathcal{B}=21.3\pm0.8 B = 21.3 ± 0.8 B = 11.1 ± 0.3 \mathcal{B}=11.1\pm0.3 B = 11.1 ± 0.3 B = 18.9 ± 0.7 \mathcal{B}=18.9\pm0.7 B = 18.9 ± 0.7
The reconstructed 1D and 2D posterior distributions for the two metastable-string models are shown in the corner plots in Figure 10 . From these plots, we read off that the NG15 data prefer values of the decay parameter κ \sqrt{\kappa} κ κ ∼ 8 \sqrt{\kappa}\sim8 κ ∼ 8 log 10 G μ ∼ − ( 7 ⋯ 4 ) \log_{10}G\mu\sim-(7\cdots4) log 10 G μ ∼ − ( 7 ⋯ 4 ) κ ∼ 8 \sqrt{\kappa}\sim8 κ ∼ 8 (50) . This suppression results in a decreased signal strength and a larger spectral tilt. The decrease in signal strength can, however, be compensated by a large cosmic-string tension, such that the final spectrum still has the right amplitude but is now more blue-tilted than in the stable-string case. This interplay of the different factors affecting the GW spectrum from metastable strings effectively leads to a better fit.
Figure 10. Same as in Figure 5 but for the metastable-string models META-L (GW emissions from loops only) in the left panel and META-LS (GW emission from loops and segments) in the right panel. The gray shaded regions are strongly disfavored by the NG15 data, as they result in a K K K 1 / 10 1/10 1/10 (10) ); the teal shaded regions are ruled out by the CMB bound on cosmic strings ([9 , 77 , 218 ] ); and the regions to the right of the yellow contour lines violate the LVK bound in Eq. (24) . Figure 29 in Appendix § 18 .§ 4 shows extended versions of the two plots that include the SMBHB parameters A B H B A_{\mathrm{BHB}} A BHB γ B H B \gamma_{\mathrm{BHB}} γ BHB
Figure 29. Same as Figure 10 but including the SMBHB parameters A B H B A_{\mathrm{BHB}} A BHB γ B H B \gamma_{\mathrm{BHB}} γ BHB
Let us comment on a few characteristic features of the posterior distributions shown in Figure 10 . First, we find that the log 10 G μ \log_{10}G\mu log 10 G μ log 10 G μ \log_{10} G\mu log 10 G μ log 10 G μ ∼ − ( 11 ⋯ 10 ) \log_{10} G\mu\sim-(11\cdots10) log 10 G μ ∼ − ( 11 ⋯ 10 ) κ ≳ 9 \sqrt{\kappa} \gtrsim9 κ ≳ 9 log 10 G μ \log_{10} G\mu log 10 G μ log 10 G μ ∼ − 5 \log_{10} G\mu\sim-5 log 10 G μ ∼ − 5 log 10 G μ \log_{10} G\mu log 10 G μ log 10 G μ ∼ − 5 \log_{10} G\mu\sim-5 log 10 G μ ∼ − 5 (46) , which ensures that no more new loops are formed during the decay regime of the network. Because of this Heaviside theta function, the loop contribution to the GW spectrum moves to frequencies above the PTA band if we raise log 10 G μ \log_{10} G\mu log 10 G μ log 10 G μ ∼ − 5 \log_{10} G\mu\sim-5 log 10 G μ ∼ − 5 log 10 G μ \log_{10} G\mu log 10 G μ log 10 G μ \log_{10} G\mu log 10 G μ Figure 10 can be attributed to the fact that, in the case of the metastable-string models, we have to work with tabulated data for the GW spectrum, based on the numerical evaluation of Eq. (40) .
We also use the corner plots in Figure 10 to highlight the relevant bounds on the parameter space spanned by log 10 G μ \log_{10} G\mu log 10 G μ κ \sqrt{\kappa} κ K K K log 10 G μ ∼ − 7 \log_{10} G\mu\sim-7 log 10 G μ ∼ − 7 [9 , 77 , 218 ] ). In order to derive the CMB bound, we estimate that a decaying network completely disappears because of GW emission at times around t e ≃ ( 2 / ( Γ G μ ) ) 1 / 2 t s t_e \simeq(2/(\Gamma G\mu))^{1/2}t_s t e ≃ ( 2/ ( Γ G μ ) ) 1/2 t s (50) becomes large. The teal shaded areas in Figure 10 indicate where the conditions log 10 G μ > − 7 \log_{10} G\mu> -7 log 10 G μ > − 7 t e > t r e c ≃ 370,000 y r t_e > t_{\mathrm{rec}} \simeq370{,}000\,\mathrm{yr} t e > t rec ≃ 370 , 000 yr
Last but not least, the yellow solid lines represent the LVK bound in Eq. (24) , which translates to the upper limit G μ ≲ 2 × 10 − 7 G\mu\lesssim2 \times10^{-7} G μ ≲ 2 × 1 0 − 7 7
The LVK bound on the amplitude of the stochastic GWB appears to be in tension with most of the parameter space preferred by the NANOGrav data. In particular, the 68% credible regions in the log 10 G μ − κ \log_{10} G\mu-\sqrt{\kappa} log 10 G μ − κ Figure 10 relies on an extrapolation of the GW signal across 10 orders of magnitude in frequency, from PTA frequencies, f ∼ f e w × 10 − 9 f \sim\mathrm{few} \times10^{-9} f ∼ few × 1 0 − 9 f ∼ f e w × 10 f \sim\mathrm{few} \times10 f ∼ few × 10 [22 ] ), which would suppress the GW signal from metastable strings at high frequencies and thus allow for the possibility of cosmic-string tensions as large as those in Figure 10 . On the other hand, if we trust the extrapolation to higher frequencies, we conclude that some parts of the 95% credible regions in Figure 10 are in fact not ruled out. In this case, metastable strings with κ ∼ 8 \sqrt{\kappa} \sim8 κ ∼ 8 log 10 G μ ∼ − 7 \log_{10} G\mu\sim-7 log 10 G μ ∼ − 7 log 10 G μ ∼ − 7 \log_{10} G\mu\sim-7 log 10 G μ ∼ − 7
Finally, we discuss our results for cosmic superstrings, for which we obtain the largest Bayes factors among all cosmic-string models considered in this work: B = 46 ± 2 \mathcal{B}=46\pm2 B = 46 ± 2 B = 37 ± 2 \mathcal{B}=37\pm2 B = 37 ± 2
Superstrings result in a good fit of the data because, unlike ordinary stable strings, they permit a GW signal that is neither too weak nor too flat. To see this, note that small cosmic-string tensions, log 10 G μ ≲ − 11 \log_{10} G\mu\lesssim-11 log 10 G μ ≲ − 11 h 2 Ω G W h^2\Omega_{\mathrm{GW}} h 2 Ω GW log 10 G μ ≲ − 11 \log_{10} G\mu\lesssim-11 log 10 G μ ≲ − 11 1 / P 1/P 1/ P (52) . By choosing P P P Figure 11 . The 2D posterior distribution for log 10 G μ \log_{10} G\mu log 10 G μ P P P
Figure 11. Same as in Figure 10 , but for the SUPER model. Figure 30 in Appendix § 18 .§ 4 shows an extended version of this plot that includes the SMBHB parameters A B H B A_{\mathrm{BHB}} A BHB γ B H B \gamma_{\mathrm{BHB}} γ BHB
Figure 30. Same as Figure 11 but including the SMBHB parameters A B H B A_{\mathrm{BHB}} A BHB γ B H B \gamma_{\mathrm{BHB}} γ BHB
The highest posterior density is achieved at small intercommutation probabilities and cosmic-string tensions, log 10 P ∼ − 3 \log_{10} P \sim-3 log 10 P ∼ − 3 log 10 G μ ∼ − 12 \log_{10} G\mu\sim-12 log 10 G μ ∼ − 12 log 10 P = − 3 \log_{10} P=-3 log 10 P = − 3 P P P (24) (see the yellow solid line in Figure 11 ) but may be viable assuming a nonstandard expansion history at high temperatures. Similarly to what is found in Figure 10 , we highlight the region ruled out by the NG15 data because it leads to a K K K log 10 P \log_{10} P log 10 P [77 ] . Existing upper limits are, however, located at larger G μ G\mu G μ Figure 11 .
In conclusion, we stress again that the GW spectrum from cosmic superstrings is not well understood at present, and more work is needed to arrive at a reliable prediction. Such an effort will be important to confirm that GWs from cosmic superstrings do indeed represent a realistic and viable interpretation of the PTA signal.
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 [188 ] . 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 H^{-3} H − 3 A ∼ O ( 1 ) A \sim O(1) A ∼ O ( 1 ) [254 ] ; Hindmarsh [151 ] (2003); Garagounis & Hindmarsh [128 ] ). In this regime, the energy density stored in domain walls is given by
ρ D W ∼ A σ H (53) \rho_{\mathrm{DW}} \sim A \sigma H \tag*{(53)} ρ DW ∼ A σ H (53) where A A A σ \sigma σ ρ c ∝ H 2 \rho_{\mathrm{c}} \propto H^2 ρ c ∝ H 2 Ω D W = ρ D W / ρ c \Omega_{\mathrm{DW}} = \rho_{\mathrm{DW}}/\rho_{\mathrm{c}} Ω DW = ρ DW / ρ c Ω D W ∝ H − 1 \Omega_{\mathrm{DW}} \propto H^{-1} Ω DW ∝ H − 1 [320 ] . A possible solution to this problem is to have domain walls decay at some temperature T ∗ T_* T ∗ [297 , 132 , 207 ] .
During the scaling regime, domain walls continuously change their configuration and shrink in order to maintain the scaling relation in Eq. (53) . While these processes 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 [297 , 253 , 135 , 76 , 159 , 184 ]
h 2 Ω G W ( f ) = 3 32 π D ϵ ~ α ∗ 2 S ( f / f p ) (54) h^2\Omega_{\mathrm{GW}}(f) = \frac{3}{32\pi}\mathcal{D}\,\tilde{\epsilon}\,\alpha_*^2 S(f/f_p) \tag*{(54)} h 2 Ω GW ( f ) = 32 π 3 D ϵ ~ α ∗ 2 S ( f / f p ) (54) where D \mathcal{D} D (35) , ϵ ~ = 0.7 \tilde{\epsilon} = 0.7 ϵ ~ = 0.7 α ∗ \alpha_* α ∗ T ∗ T_* T ∗
α ∗ ≡ Ω D W ( T ∗ ) = ρ D W 3 H 2 M P l 2 ∣ T = T ∗ (55) \alpha_* \equiv\Omega_{\mathrm{DW}}(T_*) = \left.\frac{\rho_{\mathrm{DW}}}{3H^2 M_{\mathrm{Pl}}^2}\right|_{T=T_*} \tag*{(55)} α ∗ ≡ Ω DW ( T ∗ ) = 3 H 2 M Pl 2 ρ DW T = T ∗ (55) The 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 ∗ H_* H ∗ T ∗ T_* T ∗
f p = 1.14 n H z ( 10.75 g ∗ , s ) 1 / 3 ( g ∗ 10.75 ) 1 / 2 ( T ∗ 10 M e V ) (56) f_p = 1.14\,\mathrm{nHz}\left(\frac{10.75}{g_{*,s}}\right)^{1/3}\left(\frac{g_*}{10.75}\right)^{1/2}\left(\frac{T_*}{10\,\mathrm{MeV}}\right) \tag*{(56)} f p = 1.14 nHz ( g ∗ , s 10.75 ) 1/3 ( 10.75 g ∗ ) 1/2 ( 10 MeV T ∗ ) (56) with both g ∗ g_* g ∗ g ∗ , s g_{*,s} g ∗ , s T ∗ T_* T ∗ S ( x ) S(x) S ( x ) § 7.1 .8
S ( x ) = ( a + b ) c ( b x − a / c + a x b / c ) c (57) S(x) = \frac{(a+b)^c}{(bx^{-a/c}+ax^{b/c})^c} \tag*{(57)} S ( x ) = ( b x − a / c + a x b / c ) c ( a + b ) c (57) Causality fixes a = 3 a=3 a = 3 b ≃ c ≃ 1 b \simeq c \simeq1 b ≃ c ≃ 1 [122 ] , we allow b b b c c c b ∈ [ 0.5 , 1 ] b \in[0.5,1] b ∈ [ 0.5 , 1 ] c ∈ [ 0.3 , 3 ] c \in[0.3,3] c ∈ [ 0.3 , 3 ]
In our analysis we consider two possible decay channels for the domain wall network: dark radiation (DW-DR) and SM particles (DW-SM). In the DW-DR model, instead of using α ∗ \alpha_* α ∗ ρ D R \rho_{\mathrm{DR}} ρ DR Δ N e f f \Delta N_{\mathrm{eff}} Δ N eff § 5.1 and which is bounded from above by BBN and CMB observations, Δ N e f f max ∼ f e w × 0.1 \Delta N_{\mathrm{eff}}^{\max} \sim\mathrm{few} \times0.1 Δ N eff m a x ∼ few × 0.1 [250 Yeh] et al. (2021). In our MCMC analysis of the DW-DR model, we follow Ferreira et al. [122 ] and use Δ N e f f max = 0.39 \Delta N_{\mathrm{eff}}^{\max}=0.39 Δ N eff m a x = 0.39 Δ N e f f \Delta N_{\mathrm{eff}} Δ N eff T ∗ T_{*} T ∗ Δ N e f f \Delta N_{\mathrm{eff}} Δ N eff α ∗ \alpha_{*} α ∗
Δ N e f f ≃ 0.62 ( 10.75 g ∗ ) 1 / 3 ( g ∗ , s g ∗ ) ( α ∗ 0.1 ) , (58) \Delta N_{\mathrm{eff}} \simeq0.62 \left(\frac{10.75}{g_{*}}\right)^{1/3}\left(\frac{g_{*,s}}{g_{*}}\right)\left(\frac{\alpha_{*}}{0.1}\right), \tag*{(58)} Δ N eff ≃ 0.62 ( g ∗ 10.75 ) 1/3 ( g ∗ g ∗ , s ) ( 0.1 α ∗ ) , (58) where, as before, both g ∗ g_{*} g ∗ g ∗ , s g_{*,s} g ∗ , s T ∗ T_{*} T ∗ T ∗ ≳ 2.7 T_{*} \gtrsim2.7 T ∗ ≳ 2.7 [174 ] ; Bai & Korwar [43 ] ) for any detectable value of α ∗ \alpha_{*} α ∗ [122 ] , we also impose α ∗ < 0.3 \alpha_{*}<0.3 α ∗ < 0.3 Δ N e f f \Delta N_{\mathrm{eff}} Δ N eff
Results and Discussion
The reconstructed posterior distributions for the parameters T ∗ T_{*} T ∗ α ∗ \alpha_{*} α ∗ T ∗ T_{*} T ∗ Δ N e f f \Delta N_{\mathrm{eff}} Δ N eff Figure 12 , for both the case where the domain walls are assumed to be the only source of GWs (blue contours) and the scenario where we consider the superposition of the domain wall and SMBHB signals (red contours). Full corner plots including the posterior distributions of the spectral shape parameters b b b c c c Figure 31 in Appendix C.5.
Figure 12. Same as in Figure 5 , but for domain walls decaying to SM particles (left panel) and a dark sector (right panel). Figure 31 in Appendix C.5 shows extended versions of the two plots that include the spectral shape parameters b b b c c c A B H B A_{\mathrm{BHB}} A BHB γ B H B \gamma_{\mathrm{BHB}} γ BHB
Figure 31. Same as Figure 12 but including the spectral shape parameters b b b c c c A B H B A_{\mathrm{BHB}} A BHB γ B H B \gamma_{\mathrm{BHB}} γ BHB
For both the DW-SM and DW-DR models, with and without the inclusion of the SMBHB signal, we find that the GWB produced by domain walls peaks around 10 − 8 10^{-8} 1 0 − 8 Figure 3 and Figure 19 ). Specifically, for the DW-SM (DW-DR) model, we find that T ∗ ∈ [ 110 , 275 ] T_{*} \in[110, 275] T ∗ ∈ [ 110 , 275 ] T ∗ ∈ [ 76 , 505 ] T_{*} \in[76, 505] T ∗ ∈ [ 76 , 505 ] T ∗ ∈ [ 108 , 309 ] T_{*} \in[108, 309] T ∗ ∈ [ 108 , 309 ] T ∗ ∈ [ 67 , 843 ] T_{*} \in[67, 843] T ∗ ∈ [ 67 , 843 ] T ∗ T_{*} T ∗
For the DW-DR model, the posterior distribution for Δ N e f f \Delta N_{\mathrm{eff}} Δ N eff Δ N e f f ≳ 0.32 \Delta N_{\mathrm{eff}} \gtrsim0.32 Δ N eff ≳ 0.32 Δ N e f f ≳ 0.25 \Delta N_{\mathrm{eff}} \gtrsim0.25 Δ N eff ≳ 0.25 Δ N e f f \Delta N_{\mathrm{eff}} Δ N eff Δ N e f f ≳ 0.23 \Delta N_{\mathrm{eff}} \gtrsim0.23 Δ N eff ≳ 0.23 Δ N e f f \Delta N_{\mathrm{eff}} Δ N eff Δ N e f f \Delta N_{\mathrm{eff}} Δ N eff Δ N e f f \Delta N_{\mathrm{eff}} Δ N eff [2 ] ), which promises to be sensitive to Δ N e f f \Delta N_{\mathrm{eff}} Δ N eff Δ N e f f ≃ 0.06 \Delta N_{\mathrm{eff}} \simeq0.06 Δ N eff ≃ 0.06
For the DW-SM model, we find α ∗ ∈ [ 0.080 , 0.19 ] \alpha_\ast\in[0.080, 0.19] α ∗ ∈ [ 0.080 , 0.19 ] α ∗ ∈ [ 0.053 , 0.35 ] \alpha_\ast\in[0.053, 0.35] α ∗ ∈ [ 0.053 , 0.35 ] α ∗ \alpha_\ast α ∗ α ∗ ∈ [ 0.079 , 0.22 ] \alpha_\ast\in[0.079, 0.22] α ∗ ∈ [ 0.079 , 0.22 ] α ∗ ∈ [ 0.047 , 0.61 ] \alpha_\ast\in[0.047, 0.61] α ∗ ∈ [ 0.047 , 0.61 ] T ∗ T_\ast T ∗ α ∗ \alpha_\ast α ∗ α ∗ \alpha_\ast α ∗
For both DW models, we identify regions of the parameter space where the GW signal from domain walls is too strong to be compatible with NG15 data. These regions, shaded in gray in Figure 12 , illustrate for the first time how PTA data can be used to constrain models of domain walls.
Finally, we state the Bayes factors for the model comparison between the domain wall models and the SMBHB reference model. The Bayes factors for DW-SM and DW-DR versus SMBHB are B = 14.8 ± 0.5 \mathcal{B} = 14.8 \pm0.5 B = 14.8 ± 0.5 B = 1.62 ± 0.05 \mathcal{B} = 1.62 \pm0.05 B = 1.62 ± 0.05 B = 21.1 ± 0.9 \mathcal{B} = 21.1 \pm0.9 B = 21.1 ± 0.9 B = 2.53 ± 0.10 \mathcal{B} = 2.53 \pm0.10 B = 2.53 ± 0.10 § 7.1 : adding the SMBHB contribution to the GWB signal results in a plateau region in the posterior distribution of the DW-DR model parameters that manifests itself as part of the 95% credible region in Figure 12 .
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 h h
Ultralight dark matter
MODEL DESCRIPTION
While DM constitutes roughly 27% ([324 ] ) 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 cosmological 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 ([286 , 98 , 270 , 24 , 92 ] ). 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 Λ \Lambda Λ [63 ] ). However, too much suppression on large scales would be in conflict with CMB measurements ([161 , 162 , 163 ] ), which sets a lower bound on the DM mass, 10 − 24 e V < m ϕ 10^{-24}\,\mathrm{eV} < m_\phi 1 0 − 24 eV < m ϕ f f f 2 π f ∼ m ϕ 2\pi f \sim m_\phi 2 π f ∼ m ϕ 10 − 23 e V ≲ m ϕ ≲ 10 − 20 e V 10^{-23}\,\mathrm{eV} \lesssim m_\phi\lesssim10^{-20}\,\mathrm{eV} 1 0 − 23 eV ≲ m ϕ ≲ 1 0 − 20 eV
Other astrophysical constraints, such as measurements of the Lyman-α \alpha α [32 , 190 , 264 ] (Iršič et al. 2017)), galactic subhalo mass functions ([272 , 45 , 237 ] ), and stellar kinematics ([89 ] ), can push this bound further up, ranging from 10 − 21 e V 10^{-21}\,\mathrm{eV} 1 0 − 21 eV 10 − 19 e V 10^{-19}\,\mathrm{eV} 1 0 − 19 eV [33 , 281 , 288 ] ) pushes the bound up to 10 − 21 … 10 − 17 e V 10^{-21}\ldots10^{-17}\,\mathrm{eV} 1 0 − 21 … 1 0 − 17 eV 10 − 23 e V ≲ m ϕ ≲ 10 − 20 e V 10^{-23}\,\mathrm{eV} \lesssim m_\phi\lesssim10^{-20}\,\mathrm{eV} 1 0 − 23 eV ≲ m ϕ ≲ 1 0 − 20 eV [321 , 322 ] ) 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.
While there is a large phenomenology of ULDM signals that can be produced in PTAs, the timing residuals for the I I I h I h_I h I
h I ( t ) = ∑ i A E , I i ( x E ) sin ( ω t + γ E i ) + A P , I i ( x P , I ) sin ( ω t + γ P , I i ) , (59) h_I(t) = \sum_i A^i_{E,I}(x_{E}) \sin(\omega t + \gamma^i_{E}) + A^i_{P,I}(x_{P,I}) \sin(\omega t + \gamma^i_{P,I}), \tag*{(59)} h I ( t ) = i ∑ A E , I i ( x E ) sin ( ω t + γ E i ) + A P , I i ( x P , I ) sin ( ω t + γ P , I i ) , (59) where i i i 9 ω \omega ω A E , I i ( x E ) A^{i}_{E,I}(\mathbf{x}_{E}) A E , I i ( x E ) γ E i \gamma^{i}_{E} γ E i A P , I i ( x P , I ) A^{i}_{P,I}(\mathbf{x}_{P,I}) A P , I i ( x P , I ) γ P , I i \gamma^{i}_{P,I} γ P , I i
γ E i = α i ( x E ) , (60) \gamma^{i}_{E} = \alpha^{i}(\mathbf{x}_{E}), \tag*{(60)} γ E i = α i ( x E ) , (60) γ P , I i = ω ∣ x E − x P , I ∣ + α i ( x P , I ) , (61) \gamma^{i}_{P,I} = \omega\lvert\mathbf{x}_{E}-\mathbf{x}_{P,I}\rvert+\alpha^{i}(\mathbf{x}_{P,I}), \tag*{(61)} γ P , I i = ω ∣ x E − x P , I ∣ + α i ( x P , I ) , (61) where α i ( x ) \alpha^{i}(\mathbf{x}) α i ( x ) x \mathbf{x} x γ P , I i \gamma^{i}_{P,I} γ P , I i
An important scale in understanding these signals is given by the correlation length of the ULDM field, ℓ c \ell_{c} ℓ c
ℓ c ≃ 2 π m ϕ v ϕ ∼ 0.4 k p c ( 10 − 22 e V m ϕ ) , (62) \ell_{c} \simeq\frac{2\pi}{m_{\phi}v_{\phi}} \sim0.4\ \mathrm{kpc}\left(\frac{10^{-22}\ \mathrm{eV}}{m_{\phi}}\right), \tag*{(62)} ℓ c ≃ m ϕ v ϕ 2 π ∼ 0.4 kpc ( m ϕ 1 0 − 22 eV ) , (62) where v ϕ ∼ 10 − 3 v_{\phi}\sim10^{-3} v ϕ ∼ 1 0 − 3 ℓ c ≫ ∣ x E − x P , I ∣ \ell_{c}\gg\lvert\mathbf{x}_{E}-\mathbf{x}_{P,I}\rvert ℓ c ≫ ∣ x E − x P , I ∣
be position independent:
A E , I i ( x E ) → A E i , A P , I i ( x P , I ) → A P , I i , α i ( x ) → α i . (63) A^{i}_{E,I}(\mathbf{x}_{E}) \to A^{i}_{E},\quad A^{i}_{P,I}(\mathbf{x}_{P,I}) \to A^{i}_{P,I},\quad\alpha^{i}(\mathbf{x}) \to\alpha^{i}. \tag*{(63)} A E , I i ( x E ) → A E i , A P , I i ( x P , I ) → A P , I i , α i ( x ) → α i . (63) Generally, 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 A_{P} A P m ϕ ∣ x E − x P , I ∣ ≫ 1 m_{\phi}\lvert\mathbf{x}_{E}-\mathbf{x}_{P,I}\rvert\gg1 m ϕ ∣ x E − x P , I ∣ ≫ 1 γ P , I \gamma_{P,I} γ P , I γ E \gamma_{E} γ E ℓ c ≪ ∣ x E − x P , I ∣ \ell_{c}\ll\lvert\mathbf{x}_{E}-\mathbf{x}_{P,I}\rvert ℓ c ≪ ∣ x E − x P , I ∣
In 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:
Metric fluctuations (Khmelnitsky & Rubakov [187 ] ; Porayko & Postnov [251 ] ; Porayko et al. [252 ] ; Kato & Soda [182 ] ; Nomura et al. [243 ] ; Unal et al. [289 ] ; Wu et al. (2022)): Oscillations in the ULDM field generate fluctuations in the local stress-energy tensor, T μ ν T_{\mu\nu} T μν T μ ν T_{\mu\nu} T μν
specific 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. [243 ] ), 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. [243 ] ).
Doppler–U ( 1 ) U(1) U ( 1 ) ([139 , 312 ] ): 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 U(1)_B U ( 1 ) B Q B = N B × q B Q_B = N_B \times q_B Q B = N B × q B N B N_B N B q B q_B q B U ( 1 ) U(1) U ( 1 ) B − L B-L B − L ϕ n n ˉ \phi n\bar{n} ϕ n n ˉ ϕ \phi ϕ n n n ∇ ϕ \nabla\phi ∇ ϕ [139 ] ), which causes them to be velocity suppressed compared to the vector DM scenario.
The 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,
L ⊃ ϕ Λ [ d γ 4 e 2 F μ ν F μ ν + d g β 3 2 g 3 G μ ν A G A μ ν − ∑ f = e , μ d f m f f ˉ f − ∑ q = u , d ( d q + γ q d g ) m q q ˉ q ] (64) \mathcal{L} \supset\frac{\phi}{\Lambda}\left[\frac{d_\gamma}{4e^2}F_{\mu\nu}F^{\mu\nu}+\frac{d_g\beta_3}{2g_3}G^A_{\mu\nu}G_A^{\mu\nu}-\sum_{f=e,\mu}d_fm_f\bar{f}f-\sum_{q=u,d}(d_q+\gamma_qd_g)m_q\bar{q}q\right] \tag*{(64)} L ⊃ Λ ϕ 4 e 2 d γ F μν F μν + 2 g 3 d g β 3 G μν A G A μν − f = e , μ ∑ d f m f f ˉ f − q = u , d ∑ ( d q + γ q d g ) m q q ˉ q (64) where F μ ν F_{\mu\nu} F μν G μ ν G_{\mu\nu} G μν Λ = M P l / 4 π \Lambda= M_{\mathrm{Pl}}/\sqrt{4\pi} Λ = M Pl / 4 π β 3 \beta_3 β 3 γ q \gamma_q γ q m f m_f m f d d d
Pulsar spin fluctuations (Kaplan et al. [181 ] ): Particle mass fluctuations change the moment of inertia of the pulsar, which, by conservation of angular momentum, leads to pulsar spin fluctuations. The amplitude of the signal is given in Table 1. The γ \gamma γ γ e \gamma_e γ e γ m \gamma_m γ m d m ^ = ( m u d u + m d d d ) / ( m u + m d ) d_{\hat{m}}=(m_ud_u+m_dd_d)/(m_u+m_d) d m ^ = ( m u d u + m d d d ) / ( m u + m d ) [181 ] , we employ in our analysis
γ g = − 5 , γ m ^ = − 2.4 × 10 − 1 , γ μ = 2 × 10 − 3 , γ e = 1.7 × 10 − 5 , (65) \begin{aligned}
\gamma_g&=-5, & \gamma_{\hat{m}}&=-2.4\times10^{-1},\\
\gamma_\mu&=2\times10^{-3}, & \gamma_e&=1.7\times10^{-5}, \tag*{(65)}
\end{aligned} γ g γ μ = − 5 , = 2 × 1 0 − 3 , γ m ^ γ e = − 2.4 × 1 0 − 1 , = 1.7 × 1 0 − 5 , (65) for each pulsar, which assumes the simplest possible model for the pulsars’ moment of inertia.
Reference clock shifts ([139 , 181 ] ): PTA timing residuals are measured with respect to a collection of mostly cesium (McCarthy [227 ] ) 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 γ \gamma γ [181 ] , we take these parameters to be
γ g = 1 γ m ^ = 0.158 γ e = 2 γ γ = 4.83 , (66) \begin{aligned}
\gamma_g&=1 & \gamma_{\hat{m}}&=0.158\\
\gamma_e&=2 & \gamma_\gamma&=4.83, \tag*{(66)}
\end{aligned} γ g γ e = 1 = 2 γ m ^ γ γ = 0.158 = 4.83 , (66) assuming cesium clocks to be used as a reference.
Results and Discussion
We search for ULDM signals by performing a variety of Bayesian fit analyses on the NG15 data set. The priors for the ULDM model parameters are summarized in Table 3, while the priors for the intrinsic red-noise parameters are reported in Table 2. Since in this analysis we want to remain agnostic about the origin of the GWB, we model the GWB as a power law and allow the values of the amplitude and spectral index to float within the following prior ranges: log 10 A G W B ∈ [ − 18 , − 11 ] \log_{10} A_{\mathrm{GWB}}\in[-18,-11] log 10 A GWB ∈ [ − 18 , − 11 ] γ G W B ∈ [ 0 , 7 ] \gamma_{\mathrm{GWB}}\in[0,7] γ GWB ∈ [ 0 , 7 ]
Over most of the ULDM mass range, 10 − 24 e V ≲ m ϕ ≲ 10 − 20 e V 10^{-24}\,\mathrm{eV}\lesssim m_\phi\lesssim10^{-20}\,\mathrm{eV} 1 0 − 24 eV ≲ m ϕ ≲ 1 0 − 20 eV B ∼ 1 \mathcal{B}\sim1 B ∼ 1 f ≃ 4 n H z f\simeq4\,\mathrm{nHz} f ≃ 4 nHz f ∈ [ 3.05 , 6.09 ] n H z f\in[3.05,6.09]\,\mathrm{nHz} f ∈ [ 3.05 , 6.09 ] nHz B ∼ 2 \mathcal{B}\sim2 B ∼ 2 B ∼ 1.5 \mathcal{B}\sim1.5 B ∼ 1.5 m ϕ ∼ 2 × 10 − 23 e V m_{\phi} \sim2 \times10^{-23}\,\mathrm{eV} m ϕ ∼ 2 × 1 0 − 23 eV α \alpha α [32 , 190 , 264 ] (Iršić et al. 2017)), galactic subhalo mass functions ([272 , 45 , 237 ] ), and stellar kinematics ([89 ] ).
Without convincing evidence for a signal, we compute constraints on the ULDM model parameters, shown in Figure 13 -Figure 15 . All constraints are the 95th percentile of the marginalized posterior distribution for the parameter on the vertical axis. The curves labeled “(un)correlated” correspond to the analysis done in the (un)correlated limit, discussed in the previous section.
Figure 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é et al. [51 ] ), and the constraints previously set by the PPTA Collaboration are reported in yellow (Xue et al. [312 ] ).
In Figure 13 , we show the constraints on the local ULDM energy density that can be derived assuming only gravitational coupling between the ULDM and SM fields. The strongest bounds are obtained in the mass range m ϕ ≲ 10 − 23 e V m_{\phi} \lesssim10^{-23}\,\mathrm{eV} m ϕ ≲ 1 0 − 23 eV m ϕ = 10 − 24 e V m_{\phi} = 10^{-24}\,\mathrm{eV} m ϕ = 1 0 − 24 eV m ∼ 10 − 26 e V m \sim10^{-26}\,\mathrm{eV} m ∼ 1 0 − 26 eV 1 / m ϕ 1/m_{\phi} 1/ m ϕ [187 , 289 ] ). While future PTA analyses will be able to improve on these constraints, constraining ULDM with m ϕ ≳ 10 − 23 e V m_{\phi} \gtrsim10^{-23}\,\mathrm{eV} m ϕ ≳ 1 0 − 23 eV ρ ϕ = 0.4 G e V c m − 3 \rho_{\phi} = 0.4\,\mathrm{GeV\,cm}^{-3} ρ ϕ = 0.4 GeV cm − 3 ρ ϕ l i m ∝ m ϕ 3 \rho_{\phi}^{\mathrm{lim}} \propto m_{\phi}^{3} ρ ϕ lim ∝ m ϕ 3 m ϕ ≳ 1 / T o b s m_{\phi} \gtrsim1/T_{\mathrm{obs}} m ϕ ≳ 1/ T obs
In Figure 14 , we show the constraints for all the d d d (64) . Note that, for each panel, we assume that the parameter constrained is the only nonvanishing one. Overall, the limits are in rough agreement with the projections from Kaplan et al. [181 ] and competitive with laboratory constraints ([149 , 51 , 204 , 101 , 244 ] ). The strongest constraints, relative to laboratory bounds, are for ULDM models coupled to the electron, d e d_e d e d μ d_{\mu} d μ d e d_e d e [294 , 181 ] ). Such an insensitivity is not a problem for the PTA observable, though, since the pulsar phase evolution is not affected in the same way as the atomic energy levels. The lack of laboratory constraints on d μ d_{\mu} d μ d ∝ 1 / m ϕ d \propto1/m_{\phi} d ∝ 1/ m ϕ m ϕ ∼ 10 − 26 e V m_{\phi} \sim10^{-26}\,\mathrm{eV} m ϕ ∼ 1 0 − 26 eV 1 / m ϕ 1/m_{\phi} 1/ m ϕ
Figure 14. The red solid (dashed) lines show our constraints at the 95% credible level on the d d d (64) , derived by searching for an (un)correlated signal. We assume that all reference clocks use the transition between the two hyperfine levels of the 133 C s {}^{133}\mathrm{Cs} 133 Cs d d d [149 ] , “Al/Hg atomic clocks” (turquoise) are from [244 ] , “MICROSCOPE” (teal) are from [51 ] , “H/Si clock shift” (orange) are from [244 ] , and “NS binary system” are from [204 ] and [101 ] .
In Figure 15 , we show constraints on the gauge coupling of models where the ULDM is the gauge boson of either U ( 1 ) B U(1)_B U ( 1 ) B U ( 1 ) B − L U(1)_{B-L} U ( 1 ) B − L [252 ] ). This result is somewhat expected: while NANOGrav observes more pulsars, the average observation time is longer in PPTA, so roughly similar bounds are expected.
The constraints presented in Figure 14 and Figure 15 assume that ULDM makes up the entire DM content of the universe. However, if ULDM is only a subcomponent of the total DM abundance, these constraints can be easily rescaled. Indeed, from Table 1 , we see that the amplitudes for the direct coupling signals scale as d ρ ϕ d\sqrt{\rho_\phi} d ρ ϕ Q V ρ ϕ Q_{\mathrm{V}}\sqrt{\rho_\phi} Q V ρ ϕ f ϕ f_\phi f ϕ f ϕ \sqrt{f_\phi} f ϕ
Table 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 A^{i} A i ω \omega ω ρ ϕ \rho_{\phi} ρ ϕ 3 ^{3} 3 m ϕ m_{\phi} m ϕ g g g Q Q Q y y y d d d (64) , and Λ = M P l / 4 π \Lambda= M_{\mathrm{Pl}}/\sqrt{4\pi} Λ = M Pl / 4 π ϕ ^ ( x ) \hat{\phi}(\mathbf{x}) ϕ ^ ( x ) ϕ ^ ( x ) → ϕ ^ \hat{\phi}(\mathbf{x}) \to\hat{\phi} ϕ ^ ( x ) → ϕ ^ n ^ I \hat{\mathbf{n}}_{I} n ^ I ϵ i \boldsymbol{\epsilon}_{i} ϵ i
Lastly, we comment on the prior choice for ϕ ^ \hat{\phi} ϕ ^ [251 , 182 ] assumed that ϕ ^ \hat{\phi} ϕ ^ [75 ] , this assumption is not appropriate. Since the observation time of PTAs is much smaller than the coherence time of the DM, τ ∼ ( m ϕ v 2 ) − 1 \tau\sim(m_\phi v^2)^{-1} τ ∼ ( m ϕ v 2 ) − 1 ρ ϕ \rho_\phi ρ ϕ [252 ] find the distribution of ϕ ^ \hat{\phi} ϕ ^ [125 ] .
Dark matter substructures
Model description
In the Λ \Lambda Λ k ≃ M p c − 1 k \simeq\mathrm{Mpc}^{-1} k ≃ Mpc − 1
PBHs are perhaps the simplest example of such small-scale DM substructures. They can be formed in inflationary theories that create large density fluctuations on small scales (like the ones described in § 6.1 ). Several studies investigated the possibility of identifying a galactic PBH population by analyzing the Doppler and Shapiro signals they can leave in PTA data ([276 , 278 , 102 , 259 , 209 b] ). In this analysis, we will closely follow the method outlined by Lee et al. [210 ] to constrain the local PBH abundance.10
The 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, τ \tau τ dynamic (static ) if τ \tau τ T o b s T_{\mathrm{obs}} T obs
by ([259 , 209 b] )
h D , s t a ( t ) = A D , s t a γ r 2 t 3 , (67) h_{\mathrm{D,sta}}(t) = \frac{A_{\mathrm{D,sta}}}{\gamma r^{2}} t^{3}, \tag*{(67)} h D , sta ( t ) = γ r 2 A D , sta t 3 , (67) where A D , s t a A_{\mathrm{D,sta}} A D , sta § 18 .§ 11 for more details). In the dynamic limit, and assuming that the signal is dominated by the closest transiting PBH, we get
h D , d y n ( t ) = A D , d y n ( t − t 0 ) Θ ( t − t 0 ) , (68) h_{\mathrm{D,dyn}}(t) = A_{\mathrm{D,dyn}} (t-t_{0}) \Theta(t-t_{0}), \tag*{(68)} h D , dyn ( t ) = A D , dyn ( t − t 0 ) Θ ( t − t 0 ) , (68) where Θ \Theta Θ A D , d y n A_{\mathrm{D,dyn}} A D , dyn t 0 t_{0} t 0 § 18 .§ 11 for more details).
The 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 leading-order term of this signal can be parameterized as
h S , s t a ( t ) = A S , s t a γ r 2 t 3 , (69) h_{\mathrm{S,sta}}(t) = \frac{A_{\mathrm{S,sta}}}{\gamma r^{2}} t^{3}, \tag*{(69)} h S , sta ( t ) = γ r 2 A S , sta t 3 , (69) where, as for the Doppler case, A S , s t a A_{\mathrm{S,sta}} A S , sta § 18 .§ 11 for more details). In the dynamic limit, there is no simple parametrization of the Shapiro signal; therefore, we do not search for this signal.
Assuming a monochromatic PBH population, our goal is to derive a posterior distribution for the PBH mass fraction, f P B H ≡ Ω P B H / Ω D M f_{\mathrm{PBH}}\equiv\Omega_{\mathrm{PBH}}/\Omega_{\mathrm{DM}} f PBH ≡ Ω PBH / Ω DM M M M p ( f P B H ∣ δ t , M ) p(f_{\mathrm{PBH}}|\delta t,M) p ( f PBH ∣ δ t , M )
For each given value of f P B H f_{\mathrm{PBH}} f PBH M M M [209 ] 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 2.5\times10^3 2.5 × 1 0 3 p ( log 10 A I ∣ f P B H ) p(\log_{10} A_I|f_{\mathrm{PBH}}) p ( log 10 A I ∣ f PBH ) I I I A A A (67) , (68) , and (69) . In Figure 16 we report some of the distributions derived in this way.11
Figure 16. The black solid (dashed) lines show the posterior distributions p ( log 10 A s t a ∣ δ t ) p(\log_{10} A_{\mathrm{sta}}|\delta t) p ( log 10 A sta ∣ δ t ) p ( log 10 A d y n ∣ δ t ) p(\log_{10} A_{\mathrm{dyn}}|\delta t) p ( log 10 A dyn ∣ δ t ) p ( log 10 A ∣ f P B H ) p(\log_{10} A|f_{\mathrm{PBH}}) p ( log 10 A ∣ f PBH ) f P B H f_{\mathrm{PBH}} f PBH M = 10 − 6 M=10^{-6} M = 1 0 − 6 10 − 10 10^{-10} 1 0 − 10 M ⊙ M_\odot M ⊙ M = 10 M ⊙ M=10M_\odot M = 10 M ⊙
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 ) p(\log_{10} A|\delta t) p ( log 10 A ∣ δ t )
p ( log 10 A ∣ δ t ) = ∏ I = 1 N P p ( log 10 A I ∣ δ t ) . (70) p(\log_{10} A|\delta t)=\prod_{I=1}^{N_P}p(\log_{10} A_I|\delta t). \tag*{(70)} p ( log 10 A ∣ δ t ) = I = 1 ∏ N P p ( log 10 A I ∣ δ t ) . (70) Some of the p ( log 10 A I ∣ δ t ) p(\log_{10} A_I|\delta t) p ( log 10 A I ∣ δ t ) Figure 16 .
p ( f P B H ∣ δ t ) = ∏ I = 1 N P ∫ d log 10 A I p ( f P B H ∣ log 10 A I ) p ( log 10 A I ∣ δ t ) ∝ ∏ I = 1 N P ∫ d log 10 A I p ( log 10 A I ∣ f P B H ) p ( log 10 A I ∣ δ t ) (71) p(f_{\mathrm{PBH}}|\delta t)=\prod_{I=1}^{N_P}\int d\log_{10} A_I\,p(f_{\mathrm{PBH}}|\log_{10} A_I)p(\log_{10} A_I|\delta t)\propto\prod_{I=1}^{N_P}\int d\log_{10} A_I\,p(\log_{10} A_I|f_{\mathrm{PBH}})p(\log_{10} A_I|\delta t) \tag*{(71)} p ( f PBH ∣ δ t ) = I = 1 ∏ N P ∫ d log 10 A I p ( f PBH ∣ log 10 A I ) p ( log 10 A I ∣ δ t ) ∝ I = 1 ∏ N P ∫ d log 10 A I p ( log 10 A I ∣ f PBH ) p ( log 10 A I ∣ δ t ) (71) where, in the second step, we used Bayes theorem and assumed uniform priors on log 10 A I \log_{10} A_I log 10 A I f P B H f_{\mathrm{PBH}} f PBH [210 ] .
DM substructures can also possess macroscopic charges and interact with baryonic matter via long-range Yukawa interactions. These interactions can be modeled by a potential of the form
V f i f t h ( r ) = − α ~ G M M P r e − r / λ , (72) V_{\mathrm{fifth}}(r)=-\tilde{\alpha}\frac{GMM_P}{r}e^{-r/\lambda}, \tag*{(72)} V fifth ( r ) = − α ~ r GM M P e − r / λ , (72) where M M M M P M_P M P λ \lambda λ α ~ \tilde{\alpha} α ~ [97 , 304 , 147 , 203 , 141 ] (); Gresham et al. 2018)). These Yukawa interactions can arise from an effective Lagrangian of the form L ⊃ \mathcal{L}\supset L ⊃ g X ϕ X ˉ X + g n ϕ n ˉ n g_X\phi\bar{X}X + g_n\phi\bar{n}n g X ϕ X ˉ X + g n ϕ n ˉ n X X X n n n ϕ \phi ϕ α ~ ≈ g n g X 4 π G m X m n \tilde{\alpha} \approx\frac{g_ng_X}{4\pi Gm_Xm_n} α ~ ≈ 4 π G m X m n g n g X m n m_n m n M ≲ 100 G e V M \lesssim100\ \mathrm{GeV} M ≲ 100 GeV α \alpha α [220 ] ) 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 G e V M \gg100\ \mathrm{GeV} M ≫ 100 GeV X X X
If 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 § 18 .§ 11 ). Therefore, as shown by Gresham et al. [140 ] , following a procedure similar to the one used to constrain the abundance of PBHs, we can constrain the value of the Yukawa parameter, α ~ \tilde{\alpha} α ~ α ~ \tilde{\alpha} α ~ M M M [209 ] 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 ∣ α ~ ) p(\log_{10} A_I|\tilde{\alpha}) p ( log 10 A I ∣ α ~ ) (71) , we can derive p ( α ~ ∣ δ t ) p(\tilde{\alpha}|\delta t) p ( α ~ ∣ δ t ) α ~ \tilde{\alpha} α ~
Results and Discussion
We start by searching for PBH signals on top of a GWB that we model as a power law with amplitude and spectral index allowed to float within the following prior ranges: log 10 A G W B ∈ [ − 11 , − 18 ] \log_{10} A_{\mathrm{GWB}} \in[-11,-18] log 10 A GWB ∈ [ − 11 , − 18 ] γ G W B ∈ [ 0 , 7 ] \gamma_{\mathrm{GWB}} \in[0,7] γ GWB ∈ [ 0 , 7 ] Table 3 .
The 95% upper limits on f P B H f_{\mathrm{PBH}} f PBH Figure 17 . The dynamic Doppler signal is too weak to produce any detectable signal for any of the f P B H f_{\mathrm{PBH}} f PBH f P B H f_{\mathrm{PBH}} f PBH
Figure 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.
real PTA data. As expected, our constraints are much weaker compared to the projections that were derived by Lee et al. [210 ] using mock data and including only white noise. Indeed, as already discussed by Lee et al. [210 ] , the presence of a common red-noise process significantly reduces the sensitivity to PBH signals.
Finally, in Figure 18 we show the constraints on α ~ \tilde{\alpha} α ~ α ~ \tilde{\alpha} α ~ [300 ] ; Shao et al. [277 ] ; Sun et al. [282 ] ) (properly rescaled to take into account the finite range of the interaction (Gresham et al. [140 ] )) 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. [140 ] ), derived from the coldest known NS to date - PSR J2144-3933 (Guillot et al. [144 ] ). And in gray we report the indirect constraints that can be derived by combining the fifth-force constraints on baryon-baryon interactions (Bergé et al. [51 ] ; Fayet [120 ] ), and Bullet Cluster constraints on DM-DM interactions (Spergel & Steinhardt [279 ] ; Kahlhoefer et al. [179 ] ) (see Coskuner et al. [85 ] ; Gresham et al. [140 ] ). We find that the NG15 constraints are competitive with WEP and NS kinetic heating, especially in the M ≳ 10 − 12 M ⊙ M \gtrsim10^{-12} M_{\odot} M ≳ 1 0 − 12 M ⊙ L ⊃ g X ϕ X ˉ X + g n ϕ n ˉ n \mathcal{L} \supset g_X\phi\bar{X}X + g_n\phi\bar{n}n L ⊃ g X ϕ X ˉ X + g n ϕ n ˉ n
Figure 18. The 95% credible level for the fifth-force strength α ~ \tilde{\alpha} α ~ λ = 1 \lambda= 1 λ = 1 λ = 10 − 1 \lambda= 10^{-1} λ = 1 0 − 1 λ = 10 − 2 \lambda= 10^{-2} λ = 1 0 − 2 λ = 10 − 3 \lambda= 10^{-3} λ = 1 0 − 3
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 the signals produced by nonminimal inflation scenarios, scalar-induced GWs, cosmological phase transitions, several cosmic-string models, and domain walls.
For 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 Figure 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.
Despite 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.
Future 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 spectral 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.
It is worth emphasizing that in all parts of our analysis we described cosmological signals using effective parameters, e.g., T ∗ T_{*} T ∗ α ∗ \alpha_{*} α ∗ H ∗ R ∗ H_{*}R_{*} H ∗ R ∗
Along 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 e V ≲ m ϕ ≲ 10 − 20 e V 10^{-23}\,\mathrm{eV} \lesssim m_{\phi} \lesssim10^{-20}\,\mathrm{eV} 1 0 − 23 eV ≲ m ϕ ≲ 1 0 − 20 eV
The 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.
A.Mi. proposed and initiated the project. A.Mi. and K.Sc. coordinated and led the search. The enterprise wrapper, PTArcade [230 ] , 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 K K § 11.1 , which were prepared by T.T.; Figure 18 , which was prepared by V.L.; and Figure 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.
Contributions 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 e f f N_{\rm eff} N eff
Acknowledgments. 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éndez (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.
Part of this work was conducted using the High Performance Computing Cluster PALMA II at the University of Münster (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.
J.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 (IT-1628-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 NAS5-26555. 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 Caltech 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 Eötvös Loránd University, funded by the Eötvös Loránd 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 DoE under contract DE-SC0011632, and by a Simons Investigator award.
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 r e d A_{\mathrm{red}} A red γ r e d \gamma_{\mathrm{red}} γ red A a A_a A a γ a \gamma_a γ a (4) ) as well as for the SMBHB parameters A B H B A_{\mathrm{BHB}} A BHB γ B H B \gamma_{\mathrm{BHB}} γ BHB ( log 10 A B H B , γ B H B ) (\log_{10} A_{\mathrm{BHB}},\gamma_{\mathrm{BHB}}) ( log 10 A BHB , γ BHB )
μ B H B = ( − 15.6 4.7 ) , σ B H B = 10 − 1 × ( 2.8 − 0.026 − 0.026 1.2 ) , (73) \mu_{\mathrm{BHB}} =
\begin{pmatrix}
-15.6 \\
4.7
\end{pmatrix},
\qquad
\sigma_{\mathrm{BHB}} = 10^{-1} \times
\begin{pmatrix}
2.8 & -0.026 \\
-0.026 & 1.2
\end{pmatrix}, \tag*{(73)} μ BHB = ( − 15.6 4.7 ) , σ BHB = 1 0 − 1 × ( 2.8 − 0.026 − 0.026 1.2 ) , (73) which we obtain by fitting the log 10 A B H B \log_{10} A_{\mathrm{BHB}} log 10 A BHB γ B H B \gamma_{\mathrm{BHB}} γ BHB § 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 K K (10) should not be smaller than 1 / 10 1/10 1/10
The Bayes estimator ⟨ θ ⟩ \langle\theta\rangle ⟨ θ ⟩ θ \theta θ P ( θ ∣ D , H ) P(\theta|D,\mathcal{H}) P ( θ ∣ D , H ) P ( θ ∣ D , H ) P(\theta|D,\mathcal{H}) P ( θ ∣ D , H ) σ θ \sigma_\theta σ θ σ θ 2 \sigma_\theta^2 σ θ 2
⟨ θ ⟩ = ∫ d θ θ P ( θ ∣ D , H ) , σ θ 2 = ⟨ θ 2 ⟩ − ⟨ θ ⟩ 2 = [ ∫ d θ θ 2 P ( θ ∣ D , H ) ] − [ ∫ d θ θ P ( θ ∣ D , H ) ] 2 . (74) \langle\theta\rangle= \int d\theta\,\theta\,P(\theta|D,\mathcal{H}), \qquad
\sigma_\theta^2 = \langle\theta^2\rangle-\langle\theta\rangle^2
= \left[\int d\theta\,\theta^2\,P(\theta|D,\mathcal{H})\right]
-\left[\int d\theta\,\theta\,P(\theta|D,\mathcal{H})\right]^2. \tag*{(74)} ⟨ θ ⟩ = ∫ d θ θ P ( θ ∣ D , H ) , σ θ 2 = ⟨ θ 2 ⟩ − ⟨ θ ⟩ 2 = [ ∫ d θ θ 2 P ( θ ∣ D , H ) ] − [ ∫ d θ θ P ( θ ∣ D , H ) ] 2 . (74) In 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 \theta_{\max} θ m a x θ \theta θ P ( θ ∣ D , H ) P(\theta|D,\mathcal{H}) P ( θ ∣ D , H ) θ \theta θ P ( θ ∣ D , H ) P(\theta|D,\mathcal{H}) P ( θ ∣ D , H )
P ( θ max ∣ D , H ) = max θ P ( θ ∣ D , H ) , (75) P(\theta_{\max}|D,\mathcal{H}) = \max_{\theta} P(\theta|D,\mathcal{H}), \tag*{(75)} P ( θ m a x ∣ D , H ) = θ max P ( θ ∣ D , H ) , (75) and the 68% Bayesian credible interval [ θ 68 min , θ 68 max ] [\theta_{68}^{\min},\theta_{68}^{\max}] [ θ 68 m i n , θ 68 m a x ] P ( θ ∣ D , H ) P(\theta|D,\mathcal{H}) P ( θ ∣ D , H )
∫ θ 68 min θ 68 max d θ P ( θ ∣ D , H ) = 0.68 , (76) \int_{\theta_{68}^{\min}}^{\theta_{68}^{\max}} d\theta\,P(\theta|D,\mathcal{H}) = 0.68, \tag*{(76)} ∫ θ 68 m i n θ 68 m a x d θ P ( θ ∣ D , H ) = 0.68 , (76) where P ( θ ∣ D , H ) > P 68 P(\theta|D,\mathcal{H}) > P_{68} P ( θ ∣ D , H ) > P 68 θ ∈ [ θ 68 min , θ 68 max ] \theta\in[\theta_{68}^{\min},\theta_{68}^{\max}] θ ∈ [ θ 68 m i n , θ 68 m a x ] P 68 P_{68} P 68 K K K Table 4 indicates the value θ K \theta_K θ K θ \theta θ K = 1 / 10 K = 1/10 K = 1/10 § 3.2 ),
K ( θ K ) = P ( D ∣ θ K , H ) P ( D ∣ θ 0 , H ) = 1 10 . (77) K(\theta_K) = \frac{P(D|\theta_K,\mathcal{H})}{P(D|\theta_0,\mathcal{H})} = \frac{1}{10}. \tag*{(77)} K ( θ K ) = P ( D ∣ θ 0 , H ) P ( D ∣ θ K , H ) = 10 1 . (77) Note that, unlike all other quantities discussed in this section, the K K K
Table 5.
Table 6.
caption:
:::
Median GW power spectra
In Figure 3 and Figure 4 in the main text, as well as in Figure 19 and Figure 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 Ω G W h^{2}\Omega_{\mathrm{GW}} h 2 Ω GW f = 1 / T o b s f = 1/T_{\mathrm{obs}} f = 1/ T obs f = 30 / T o b s f = 30/T_{\mathrm{obs}} f = 30/ T obs h 2 Ω G W h^{2}\Omega_{\mathrm{GW}} h 2 Ω GW h 2 Ω G W h^{2}\Omega_{\mathrm{GW}} h 2 Ω GW f f f Figure 19 and Figure 20 .
We stress that the median GW spectra in Figure 3 , Figure 4 , Figure 19 , and Figure 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 Figure 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.
Another 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)
individual 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 Figure 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 Figure 11 ). Another example is the median GWB spectrum of the IGW model, which appears to violate the N e f f N_{\rm eff} N eff Figure 4 ). However, the IGW model as a whole is not ruled out, as one can still choose the number of e e e N r h N_{\mathrm{rh}} N rh N e f f N_{\mathrm{eff}} N eff N r h → 0 N_{\mathrm{rh}} \to0 N rh → 0 § 18 .§ 1 and Figure 22 ).
Supplementary material
Cosmic inflation
In Figure 22 , we present constraints on the parameter space of inflationary GWs, i.e., the IGW model discussed in § 5.1 . This parameter space is spanned by four quantities: the tensor-to-scalar ratio r r r n t n_{\mathrm{t}} n t T r h T_{\mathrm{rh}} T rh e e e N r h N_{\mathrm{rh}} N rh n t n_{\mathrm{t}} n t r r r Figure 5 , we are able to eliminate r r r r r r (22) . The effective parameter space spanned by T r h T_{\mathrm{rh}} T rh n t n_{\mathrm{t}} n t N r h N_{\mathrm{rh}} N rh Δ N e f f \Delta N_{\mathrm{eff}} Δ N eff (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 e f f N_{\mathrm{eff}} N eff f e n d f_{\mathrm{end}} f end T r h T_{\mathrm{rh}} T rh H e n d H_{\mathrm{end}} H end N e f f N_{\mathrm{eff}} N eff f e n d max f_{\mathrm{end}}^{\max} f end m a x H e n d max H_{\mathrm{end}}^{\max} H end m a x N r h max N_{\mathrm{rh}}^{\max} N rh m a x e e e f l v k = 25 H z f_{\mathrm{lvk}}=25\,\mathrm{Hz} f lvk = 25 Hz H e n d max H_{\mathrm{end}}^{\max} H end m a x N r h e n d N_{\mathrm{rh}}^{\mathrm{end}} N rh end f e n d f_{\mathrm{end}} f end f e n d max = 20 H z f_{\mathrm{end}}^{\max}=20\,\mathrm{Hz} f end m a x = 20 Hz
Our results for N r h max N_{\mathrm{rh}}^{\max} N rh m a x Figure 22 in terms of black and maroon contour lines, where the black solid lines refer to the N e f f N_{\mathrm{eff}} N eff N r h max N_{\mathrm{rh}}^{\max} N rh m a x N r h max N_{\mathrm{rh}}^{\max} N rh m a x (15) and (16) ),
N r h max = 2 3 ln ( H e n d max H e n d min ) , N r h n a i v e = 2 3 ln ( H e n d n a i v e H e n d min ) , H e n d n a i v e = ( π 2 2 r A s ) 1 / 2 M P l , (78) N_{\rm rh}^{\max} = \frac{2}{3}\ln\left(\frac{H_{\rm end}^{\max}}{H_{\rm end}^{\min}}\right), \qquad N_{\rm rh}^{\rm naive} = \frac{2}{3}\ln\left(\frac{H_{\rm end}^{\rm naive}}{H_{\rm end}^{\min}}\right), \qquad H_{\rm end}^{\rm naive} = \left(\frac{\pi^2}{2}r A_{\rm s}\right)^{1/2} M_{\rm Pl}, \tag*{(78)} N rh m a x = 3 2 ln ( H end m i n H end m a x ) , N rh naive = 3 2 ln ( H end m i n H end naive ) , H end naive = ( 2 π 2 r A s ) 1/2 M Pl , (78) where H e n d min H_{\rm end}^{\min} H end m i n T r h T_{\rm rh} T rh
H e n d min = ( π 2 g ∗ r h 90 ) 1 / 2 T r h 2 M P l . (79) H_{\rm end}^{\min} = \left(\frac{\pi^2 g_*^{\rm rh}}{90}\right)^{1/2}\frac{T_{\rm rh}^2}{M_{\rm Pl}}. \tag*{(79)} H end m i n = ( 90 π 2 g ∗ rh ) 1/2 M Pl T rh 2 . (79) In other words, for given T r h T_{\rm rh} T rh H e n d min H_{\rm end}^{\min} H end m i n N r h = 0 N_{\rm rh}=0 N rh = 0 N r h max N_{\rm rh}^{\max} N rh m a x N r h n a i v e N_{\rm rh}^{\rm naive} N rh naive
Δ N r h = N r h max − N r h n a i v e = 2 3 ln ( H e n d max H e n d n a i v e ) , (80) \Delta N_{\rm rh}=N_{\rm rh}^{\max}-N_{\rm rh}^{\rm naive}=\frac{2}{3}\ln\left(\frac{H_{\rm end}^{\max}}{H_{\rm end}^{\rm naive}}\right), \tag*{(80)} Δ N rh = N rh m a x − N rh naive = 3 2 ln ( H end naive H end m a x ) , (80) which is the quantity shown in the right panel of Figure 22 . In addition to the contour lines for N r h max N_{\rm rh}^{\max} N rh m a x Δ N r h \Delta N_{\rm rh} Δ N rh T r h − n t T_{\rm rh}-n_{\rm t} T rh − n t T r h T_{\rm rh} T rh n t n_{\rm t} n t
In view of the results in Figure 22 , we can draw several conclusions: (i) Both the N e f f N_{\rm eff} N eff T r h ∼ 10 − 3 G e V T_{\rm rh}\sim10^{-3}\ {\rm GeV} T rh ∼ 1 0 − 3 GeV n t ≳ 2.5 n_{\rm t}\gtrsim2.5 n t ≳ 2.5 T r h ∼ 10 3 G e V T_{\rm rh}\sim10^3\ {\rm GeV} T rh ∼ 1 0 3 GeV N r h N_{\rm rh} N rh n t ≳ 1 n_{\rm t}\gtrsim1 n t ≳ 1 f ≫ f r h f\gg f_{\rm rh} f ≫ f rh n t n_{\rm t} n t α \alpha α α ∼ − 1 \alpha\sim-1 α ∼ − 1 α ∼ 0.5 \alpha\sim0.5 α ∼ 0.5 (25) to be valid (see also the discussion in Kuroyanagi et al. [205 ] ). (ii) The LVK bound only depends on T r h T_{\rm rh} T rh f e n d max = 20 H z f_{\rm end}^{\max}=20\ {\rm Hz} f end m a x = 20 Hz r r r n t n_{\rm t} n t N r h N_{\rm rh} N rh N e f f N_{\rm eff} N eff N r h N_{\mathrm{rh}} N rh N e f f N_{\mathrm{eff}} N eff N r h max → ∞ N_{\mathrm{rh}}^{\max} \to\infty N rh m a x → ∞ N e f f N_{\mathrm{eff}} N eff N r h max N_{\mathrm{rh}}^{\max} N rh m a x N e f f N_{\mathrm{eff}} N eff Figure 5 ). In the excluded region, we also find Δ N r h ≪ 0 \Delta N_{\mathrm{rh}} \ll0 Δ N rh ≪ 0 N r h n a i v e N_{\mathrm{rh}}^{\mathrm{naive}} N rh naive N e f f N_{\mathrm{eff}} N eff N r h N_{\mathrm{rh}} N rh N r h max ≫ 1 N_{\mathrm{rh}}^{\max} \gg1 N rh m a x ≫ 1 N r h ∼ O ( 1 ⋯ 10 ) N_{\mathrm{rh}} \sim\mathcal{O}(1 \cdots10) N rh ∼ O ( 1 ⋯ 10 ) Figure 22 can host viable realizations of the igw model. This includes, in particular, scenarios with a low reheating temperature, T r h ∼ 10 − ( 3 ⋯ 0 ) G e V T_{\mathrm{rh}} \sim10^{-(3\cdots0)}\,\mathrm{GeV} T rh ∼ 1 0 − ( 3 ⋯ 0 ) GeV n t ∼ 3 ⋯ 4 n_t \sim3\cdots4 n t ∼ 3 ⋯ 4
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 P B H ≤ 1 f_{\mathrm{PBH}} \le1 f PBH ≤ 1 (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 [255 ] ; Carr [74 ] ) and mainly follow Inomata et al. [168 ] and Ando et al. [28 ] . The main quantity of interest is f P B H f_{\mathrm{PBH}} f PBH
f P B H = ∫ d M M f ( M ) . (81) f_{\mathrm{PBH}} = \int\frac{\mathrm{d}M}{M}\,f(M). \tag*{(81)} f PBH = ∫ M d M f ( M ) . (81) with PBH mass function f f f [28 ] )
f ( M ) ≃ γ 3 / 2 ( β ( M ) 1.6 × 10 − 9 ) ( 10.75 g ∗ ( T ) ) 1 / 4 ( 0.12 Ω D M h 2 ) ( M ⊙ M ) 1 / 2 . (82) f(M) \simeq\gamma^{3/2}\left(\frac{\beta(M)}{1.6 \times10^{-9}}\right)\left(\frac{10.75}{g_*(T)}\right)^{1/4}\left(\frac{0.12}{\Omega_{\mathrm{DM}}h^2}\right)\left(\frac{M_\odot}{M}\right)^{1/2}. \tag*{(82)} f ( M ) ≃ γ 3/2 ( 1.6 × 1 0 − 9 β ( M ) ) ( g ∗ ( T ) 10.75 ) 1/4 ( Ω DM h 2 0.12 ) ( M M ⊙ ) 1/2 . (82) Here M M M T T T γ \gamma γ β ( M ) \beta(M) β ( M ) γ = 0.356 \gamma= 0.356 γ = 0.356 [79 ] , Koike et al. [195 ] and Niemeyer & Jedamzik [242 ] and assume that the mass of the PBH is proportional to the horizon mass when k = a H k = aH k = a H
M = γ M H ≃ γ M ⊙ ( g ∗ ( T ) 10.75 ) − 1 / 6 ( 4.2 × 10 6 M p c − 1 k ) 2 . (83) M = \gamma M_H \simeq\gamma M_\odot\left(\frac{g_*(T)}{10.75}\right)^{-1/6} \left(\frac{4.2 \times10^6\,\mathrm{Mpc}^{-1}}{k}\right)^2 . \tag*{(83)} M = γ M H ≃ γ M ⊙ ( 10.75 g ∗ ( T ) ) − 1/6 ( k 4.2 × 1 0 6 Mpc − 1 ) 2 . (83) We can then express the PBH mass as a function of the temperature using the standard temperature–wavenumber relation in Λ \Lambda Λ β ( M ) \beta(M) β ( M )
β ( M ) = ∫ δ c ∞ d δ 2 π σ ( M ) exp ( − δ 2 2 σ 2 ( M ) ) ≃ σ ( M ) 2 π δ c exp ( − δ c 2 2 σ 2 ( M ) ) , (84) \beta(M) = \int_{\delta_c}^{\infty} \frac{d\delta}{\sqrt{2\pi}\,\sigma(M)} \exp\left(-\frac{\delta^2}{2\sigma^2(M)}\right) \simeq\frac{\sigma(M)}{\sqrt{2\pi}\,\delta_c} \exp\left(-\frac{\delta_c^2}{2\sigma^2(M)}\right), \tag*{(84)} β ( M ) = ∫ δ c ∞ 2 π σ ( M ) d δ exp ( − 2 σ 2 ( M ) δ 2 ) ≃ 2 π δ c σ ( M ) exp ( − 2 σ 2 ( M ) δ c 2 ) , (84) where δ c \delta_c δ c 0.4 ≲ δ c ≲ 0.6 0.4 \lesssim\delta_c \lesssim0.6 0.4 ≲ δ c ≲ 0.6 [236 , 234 ] ), and σ 2 ( M ) \sigma^2(M) σ 2 ( M )
σ 2 ( k ) = 16 81 ∫ d q q ( q k ) 4 W 2 ( q k ) T 2 ( q , k − 1 ) P R ( q ) . (85) \sigma^2(k) = \frac{16}{81} \int\frac{dq}{q} \left(\frac{q}{k}\right)^4 W^2\left(\frac{q}{k}\right) \mathcal{T}^2(q,k^{-1}) \mathcal{P}_{\mathcal{R}}(q) . \tag*{(85)} σ 2 ( k ) = 81 16 ∫ q d q ( k q ) 4 W 2 ( k q ) T 2 ( q , k − 1 ) P R ( q ) . (85) Here W ( x ) W(x) W ( x ) k − 1 k^{-1} k − 1
T ( q , k − 1 ) = 3 [ sin ( q 3 k ) − ( q 3 k ) cos ( q 3 k ) ] / ( q 3 k ) 3 . (86) \mathcal{T}(q,k^{-1}) = 3\left[\sin\left(\frac{q}{\sqrt{3}k}\right) - \left(\frac{q}{\sqrt{3}k}\right)\cos\left(\frac{q}{\sqrt{3}k}\right)\right]\left/\left(\frac{q}{\sqrt{3}k}\right)^3\right. . \tag*{(86)} T ( q , k − 1 ) = 3 [ sin ( 3 k q ) − ( 3 k q ) cos ( 3 k q ) ] / ( 3 k q ) 3 . (86) is the transfer function for the radiation-dominated era. In our work, we choose δ c = 0.45 \delta_c = 0.45 δ c = 0.45 [28 , 301 ] ) and a Gaussian window function W ( x ) = e − x 2 / 2 W(x) = e^{-x^2/2} W ( x ) = e − x 2 /2
Our analysis is subject to a few caveats. The correct value of the critical density δ c \delta_c δ c [28 , 233 , 315 , 112 , 114 , 138 , 234 , 91 ] ), and the nonlinear relation between density perturbations and density contrast ([93 , 317 , 112 ] ). In this work, we use a single value of δ c \delta_c δ c [6 , 113 , 178 , 91 , 235 ] ).
The statistics of the primordial scalar perturbations also strongly affect the abundance of PBHs since, in general, enhanced scalar perturbations come with different levels of non-Gaussianity. In this work, we focus on Gaussian primordial scalar perturbations. With non-Gaussianities, the same PBH abundance can be produced with scalar perturbations whose amplitudes are orders of magnitude smaller than those produced with Gaussian scalar perturbations. Therefore, primordial non-Gaussianity can dramatically affect the PBH bounds ([316 , 239 , 246 , 131 , 126 , 65 , 93 , 317 , 167 ] ). 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 ([21 , 257 , 295 , 94 ] ).
In 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 P B H f_{\rm PBH} f PBH
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 a = 3 a = 3 Figure 25 we report the results of a phase transition analysis where we assume a = 3 a = 3 a = 3 Figure 25 shows the reconstructed posterior distributions for the parameters α ∗ \alpha_\ast α ∗ T ∗ T_\ast T ∗ H ∗ R ∗ H_\ast R_\ast H ∗ R ∗ A B H B A_\mathrm{BHB} A BHB γ B H B \gamma_\mathrm{BHB} γ BHB T ∗ T_\ast T ∗
Figure 25. Same as in Figure 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 a = 3 a = 3
In Figure 26 and Figure 27 , we report the posterior distributions for all the parameters of the phase transition models, including the spectral shape parameters a a a b b b c c c Figure 8 . As expected for the PT-BUBBLE model, the low-frequency slope is peaked around a ∼ 2 a \sim2 a ∼ 2 b b b c c c a a a a = 3 a = 3 a = 3
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 μ G\mu G μ Figure 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 μ G\mu G μ A B H B A_\mathrm{BHB} A BHB γ B H B \gamma_\mathrm{BHB} γ BHB § 8.1 and show the corner plots for these parameters (see Figure 28 ). A notable feature of these corner plots is that the plateau region at small values of G μ G\mu G μ Figure 9 , now also appears in the form of flat directions in the 2D posterior distributions for G μ G\mu G μ A B H B A_\mathrm{BHB} A BHB G μ G\mu G μ γ B H B \gamma_\mathrm{BHB} γ BHB A B H B A_\mathrm{BHB} A BHB γ B H B \gamma_\mathrm{BHB} γ BHB A B H B A_\mathrm{BHB} A BHB γ B H B \gamma_\mathrm{BHB} γ BHB Figure 28 are qualitatively identical and only display slight numerical differences.
Figure 28. Posterior distributions for the STABLE+SMBHB strings models including the SMBHB parameters A B H B A_{\mathrm{BHB}} A BHB γ B H B \gamma_{\mathrm{BHB}} γ BHB log 10 G μ \log_{10} G\mu log 10 G μ
Next, 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 (42) , where we replace the loop number density n l n_l n l n s n_s n s P k P_k P k (41) by P k = 4 / k P_k = 4/k P k = 4/ k [221 ] in the approximation of a straight string segment in Martin & Vilenkin [221 ] . Meanwhile, all relevant expressions for the segment number density n s n_s n s [208 ] and Buchmuller et al. [62 ] , which we shall summarize in this appendix. For more details, we refer to Leblond et al. [208 ] and Buchmuller et al. [62 ] .
First, we consider segments that form when long strings break apart. We denote their number density by n ˉ s \bar{n}_s n ˉ s
n ˉ s ( ℓ , t ) = n ˉ s r r ( ℓ , t ) + n ˉ s r m ( ℓ , t ) + n ˉ s m m ( ℓ , t ) . (87) \bar{n}_s(\ell,t) = \bar{n}_s^{\mathrm{rr}}(\ell,t) + \bar{n}_s^{\mathrm{rm}}(\ell,t) + \bar{n}_s^{\mathrm{mm}}(\ell,t) . \tag*{(87)} n ˉ s ( ℓ , t ) = n ˉ s rr ( ℓ , t ) + n ˉ s rm ( ℓ , t ) + n ˉ s mm ( ℓ , t ) . (87) Here n ˉ s r r ( ℓ , t ) \bar{n}^{\mathrm{rr}}_s(\ell,t) n ˉ s rr ( ℓ , t )
n ˉ s r r ( ℓ , t ) = Θ ( t e q − t ) Θ ( t − t s ) Γ d 2 ( t + t s ) 2 ξ r 2 ( 3 t s ) 1 / 2 exp [ − Γ d ( ℓ ( t + t s ) + 1 2 Γ G μ ( t − t s ) ( t + 3 t s ) ) ] , (88) \bar{n}^{\mathrm{rr}}_s(\ell,t) = \Theta(t_{\mathrm{eq}}-t)\Theta(t-t_s)\frac{\Gamma_d^2(t+t_s)^2}{\xi_r^2(3t_s)^{1/2}}\exp\left[-\Gamma_d\left(\ell(t+t_s)+\frac{1}{2}\Gamma G\mu(t-t_s)(t+3t_s)\right)\right], \tag*{(88)} n ˉ s rr ( ℓ , t ) = Θ ( t eq − t ) Θ ( t − t s ) ξ r 2 ( 3 t s ) 1/2 Γ d 2 ( t + t s ) 2 exp [ − Γ d ( ℓ ( t + t s ) + 2 1 Γ G μ ( t − t s ) ( t + 3 t s ) ) ] , (88) where n ˉ s r m ( ℓ , t ) \bar{n}^{\mathrm{rm}}_s(\ell,t) n ˉ s rm ( ℓ , t )
n ˉ s r m ( ℓ , t ) = Θ ( t − t e q ) Θ ( t e q − t s ) Γ d 2 ξ r 2 ( t e q t ) 2 ( t + t s ) 2 ( t e q 3 t s ) 1 / 2 exp [ − Γ d ( ℓ ( t + t s ) + 1 2 Γ G μ ( t − t s ) ( t + 3 t s ) ) ] . (89) \bar{n}^{\mathrm{rm}}_s(\ell,t) = \Theta(t-t_{\mathrm{eq}})\Theta(t_{\mathrm{eq}}-t_s)\frac{\Gamma_d^2}{\xi_r^2}\left(\frac{t_{\mathrm{eq}}}{t}\right)^2\frac{(t+t_s)^2}{(t_{\mathrm{eq}}^3t_s)^{1/2}}\exp\left[-\Gamma_d\left(\ell(t+t_s)+\frac{1}{2}\Gamma G\mu(t-t_s)(t+3t_s)\right)\right]. \tag*{(89)} n ˉ s rm ( ℓ , t ) = Θ ( t − t eq ) Θ ( t eq − t s ) ξ r 2 Γ d 2 ( t t eq ) 2 ( t eq 3 t s ) 1/2 ( t + t s ) 2 exp [ − Γ d ( ℓ ( t + t s ) + 2 1 Γ G μ ( t − t s ) ( t + 3 t s ) ) ] . (89) and n ˉ s m m ( ℓ , t ) \bar{n}_{s}^{\mathrm{mm}}(\ell,t) n ˉ s mm ( ℓ , t )
n ˉ s m m ( ℓ , t ) = Θ ( t − t s ) Θ ( t s − t e q ) Γ 2 ξ m 2 exp [ − Γ d ( ℓ t + 1 2 Γ G μ ( t − t s ) ( t + t s ) ) ] . (90) \bar{n}_{s}^{\mathrm{mm}}(\ell,t) = \Theta(t-t_s)\Theta(t_s-t_{\mathrm{eq}})\frac{\Gamma^2}{\xi_m^2}\exp\left[-\Gamma_d\left(\ell t+\frac{1}{2}\Gamma G\mu(t-t_s)(t+t_s)\right)\right]. \tag*{(90)} n ˉ s mm ( ℓ , t ) = Θ ( t − t s ) Θ ( t s − t eq ) ξ m 2 Γ 2 exp [ − Γ d ( ℓ t + 2 1 Γ G μ ( t − t s ) ( t + t s ) ) ] . (90) In 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 \Gamma= 50 Γ = 50 ξ r = 0.27 \xi_r = 0.27 ξ r = 0.27 ξ m = 0.63 \xi_m = 0.63 ξ m = 0.63 n ˉ s r r \bar{n}_{s}^{\mathrm{rr}} n ˉ s rr n ˉ s r m \bar{n}_{s}^{\mathrm{rm}} n ˉ s rm n ˉ s m m \bar{n}_{s}^{\mathrm{mm}} n ˉ s mm t e q > t s t_{\mathrm{eq}} > t_s t eq > t s n ˉ s r r \bar{n}_{s}^{\mathrm{rr}} n ˉ s rr n ˉ s r m \bar{n}_{s}^{\mathrm{rm}} n ˉ s rm n ˉ s r r \bar{n}_{s}^{\mathrm{rr}} n ˉ s rr t < t e q t < t_{\mathrm{eq}} t < t eq n ˉ s r m \bar{n}_{s}^{\mathrm{rm}} n ˉ s rm t > t e q t > t_{\mathrm{eq}} t > t eq t e q < t s t_{\mathrm{eq}} < t_s t eq < t s n ˉ s m m \bar{n}_{s}^{\mathrm{mm}} n ˉ s mm t > t s t > t_s t > t s
In 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 \tilde{n}_{s,1} n ~ s , 1 n ~ s , 1 \tilde{n}_{s,1} n ~ s , 1 n ~ s , 1 \tilde{n}_{s,1} n ~ s , 1
n ~ s , 1 ( ℓ , t ) = n ~ s , 1 r r ( ℓ , t ) + n ~ s , 1 r m ( ℓ , t ) + n ~ s , 1 m m ( ℓ , t ) . (91) \tilde{n}_{s,1}(\ell,t)=\tilde{n}_{s,1}^{\mathrm{rr}}(\ell,t)+\tilde{n}_{s,1}^{\mathrm{rm}}(\ell,t)+\tilde{n}_{s,1}^{\mathrm{mm}}(\ell,t). \tag*{(91)} n ~ s , 1 ( ℓ , t ) = n ~ s , 1 rr ( ℓ , t ) + n ~ s , 1 rm ( ℓ , t ) + n ~ s , 1 mm ( ℓ , t ) . (91) Here n ~ s , 1 r r ( ℓ , t ) \tilde{n}_{s,1}^{\mathrm{rr}}(\ell,t) n ~ s , 1 rr ( ℓ , t )
n ~ s , 1 r r ( ℓ , t ) = Θ ( t − t s ) Θ ( t e q − t ) Θ ( t e q − t ∗ ) Γ d [ ℓ ( t − t s ) + 1 2 Γ G μ ( t − t s ) 2 ] n ℓ m e t a ( ℓ , t ) , (92) \tilde{n}_{s,1}^{\mathrm{rr}}(\ell,t)=\Theta(t-t_s)\Theta(t_{\mathrm{eq}}-t)\Theta(t_{\mathrm{eq}}-t_*)\Gamma_d\left[\ell(t-t_s)+\tfrac{1}{2}\Gamma G\mu(t-t_s)^2\right]n_{\ell}^{\mathrm{meta}}(\ell,t), \tag*{(92)} n ~ s , 1 rr ( ℓ , t ) = Θ ( t − t s ) Θ ( t eq − t ) Θ ( t eq − t ∗ ) Γ d [ ℓ ( t − t s ) + 2 1 Γ G μ ( t − t s ) 2 ] n ℓ meta ( ℓ , t ) , (92) where n ~ s , 1 r m ( ℓ , t ) \tilde{n}_{s,1}^{\mathrm{rm}}(\ell,t) n ~ s , 1 rm ( ℓ , t )
n ~ s , 1 r m ( ℓ , t ) = Θ ( t − t s ) Θ ( t − t e q ) Θ ( t e q − t ∗ ) Γ d [ ℓ ( t − t s ) + 1 2 Γ G μ ( t − t s ) 2 ] n ℓ m e t a ( ℓ , t ) , (93) \tilde{n}_{s,1}^{\mathrm{rm}}(\ell,t)=\Theta(t-t_s)\Theta(t-t_{\mathrm{eq}})\Theta(t_{\mathrm{eq}}-t_*)\Gamma_d\left[\ell(t-t_s)+\tfrac{1}{2}\Gamma G\mu(t-t_s)^2\right]n_{\ell}^{\mathrm{meta}}(\ell,t), \tag*{(93)} n ~ s , 1 rm ( ℓ , t ) = Θ ( t − t s ) Θ ( t − t eq ) Θ ( t eq − t ∗ ) Γ d [ ℓ ( t − t s ) + 2 1 Γ G μ ( t − t s ) 2 ] n ℓ meta ( ℓ , t ) , (93) and n ~ s , 1 m m ( ℓ , t ) \tilde{n}_{s,1}^{\mathrm{mm}}(\ell,t) n ~ s , 1 mm ( ℓ , t )
n ~ s , 1 m m ( ℓ , t ) = Θ ( t − t s ) Θ ( t s − t ∗ ) Θ ( t ∗ − t e q ) Γ d e − Γ d [ ℓ ( t − t ∗ ) + 1 2 Γ G μ ( t − t ∗ ) 2 ] t 2 ( ℓ + Γ G μ t ) 2 (94) \tilde{n}_{s,1}^{\mathrm{mm}}(\ell,t)=\Theta(t-t_s)\Theta(t_s-t_*)\Theta(t_*-t_{\mathrm{eq}})\frac{\Gamma_d e^{-\Gamma_d[\ell(t-t_*)+\frac{1}{2}\Gamma G\mu(t-t_*)^2]}}{t^2(\ell+\Gamma G\mu t)^2} \tag*{(94)} n ~ s , 1 mm ( ℓ , t ) = Θ ( t − t s ) Θ ( t s − t ∗ ) Θ ( t ∗ − t eq ) t 2 ( ℓ + Γ G μ t ) 2 Γ d e − Γ d [ ℓ ( t − t ∗ ) + 2 1 Γ G μ ( t − t ∗ ) 2 ] (94) × { 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 ) ] } . (95) \times\left\{0.27\left[\ell(t-t_s)+\tfrac{1}{2}\Gamma G\mu(t-t_s)^2\right]+0.45(\ell+\Gamma G\mu t)^{1+0.31}\left[F_2(t)-F_1(t)-F_2(t_s)+F_1(t_s)\right]\right\}. \tag*{(95)} × { 0.27 [ ℓ ( t − t s ) + 2 1 Γ 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 ) ] } . (95) where F 2 F_2 F 2 F 1 F_1 F 1 2 F 1 {}_2F_1 2 F 1
F 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. (96) F_n(x)={}_2F_1\left(n-\beta,-\beta;n+1-\beta;\frac{\Gamma G\mu}{\ell+\Gamma G\mu t}x\right)\left(\frac{\Gamma G\mu}{\ell+\Gamma G\mu t}\right)^{n-1}\frac{x^{n-\beta}}{n-\beta},\qquad n=1,2,\qquad\beta=0.31. \tag*{(96)} F n ( x ) = 2 F 1 ( n − β , − β ; n + 1 − β ; ℓ + Γ G μ t Γ G μ x ) ( ℓ + Γ G μ t Γ G μ ) n − 1 n − β x n − β , n = 1 , 2 , β = 0.31. (96) These results were derived in Buchmuller et al. [62 ] 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. (94) and (C23) (see Buchmuller et al. [62 ] for more details).
Finally, the expressions in Eqs. (92) , (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 n n ( n + 1 ) (n+1) ( n + 1 ) [62 ] ). 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. [62 ] , 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. (92) , (C21), and (C22)—up to a numerical fudge factor of O ( 10 ) O(10) O ( 10 )
n ~ s r r ( ℓ , t ) → f u d g e × n ~ s , 1 r r ( ℓ , t ) , n ~ s n s ( ℓ , t ) → f u d g e × n ~ s , 1 n s ( ℓ , t ) , n ~ s m m ( ℓ , t ) → f u d g e × n ~ s , 1 m m ( ℓ , t ) . (97) \tilde n_s^{\mathrm{rr}}(\ell,t)\to\mathrm{fudge}\times\tilde n_{s,1}^{\mathrm{rr}}(\ell,t),\qquad
\tilde n_s^{\mathrm{ns}}(\ell,t)\to\mathrm{fudge}\times\tilde n_{s,1}^{\mathrm{ns}}(\ell,t),\qquad
\tilde n_s^{\mathrm{mm}}(\ell,t)\to\mathrm{fudge}\times\tilde n_{s,1}^{\mathrm{mm}}(\ell,t). \tag*{(97)} n ~ s rr ( ℓ , t ) → fudge × n ~ s , 1 rr ( ℓ , t ) , n ~ s ns ( ℓ , t ) → fudge × n ~ s , 1 ns ( ℓ , t ) , n ~ s mm ( ℓ , t ) → fudge × n ~ s , 1 mm ( ℓ , t ) . (97) We stress that the functional dependence on ℓ \ell ℓ t t t (97) , 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.
Domain walls
In Figure 31 , we report the posterior distribution for all the parameters of the domain wall models, including the spectral shape parameters b b b c c c Figure 12 . For both the DW-SM and DW-DR models we find a flat posterior for the high-frequency slope, b b b c c c α ∼ 2 \alpha\sim2 α ∼ 2 α ∼ 3 \alpha\sim3 α ∼ 3
Dark matter substructures
In this appendix we provide more details on the procedure that we followed to derive the constraints on f P B H f_{\mathrm{PBH}} f PBH § 12.1 .
Given a PBH with position relative to the pulsar given by r ( t ) = r 0 + v t \mathbf{r}(t)=\mathbf{r}_0+\mathbf{v}t r ( t ) = r 0 + v t r 0 \mathbf{r}_0 r 0 v \mathbf{v} v
h D ( t ) = G M v 2 d ^ ⋅ ( 1 + x D 2 b ^ D − sinh − 1 ( x ) v ^ ) (98) h_D(t)=\frac{GM}{v^2}\hat{\mathbf{d}}\cdot\left(\sqrt{1+x_D^2}\hat{\mathbf{b}}_D-\sinh^{-1}(x)\hat{\mathbf{v}}\right) \tag*{(98)} h D ( t ) = v 2 GM d ^ ⋅ ( 1 + x D 2 b ^ D − sinh − 1 ( x ) v ^ ) (98) h S ( t ) = 2 G M log ( 1 + x S 2 ) , (99) h_S(t)=2GM\log(1+x_S^2), \tag*{(99)} h S ( t ) = 2 GM log ( 1 + x S 2 ) , (99) where x D = ( t − t D , 0 ) / τ D x_D=(t-t_{D,0})/\tau_D x D = ( t − t D , 0 ) / τ D x S = ( t − t S , 0 ) / τ S x_S=(t-t_{S,0})/\tau_S x S = ( t − t S , 0 ) / τ S t t t [102 ] ; Lee et al. [210 ] ). For the static limit in which t ≫ T o b s t\gg T_{\mathrm{obs}} t ≫ T obs t 0 / τ t_0/\tau t 0 / τ O ( t 2 ) O(t^2) O ( t 2 ) ∝ t 3 \propto t^3 ∝ t 3
h ( t ) D ( S ) = A D ( S ) , s t a , 3 y r 2 t 3 , (100) h(t)_{\mathrm{D(S)}}=\frac{A_{\mathrm{D(S)},\mathrm{sta},3}}{\mathrm{yr}^2}t^3, \tag*{(100)} h ( t ) D ( S ) = yr 2 A D ( S ) , sta , 3 t 3 , (100) where A D ( S ) , s t a A_{\mathrm{D(S)},\mathrm{sta}} A D ( S ) , sta
A D , s t a t = γ r 2 G M 2 v 2 d ^ ⋅ [ t D , 0 τ D 4 1 ( 1 + t D , 0 2 / τ D 2 ) 5 / 2 b ^ D + 1 3 τ D 3 1 − 2 t D , 0 2 / τ D 2 ( 1 + t D , 0 2 / τ D 2 ) 5 / 2 v ^ ] (101) A_{\mathrm{D,stat}} = \gamma r^2 \frac{GM}{2v^2}\hat{\mathbf{d}}\cdot\left[\frac{t_{\mathrm{D},0}}{\tau_{\mathrm{D}}^4}\frac{1}{(1+t_{\mathrm{D},0}^2/\tau_{\mathrm{D}}^2)^{5/2}}\hat{\mathbf{b}}_{\mathrm{D}}+\frac{1}{3\tau_{\mathrm{D}}^3}\frac{1-2t_{\mathrm{D},0}^2/\tau_{\mathrm{D}}^2}{(1+t_{\mathrm{D},0}^2/\tau_{\mathrm{D}}^2)^{5/2}}\hat{\mathbf{v}}\right] \tag*{(101)} A D , stat = γ r 2 2 v 2 GM d ^ ⋅ [ τ D 4 t D , 0 ( 1 + t D , 0 2 / τ D 2 ) 5/2 1 b ^ D + 3 τ D 3 1 ( 1 + t D , 0 2 / τ D 2 ) 5/2 1 − 2 t D , 0 2 / τ D 2 v ^ ] (101) A S , s t a t = − γ r 2 4 G M 3 t S , 0 τ S 4 3 + t S , 0 2 / τ S 2 ( 1 + t S , 0 2 / τ S 2 ) 3 . (102) A_{\mathrm{S,stat}}=-\gamma r^2\frac{4GM}{3}\frac{t_{\mathrm{S},0}}{\tau_{\mathrm{S}}^4}\frac{3+t_{\mathrm{S},0}^2/\tau_{\mathrm{S}}^2}{(1+t_{\mathrm{S},0}^2/\tau_{\mathrm{S}}^2)^3}. \tag*{(102)} A S , stat = − γ r 2 3 4 GM τ S 4 t S , 0 ( 1 + t S , 0 2 / τ S 2 ) 3 3 + t S , 0 2 / τ S 2 . (102) For the Doppler case, in the dynamic limit when τ ≪ T o b s \tau\ll T_{\mathrm{obs}} τ ≪ T obs (98) where 1 + x D 2 ∝ ∣ t − t 0 ∣ \sqrt{1+x_{\mathrm{D}}^2}\propto\lvert t-t_0\rvert 1 + x D 2 ∝ ∣ t − t 0 ∣ x D x_{\mathrm{D}} x D
h D , d y n ( t ) = A D , d y n ( t − t 0 ) Θ ( t − t 0 ) , (103) h_{\mathrm{D,dyn}}(t)=A_{\mathrm{D,dyn}}(t-t_0)\Theta(t-t_0), \tag*{(103)} h D , dyn ( t ) = A D , dyn ( t − t 0 ) Θ ( t − t 0 ) , (103) where A D , d y n A_{\mathrm{D,dyn}} A D , dyn
A D , d y n = 2 G M τ v 2 d ^ ⋅ b ^ D . (104) A_{\mathrm{D,dyn}}=\frac{2GM}{\tau v^2}\hat{\mathbf{d}}\cdot\hat{\mathbf{b}}_{\mathrm{D}}. \tag*{(104)} A D , dyn = τ v 2 2 GM d ^ ⋅ b ^ D . (104) Given the expressions in Eqs. (98) and (C27) for the Doppler and Shapiro signals, we use the MC developed in Lee et al. [209 ] to derive the distributions p ( log 10 A I ∣ f P B H ) p(\log_{10} A_I\mid f_{\mathrm{PBH}}) p ( log 10 A I ∣ f PBH )
N P B H = f P B H ρ D M V M P B H , (105) N_{\mathrm{PBH}}=f_{\mathrm{PBH}}\frac{\rho_{\mathrm{DM}}V}{M_{\mathrm{PBH}}}, \tag*{(105)} N PBH = f PBH M PBH ρ DM V , (105) where the simulation volume, V V V R = v ˉ T o b s R=\bar{v}T_{\mathrm{obs}} R = v ˉ T obs v ˉ ≃ 340 k m / s \bar{v}\simeq340\ \mathrm{km/s} v ˉ ≃ 340 km/s T o b s , I T_{\mathrm{obs},I} T obs , I I I I
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 f(\mathbf{r})=1/V f ( r ) = 1/ V v 0 = 325 k m s − 1 v_0=325\ \mathrm{km\,s}^{-1} v 0 = 325 km s − 1 v e s c = 600 k m s − 1 v_{\mathrm{esc}}=600\ \mathrm{km\,s}^{-1} v esc = 600 km s − 1
The PBHs’ signals are then classified as dynamic if they satisfy the condition T o b s , I − τ > t 0 > τ T_{\mathrm{obs},I}-\tau>t_0>\tau T obs , I − τ > t 0 > τ
To evaluate A s t a t A_{\mathrm{stat}} A stat (98) and (C27), and we fit the resulting signal to a cubic polynomial in time to extract the t 3 t^3 t 3 A D , d y n A_{\mathrm{D,dyn}} A D , dyn (104) .
All the previous points are repeated for 2.5 × 10 3 2.5\times10^3 2.5 × 1 0 3 p ( log 10 A I ∣ f P B H ) p(\log_{10} A_I\mid f_{\mathrm{PBH}}) p ( log 10 A I ∣ f PBH ) p ( log 10 A I ∣ f P B H ) p(\log_{10} A_I\mid f_{\mathrm{PBH}}) p ( log 10 A I ∣ f PBH ) p ( log 10 A I ∣ δ t ) p(\log_{10} A_I\mid\delta t) p ( log 10 A I ∣ δ t )
p ( f P B H ∣ δ t ) = ∏ I = 1 N P ∫ p ( f P B H ∣ log 10 A I ) p ( log 10 A I ∣ δ t ) d log 10 A I . (106) p(f_{\mathrm{PBH}}\mid\delta t)=\prod_{I=1}^{N_{\mathrm{P}}}\int p(f_{\mathrm{PBH}}\mid\log_{10} A_I)p(\log_{10} A_I\mid\delta t)\,\mathrm{d}\log_{10} A_I. \tag*{(106)} p ( f PBH ∣ δ t ) = I = 1 ∏ N P ∫ p ( f PBH ∣ log 10 A I ) p ( log 10 A I ∣ δ t ) d log 10 A I . (106) Then, using Bayes’ theorem, we can rewrite
p ( f P B H ∣ log 10 A I ) = p ( log 10 A I ∣ f P B H ) p ( f P B H ) p ( log 10 A I ) . (107) p(f_{\mathrm{PBH}}\mid\log_{10} A_I)=\frac{p(\log_{10} A_I\mid f_{\mathrm{PBH}})\,p(f_{\mathrm{PBH}})}{p(\log_{10} A_I)}. \tag*{(107)} p ( f PBH ∣ log 10 A I ) = p ( log 10 A I ) p ( log 10 A I ∣ f PBH ) p ( f PBH ) . (107) Our priors on p ( f P B H ) p(f_{\mathrm{PBH}}) p ( f PBH ) p ( log 10 A I ) p(\log_{10} A_I) p ( log 10 A I ) (106) as
p ( f P B H ∣ δ t ) ∝ ∏ I = 1 N P ∫ p ( log 10 A I ∣ f P B H ) p ( log 10 A I ∣ δ t ) d log 10 A I , (108) p(f_{\mathrm{PBH}}\mid\delta t)\propto\prod_{I=1}^{N_{\mathrm{P}}}\int p(\log_{10} A_I\mid f_{\mathrm{PBH}})p(\log_{10} A_I\mid\delta t)\,\mathrm{d}\log_{10} A_I, \tag*{(108)} p ( f PBH ∣ δ t ) ∝ I = 1 ∏ N P ∫ p ( log 10 A I ∣ f PBH ) p ( log 10 A I ∣ δ t ) d log 10 A I , (108) where the ∝ \propto ∝ p ( f P B H ∣ δ t ) p(f_{\mathrm{PBH}}\mid\delta t) p ( f PBH ∣ δ t ) ∫ 0 ∞ p ( f P B H ∣ δ t ) d f P B H = 1 \int_0^\infty p(f_{\mathrm{PBH}}\mid\delta t)\,\mathrm{d}f_{\mathrm{PBH}}=1 ∫ 0 ∞ p ( f PBH ∣ δ t ) d f PBH = 1
In 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. [140 ]
( δ ν ν ) f i f t h = α ~ G M τ D 2 v d ^ ⋅ ∫ 1 ( 1 + x D 2 ) 3 / 2 ( 1 + b λ 1 + x D 2 ) e − ( b / λ ) 1 + x D 2 ( b ^ + x D v ^ ) d x D (109) \left(\frac{\delta\nu}{\nu}\right)_{\mathrm{fifth}} = \frac{\tilde{\alpha}GM}{\tau_D^2 v}\hat{d}\cdot\int\frac{1}{\left(1+x_D^2\right)^{3/2}}\left(1+\frac{b}{\lambda}\sqrt{1+x_D^2}\right)e^{-\left(b/\lambda\right)\sqrt{1+x_D^2}}\left(\hat{b}+x_D\hat{v}\right)\,dx_D \tag*{(109)} ( ν δ ν ) fifth = τ D 2 v α ~ GM d ^ ⋅ ∫ ( 1 + x D 2 ) 3/2 1 ( 1 + λ b 1 + x D 2 ) e − ( b / λ ) 1 + x D 2 ( b ^ + x D v ^ ) d x D (109) The integral in Eq. (109) 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 , t o t a l ( t ) = h D , f i f t h ( t ) + h D ( S ) ( t ) h_{D,\mathrm{total}}(t)=h_{D,\mathrm{fifth}}(t)+h_{D}^{(S)}(t) h D , total ( t ) = h D , fifth ( t ) + h D ( S ) ( t ) f P B H = 1 f_{\mathrm{PBH}}=1 f PBH = 1 h D , t o t a l ( t ) = A D , t o t a l y r 2 t 3 h_{D,\mathrm{total}}(t)=\frac{A_{D,\mathrm{total}}}{\mathrm{yr}^{2}}t^{3} h D , total ( t ) = yr 2 A D , total t 3 A D , t o t a l A_{D,\mathrm{total}} A D , total λ \lambda λ α ~ \tilde{\alpha} α ~ P ( log 10 A D , t o t a l ∣ λ , α ~ ) P(\log_{10}A_{D,\mathrm{total}}|\lambda,\tilde{\alpha}) P ( log 10 A D , total ∣ λ , α ~ ) α ~ \tilde{\alpha} α ~ λ \lambda λ P ( α ~ , λ ∣ δ t ) P(\tilde{\alpha},\lambda|\delta t) P ( α ~ , λ ∣ δ t ) (108)
p ( α ~ , λ ∣ δ t ) ∝ ∏ I = 1 N P ∫ p ( log 10 A I ∣ α ~ , λ ) p ( log 10 A I ∣ δ t ) d log 10 A I . (110) p(\tilde{\alpha},\lambda|\delta t)\propto\prod_{I=1}^{N_P}\int p(\log_{10}A_I|\tilde{\alpha},\lambda)p(\log_{10}A_I|\delta t)\,d\log_{10}A_I . \tag*{(110)} p ( α ~ , λ ∣ δ t ) ∝ I = 1 ∏ N P ∫ p ( log 10 A I ∣ α ~ , λ ) p ( log 10 A I ∣ δ t ) d log 10 A I . (110) 
