The NANOGrav 15 yr Data Set: Search for Signals from New Physics

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 2009). 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 GWs (IGWs; Grishchuk 1974; Starobinsky 1979; Rubakov et al. 1982; Fabbri & Pollock 1983; Abbott & Wise 1984) 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 (2021), Kuroyanagi et al. (2021), and Benetti et al. (2022).

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. For the simplest type of inflation—standard single-field slow-roll inflation—the h²Ω_GW spectrum is red-tilted at CMB scales, with the tensor spectral index n_t being given by the so-called consistency relation, n_t = − r/8 < 0. Meanwhile, r is bounded from above by current CMB observations, r ≤ 0.036 at 95% C.L. (Ade et al. 2021). A vanishing tensor spectral index, n_t ≈ 0, would imply an upper bound on the GW energy density spectrum of Ω_GW h² ∼ 10⁻¹⁶ at PTA frequencies, rendering any detection of an IGW signal in PTA observations hopeless. This conclusion, however, only applies to the standard case of single-field slow-roll inflation. Nonminimal scenarios may have significantly better detection prospects.

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 and tensor spectral index n_t at the CMB pivot scale, f_CMB = 0.05 Mpc⁻¹/(2π a₀) ≃ 7.73 × 10⁻¹⁷ Hz, which quantify the efficiency and scale dependence of GW production during inflation, and the reheating temperature T_rh and the number of e-folds during reheating N_rh, which describe the reheating process after inflation. Here the factor a₀ in the expression for f_CMB denotes the present value of the cosmological scale factor a(t) in the Robertson–Walker metric; in our convention, a₀ = 1.

We do not impose the standard consistency relation between r and n_t ; instead, we allow both parameters to vary independently across large prior ranges, ${\mathrm{log}}_{10}r\in [-40,0]$ and n_t ∈ [0, 6]. We note that blue values of the tensor spectral index can be generated, e.g., from axion–vector dynamics during inflation (Anber & Sorbo 2012; Cook & Sorbo 2012; Namba et al. 2016; Dimastrogiovanni et al. 2017; Caldwell & Devulder 2018) or in other nonminimal inflation models (see Piao & Zhang 2004; Satoh & Soda 2008; Kobayashi et al. 2010; Endlich et al. 2013; Fujita et al. 2019 for an incomplete list).

Similarly, we allow for more flexibility for T_rh and N_rh than in the standard treatment of single-field slow-roll inflation. To illustrate this point, note that the number of e-folds during reheating, N_rh, can be written as (Liddle & Lyth 2009)

$$\begin{eqnarray}&&{N}_{\mathrm{rh}}=\displaystyle \frac{1}{3\left(1+{w}_{\mathrm{rh}}\right)}\,\mathrm{ln}\left(\displaystyle \frac{3{H}_{\mathrm{end}}^{2}{M}_{\mathrm{Pl}}^{2}}{{\pi }^{2}/30\,{g}_{* }^{\mathrm{rh}}\,{T}_{\mathrm{rh}}^{4}}\right),\end{eqnarray} \tag{ 14 }$$

where H_end is the Hubble rate at the end of inflation, w_rh is the equation-of-state parameter during reheating, M_Pl ≃ 2.44 × 10¹⁸ GeV is the reduced Planck mass, and ${g}_{* }^{\mathrm{rh}}$ is defined below. We assume for definiteness that reheating is dominated by the coherent oscillations of the inflaton field, such that the equation of state is equivalent to the one of pressureless dust (i.e., matter), w_rh = 0. In typical models of single-field slow-roll inflation, one can often approximate H_end by the Hubble rate at the time of CMB horizon exit during inflation, such that

$$\begin{eqnarray}&&{N}_{\mathrm{rh}}^{\mathrm{naive}}\approx \displaystyle \frac{1}{3}\mathrm{ln}\left(\displaystyle \frac{3{H}_{\mathrm{naive}}^{2}{M}_{\mathrm{Pl}}^{2}}{{\pi }^{2}/30\,{g}_{* }^{\mathrm{rh}}\,{T}_{\mathrm{rh}}^{4}}\right),\end{eqnarray} \tag{ 15 }$$

where H_naive is fixed by the tensor-to-scalar ratio r and the amplitude of the primordial scalar power spectrum, A_s ≃ 2.10 × 10⁻⁹ (Aghanim et al. 2020),

$$\begin{eqnarray}&&{H}_{\mathrm{end}}\approx {H}_{\mathrm{naive}}={\left(\displaystyle \frac{{\pi }^{2}}{2}{{rA}}_{s}\right)}^{1/2}{M}_{\mathrm{Pl}}.\end{eqnarray} \tag{ 16 }$$

However, we already assume nonminimal dynamics in order to motivate a strongly blue-tilted h²Ω_GW spectrum, so there is no reason why we should make use of this approximation. In our analysis, we therefore treat H_end and correspondingly N_rh as independent parameters and do not fix them in terms of r and T_rh as in Equations (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 (Abbott et al. 2021a) and on the amount of dark radiation, i.e., the effective number of neutrino species, N_eff, inferred from big bang nucleosynthesis (BBN) and the CMB (Pisanti et al. 2021; Yeh et al. 2021). As for the latter constraint, we specifically work with ΔN_eff = ρ_DR/ρ_ν , where ρ_DR is the energy density of dark radiation (i.e., the integrated GW energy density in the context of the igw model) and ρ_ν denotes the energy density of a single neutrino species. ΔN_eff characterizes the excess energy in radiation beyond the SM expectation (i.e., dark radiation) after neutrino decoupling and e⁺ e⁻ annihilation, ${N}_{\mathrm{eff}}={N}_{\mathrm{eff}}^{\mathrm{SM}}+{\rm{\Delta }}{N}_{\mathrm{eff}}$, where ${N}_{\mathrm{eff}}^{\mathrm{SM}}\,\simeq 3.0440$ (Bennett et al. 2021).

Under the assumptions outlined above, we are able to model the IGW spectrum at PTA frequencies as

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{\mathrm{GW}}^{\inf }\left(f\right)=\displaystyle \frac{{{\rm{\Omega }}}_{{\rm{r}}}}{24}\left(\displaystyle \frac{{g}_{* }\left(f\right)}{{g}_{* }^{0}}\right){\left(\displaystyle \frac{{g}_{* ,s}^{0}}{{g}_{* ,s}\left(f\right)}\right)}^{4/3}{{ \mathcal P }}_{t}\left(f\right){ \mathcal T }\left(f\right).\end{eqnarray} \tag{ 17 }$$

Here ${{\rm{\Omega }}}_{{\rm{r}}}/{g}_{* }^{0}\simeq 2.72\times {10}^{-5}$ is the current radiation energy density per relativistic degree of freedom, in units of the critical (closure) density, ${g}_{* ,s}^{0}\simeq 3.93$ counts the effective number of relativistic degrees of freedom contributing to the radiation entropy today, and g_*(f) and g_*,s (f) denote the effective numbers of relativistic degrees of freedom in the early universe when GWs with comoving wavenumber k = 2π a₀ f reentered the Hubble horizon after inflation. In order to evaluate g_*(f) and g_*,s (f), we use the numbers of relativistic degrees of freedom as functions of temperature tabulated in Saikawa & Shirai (2020), g_*(T) and g_*,s (T), in combination with the standard temperature–frequency relation in ΛCDM (i.e., the cosmological Lambda cold dark matter SM) that follows from the condition k = a(T)H(T) at the time of horizon reentry. In the remainder of this paper, whenever we need g_* or g_*,s in a different part of our analysis, we will use the same functions g_*(f), g_*,s (f), g_*(T), and g_*,s (T).

The inflationary dynamics give rise to the primordial tensor power spectrum ${{ \mathcal P }}_{t}$, while the transfer function ${ \mathcal T }$ accounts for the redshifting behavior of GWs after horizon reentry. In our analysis, we assume a constant tensor spectral tilt (i.e., zero running of n_t ) from CMB to PTA frequencies, such that

$$\begin{eqnarray}&&{{ \mathcal P }}_{t}\left(f\right)=r\,{A}_{s}{\left(\displaystyle \frac{f}{{f}_{\mathrm{CMB}}}\right)}^{{n}_{t}}.\end{eqnarray} \tag{ 18 }$$

Meanwhile, the only relevant contribution to ${ \mathcal T }$ in the PTA band corresponds to the transfer function that connects the radiation-dominated era to reheating,

$$\begin{eqnarray}&&{ \mathcal T }\left(f\right)\approx \displaystyle \frac{{\rm{\Theta }}\left({f}_{\mathrm{end}}-f\right)}{1-0.22{\left(f/{f}_{\mathrm{rh}}\right)}^{1.5}+0.65{\left(f/{f}_{\mathrm{rh}}\right)}^{2}}.\end{eqnarray} \tag{ 19 }$$

Here the fit function in the denominator of this expression is taken from Kuroyanagi et al. (2015, 2021) and describes the spectral turnover, ${f}^{{n}_{t}}\to {f}^{{n}_{t}-2}$, at frequencies around f ∼ f_rh, which marks the end of the reheating period,

$$\begin{eqnarray}&&{f}_{\mathrm{rh}}=\displaystyle \frac{1}{2\pi }\ {\left(\displaystyle \frac{{g}_{* ,s}^{0}}{{g}_{* ,s}^{\mathrm{rh}}}\right)}^{1/3}{\left(\displaystyle \frac{{\pi }^{2}{g}_{* }^{\mathrm{rh}}}{90}\right)}^{1/2}\displaystyle \frac{{T}_{\mathrm{rh}}{T}_{0}}{{M}_{\mathrm{Pl}}},\end{eqnarray} \tag{ 20 }$$

with ${g}_{* }^{\mathrm{rh}}={g}_{* }({T}_{\mathrm{rh}})$ and ${g}_{* ,s}^{\mathrm{rh}}={g}_{* ,s}({T}_{\mathrm{rh}})$ and the present-day CMB temperature T₀ ≃ 2.73 K (Fixsen 2009). The Heaviside theta function in Equation (19) denotes the endpoint of the IGW spectrum at f = f_end, which marks the end of inflation and hence the onset of reheating,

$$\begin{eqnarray}&&{f}_{\mathrm{end}}=\displaystyle \frac{1}{2\pi }\ {\left(\displaystyle \frac{{g}_{* ,s}^{0}}{{g}_{* ,s}^{\mathrm{rh}}}\right)}^{1/3}{\left(\displaystyle \frac{{\pi }^{2}{g}_{* }^{\mathrm{rh}}}{90}\right)}^{1/3}\displaystyle \frac{{T}_{\mathrm{rh}}^{1/3}{H}_{\mathrm{end}}^{1/3}\,{T}_{0}}{{M}_{\mathrm{Pl}}^{2/3}}.\end{eqnarray} \tag{ 21 }$$

For fixed T_rh, the frequency f_end is solely controlled by H_end, which follows from our choice of N_rh according to Equation (14) (recall that we set w_rh = 0). In our MCMC analysis, we do not sample over f_end, since its precise value does not affect the shape of the IGW spectrum in the PTA band and hence the quality of our fit. Instead, we constrain its maximally allowed value (i.e., N_rh) after identifying the viable region in the r–n_t –T_rh parameter space, in order to make sure that the IGW spectrum does not violate the N_eff and LVK bounds. The frequency f_rh, on the other hand, can easily fall into the PTA band: from Equation (20), we have ${f}_{\mathrm{rh}}\sim 30\,{\rm{nHz}}\left({T}_{\mathrm{rh}}/1\,{\rm{GeV}}\right)$. Therefore, we sample the reheating temperature in a symmetric interval around T_rh = 1 GeV extending down to temperatures relevant for BBN, T_BBN ∼ 1 MeV. That is, we work with the log-uniform prior ${\mathrm{log}}_{10}\left({T}_{\mathrm{rh}}/1\,{\rm{GeV}}\right)\in \left[-3,+3\right]$.

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, ${ \mathcal B }=8.8\pm 0.3$ (mean value and one standard deviation), and simply follows from the larger parametric freedom of the igw model. Both the igw and smbhb models basically correspond to a power-law approximation of the GW spectrum. However, in the case of the igw model, the parameters controlling this power law are drawn from prior distributions that allow for a larger amplitude and a steeper slope of the spectrum, which improves the quality of the fit. Meanwhile, the combined GW spectrum from inflation with an additional SMBHB signal on top compared to the smbhb model alone results in a Bayes factor of ${ \mathcal B }=7.9\pm 0.3$. We thus observe a slight decrease in the Bayes factor, which accounts for the fact that adding the SMBHB signal on top of the IGW signal does not improve the quality of fit but merely increases the prior volume compared to the igw model.

The reconstructed posterior distributions for the parameters of the igw model and its igw+smbhb extension are shown in Figure 5. ⁸¹ For both models, we find a strong covariance between the spectral index n_t and the tensor-to-scalar ratio r, which is approximately fit by

$$\begin{eqnarray}&&{n}_{t}=-0.14{\mathrm{log}}_{10}r+0.58,\end{eqnarray} \tag{ 22 }$$

and which can be explained as follows: the igw interpretation of the PTA signal requires the primordial tensor power spectrum ${{ \mathcal P }}_{t}$ to take values of ${ \mathcal O }\left({10}^{-4}\right)$ at nanohertz frequencies. This requirement fixes the parameter combination $r{\left({f}_{\mathrm{PTA}}/{f}_{\mathrm{CMB}}\right)}^{{n}_{t}}$ in Equation (18) and thus allows us to estimate the coefficients in Equation (22) as $1/{\mathrm{log}}_{10}\left({f}_{\mathrm{PTA}}/{f}_{\mathrm{CMB}}\right)\approx 0.14$ and ${\mathrm{log}}_{10}\left({{ \mathcal P }}_{t}/{A}_{s}\right)/{\mathrm{log}}_{10}\left({f}_{\mathrm{PTA}}/{f}_{\mathrm{CMB}}\right)\approx 0.58$, respectively, where we used f_PTA = 1 nHz and ${{ \mathcal P }}_{t}=0.3\times {10}^{-4}$.

Zoom In Zoom Out Reset image size

In addition to the strong covariance between n_t and r, we note that the posterior probabilities of both parameters exhibit a bimodal distribution for both igw and igw+smbhb. In the 2D distributions of the parameter pairs $\left({T}_{\mathrm{rh}},{n}_{t}\right)$ and $\left({T}_{\mathrm{rh}},r\right)$, this bimodality is accompanied by an approximate reflection symmetry with respect to the points $\left({\mathrm{log}}_{10}{T}_{\mathrm{rh}}/{\rm{GeV}},{n}_{t}\right)\sim \left(-0.5,2.75\right)$ and $\left({\mathrm{log}}_{10}{T}_{\mathrm{rh}}/{\rm{GeV}},{\mathrm{log}}_{10}r\right)\sim \left(-0.5,-15\right)$, respectively. These features of the corner plot in Figure 5 indicate that the igw model can operate in two regimes: for T_rh ≫ 1 GeV, the reference frequency f_rh is larger than the frequencies in the PTA band, and the GW spectrum seen by NANOGrav is composed of tensor modes that reentered the horizon during the radiation-dominated era. For T_rh ≪ 1 GeV, on the other hand, f_rh can be pushed below PTA frequencies, and the GW spectrum in the PTA band is composed of tensor modes that reentered the horizon during reheating after inflation. In the first case, the tilt of the spectrum is directly given by n_t ; in the second case, it corresponds to n_t − 2. Clearly, the mirror symmetry in the 2D distributions of $\left({T}_{\mathrm{rh}},{n}_{t}\right)$ and $\left({T}_{\mathrm{rh}},r\right)$ is not exact. At the level of the GW spectrum, it is explicitly broken by the frequency dependence of g_* and g_*,s , as well as by the nontrivial shape of the transfer function in Equation (19).

At small or large values of the reheating temperature, the posterior distributions develop approximately flat directions along the T_rh axis at $\left({n}_{t},{\mathrm{log}}_{10}r\right)\sim \left(2,-10\right)$ in the large-T_rh regime and at $\left({n}_{t},{\mathrm{log}}_{10}r\right)\sim \left(4,-20\right)$ in the small-T_rh regime. This behavior is broadly consistent with the linear fit in Equation (22), the reflection symmetry discussed above, and the fact that, past a certain point, raising or lowering T_rh no longer influences the shape of the GW signal in the PTA band. A tensor index of n_t = 2 at large T_rh, moreover, corresponds to an index γ = 3 in the timing residual PSD, which is among the best-fitting values—see the 2D posterior for A and γ in Figure 1. The same conclusion holds at low T_rh, where γ = 3 is mapped onto n_t − 2 = 2.

Finally, we derive the N_eff and LVK bounds on the igw parameter space. For the N_eff bound, the GW spectrum, integrated from BBN scales to the cutoff frequency f_end, must not exceed a certain upper limit that is set by the allowed amount of extra relativistic degrees of freedom at the time of BBN and recombination,

$$\begin{eqnarray}&&{\int }_{{f}_{\mathrm{BBN}}}^{{f}_{\mathrm{end}}}\displaystyle \frac{{df}}{f}\ {h}^{2}{{\rm{\Omega }}}_{\mathrm{GW}}^{\inf }\left(f\right)\lesssim 5.6\times {10}^{-6}\,{\rm{\Delta }}{N}_{\mathrm{eff}}^{\max }.\end{eqnarray} \tag{ 23 }$$

Here f_BBN ∼ 10⁻¹² Hz refers to tensor modes that reentered the Hubble horizon around the onset of BBN at T ∼ 10⁻⁴ GeV (Caprini & Figueroa 2018), and ${\rm{\Delta }}{N}_{\mathrm{eff}}^{\max }\sim {\rm{few}}\times 0.1$ denotes the upper limit on the amount of dark radiation. For definiteness, we set f_BBN = 10⁻¹² Hz and ${\rm{\Delta }}{N}_{\mathrm{eff}}^{\max }=0.5$ in our analysis; the precise ${\rm{\Delta }}{N}_{\mathrm{eff}}^{\max }$ value at 95% confidence level varies across different combinations of data sets (Pisanti et al. 2021; Yeh et al. 2021). For given values of the parameters T_rh, r, and n_t , we can then ask whether there is a cutoff frequency ${f}_{\mathrm{end}}^{\max }$ that leads to the saturation of the inequality in Equation (23). If this is the case, ${f}_{\mathrm{end}}^{\max }$ yields an upper bound on the allowed number of e-folds during reheating, ${N}_{\mathrm{rh}}^{\max }$, according to Equations (14) and (21). A second constraint on N_rh follows from the LVK bound on the amplitude of the GWB (Abbott et al. 2021a),

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{\mathrm{GW}}\leqslant 1.7\times {10}^{-8}\quad {\rm{at}}\quad {f}_{\mathrm{lvk}}=25\,{\rm{Hz}},\end{eqnarray} \tag{ 24 }$$

assuming a flat GW spectrum. For a power-law spectrum, ${{\rm{\Omega }}}_{\mathrm{GW}}={{\rm{\Omega }}}_{\mathrm{GW},\alpha }{\left(f/{f}_{\mathrm{lvk}}\right)}^{\alpha }$ with α = n_t − 2 at f ≫ f_rh, this bound can be generalized following (Kuroyanagi et al. 2015, 2021)

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{\mathrm{GW}}\leqslant 1.7\times {10}^{-8}{\left(\frac{5-2\alpha }{5}\right)}^{1/2}{\left(\frac{20\,\mathrm{Hz}}{{f}_{\mathrm{lvk}}}\right)}^{-\alpha },\end{eqnarray} \tag{ 25 }$$

which approximately holds if α ≪ 5/2. For given values of T_rh, r, and n_t , we can then evaluate ${{\rm{\Omega }}}_{\mathrm{GW}}^{\inf }$ at f_lvk and check whether it exceeds the LVK bound. If so, we place an upper bound on N_rh by demanding that ${f}_{\mathrm{end}}^{\max }=20\,{\rm{Hz}}$ (the lower end of the LVK frequency band).

The outcome of our analysis is shown Figure 22 in Appendix C.1. In both plots of this figure, the dependence on the tensor-to-scalar ratio r is eliminated by means of the linear relation in Equation (22). In the left panel of Figure 22, we present contour lines of ${N}_{\mathrm{rh}}^{\max }$ derived from the N_eff bound and the LVK bound, respectively. In a realization of the igw model with a given duration of reheating, these contour lines can be thought of as bounds on the T_rh–n_t parameter plane—for fixed N_rh, the regions with ${N}_{\mathrm{rh}}^{\max }\lt {N}_{\mathrm{rh}}$ are excluded. We find that the N_eff bound rules out most values of the spectral index n_t in the large-T_rh regime. At the same time, large regions of parameter space remain viable as long as ${N}_{\mathrm{rh}}\lt {N}_{\mathrm{rh}}^{\max }$. In fact, away from the region where ${N}_{\mathrm{rh}}^{\max }$ turns negative, our upper bound is typically rather large, ${N}_{\mathrm{rh}}^{\max }\sim { \mathcal O }\left(10\cdots 100\right)$, and hence easy to satisfy in realistic models of reheating, where ${N}_{\mathrm{rh}}\sim { \mathcal O }\left(1\cdots 10\right)$. In the right panel of Figure 22, we compare our result to the naive expectation ${N}_{\mathrm{rh}}^{\mathrm{naive}}$ in single-field slow-roll inflation with a nearly constant Hubble rate (see Equation (15)). Across large parts of parameter space, we find that ${N}_{\mathrm{rh}}^{\mathrm{naive}}$ assumes unrealistically large values, ${N}_{\mathrm{rh}}^{\mathrm{naive}}\gg 10$.

In Figure 5, we highlight the constraints on T_rh and n_t (as well as T_rh and r) that we deduce from the N_eff bound assuming N_rh values of N_rh = 0, 5, and 10. Here the constraints for N_rh = 0 correspond to the assumption of instantaneous reheating after inflation and hence represent the most conservative bound on the T_rh–n_t parameter plane. A longer duration of reheating results in tighter constraints on T_rh and n_t , as illustrated by the contours for N_rh = 5 and 10. For an even larger number of e-folds during reheating, see Figure 22 in Appendix C.1.

In view of Figures 5 and 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, a tiny tensor-to-scalar ratio, r ∼ 10^−(23⋯16), a low reheating temperature, T_rh ∼ 10^−(3⋯0) GeV, and a moderate number of e-folds during reheating, N_rh ≲ 20. It remains to be seen whether it is possible to identify explicit microscopic models that realize inflation in this parametric regime.

The amplitude of the primordial scalar power spectrum is well measured by CMB observations, A_s ≃ 2.10 × 10⁻⁹ at the CMB pivot scale k_CMB = 0.05 Mpc⁻¹ (Aghanim et al. 2020). If we naively extrapolate this value down to smaller scales, assuming a fixed and slightly red-tilted h²Ω_GW spectrum with index n_s ∼ 0.96, we are led to conclude that there must be increasingly less power in scalar perturbations on shorter scales. This conclusion can, however, be easily avoided in models that deviate from the standard picture of single-field slow-roll inflation giving rise to a nearly scale-invariant spectrum of scalar perturbations. A prominent example, among many other mechanisms, consists in a stage of inflation close to an inflection point in the scalar potential, which readily amplifies the scalar perturbations leaving the horizon (see, e.g., Garcia-Bellido & Ruiz Morales 2017; Ballesteros & Taoso 2018; Ezquiaga et al. 2018). An enhanced scalar power spectrum at small scales is, therefore, a viable possibility. Moreover, it promises a rich phenomenology with regard to the production of GWs and potentially the origin of primordial black holes (PBHs; Carr et al. 2016; Garcia-Bellido et al. 2016; Inomata et al. 2017a; Inomata & Nakama 2019; Wang et al. 2019; Escrivà et al. 2022b). The possibility of having PBH formation in models of single-field inflation is the subject of ongoing debate (Kristiano & Yokoyama 2022, 2023; Choudhury et al. 2023a, 2023b, 2023c; Firouzjahi & Riotto 2023; Riotto 2023a, 2023b). 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 SIGWs (Matarrese et al. 1993, 1994, 1998; Mollerach et al. 2004; Ananda et al. 2007; Baumann et al. 2007; see Domènech 2021 for a review and more details on the history of SIGWs). At the same time, large overdensities in the tail of the distribution of scalar perturbations can collapse into PBHs upon horizon reentry. This PBH production mechanism thus results in a PBH population whose properties are strongly correlated with the spectral shape of the SIGW signal, as both phenomena are controlled by the scalar spectrum generated during inflation (Yuan & Huang 2021a). For earlier works on the PBH/SIGW interpretation of the signal in recent PTA data sets, see Vaskonen & Veermäe (2021), De Luca et al. (2021), and Kohri & Terada (2021). For earlier Bayesian searches for an SIGW signal in PTA data sets, see Chen et al. (2020), Bian et al. (2021), Zhao & Wang (2022), Yi & Fei (2023), and Dandoy et al. (2023), and for a search in LVK data, see Romero-Rodriguez et al. (2022).

In our analysis, we consider SIGWs in the PTA band and model the associated GW spectrum as follows:

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{\mathrm{GW}}^{\mathrm{ind}}\left(f\right)={{\rm{\Omega }}}_{r}\left(\frac{{g}_{* }\left(f\right)}{{g}_{* }^{0}}\right){\left(\frac{{g}_{* ,s}^{0}}{{g}_{* ,s}\left(f\right)}\right)}^{4/3}{\bar{{\rm{\Omega }}}}_{\mathrm{GW}}^{\mathrm{ind}}\left(f\right),\end{eqnarray} \tag{ 26 }$$

where the first three factors account for the cosmological redshift as in Equation (17) and the last factor denotes the GW spectrum at the time of production, which we assume to be during the radiation-dominated era,

$$\begin{eqnarray}&&{\bar{{\rm{\Omega }}}}_{\mathrm{GW}}^{\mathrm{ind}}\left(f\right)={\int }_{0}^{\infty }{dv}{\int }_{\left|1-v\right|}^{1+v}{du}\ { \mathcal K }\left(u,v\right){{ \mathcal P }}_{{R}}\left({uk}\right){{ \mathcal P }}_{{R}}\left({vk}\right).\end{eqnarray} \tag{ 27 }$$

The present-day GW frequency is related to the comoving wavenumber by $f=k/\left(2\pi {a}_{0}\right)$, ${{ \mathcal P }}_{{R}}$ denotes the primordial spectrum of the comoving curvature perturbation ${ \mathcal R }$, and the integration kernel ${ \mathcal K }$ is given by (Espinosa et al. 2018; Kohri & Terada 2018; Pi & Sasaki 2020; Gong 2022)

$$\begin{eqnarray}\begin{array}{rcl}{ \mathcal K }\left(u,v\right) & = & \displaystyle \frac{3{\left(4{v}^{2}-{\left(1+{v}^{2}-{u}^{2}\right)}^{2}\right)}^{2}{\left({u}^{2}+{v}^{2}-3\right)}^{4}}{1024\,{u}^{8}{v}^{8}}\\ & & \times \,\left[{\left(\mathrm{ln}\left|\displaystyle \frac{3-{\left(u+v\right)}^{2}}{3-{\left(u-v\right)}^{2}}\right|-\displaystyle \frac{4{uv}}{{u}^{2}+{v}^{2}-3}\right)}^{2}+{\pi }^{2}{\rm{\Theta }}(u+v-\sqrt{3})\right].\end{array}\end{eqnarray} \tag{ 28 }$$

The expression in Equation (27) illustrates the dependence of the SIGW signal on the scalar input spectrum; in particular, it shows that ${{\rm{\Omega }}}_{\mathrm{GW}}^{\mathrm{ind}}$ scales as ${{\rm{\Omega }}}_{\mathrm{GW}}^{\mathrm{ind}}\propto {{ \mathcal P }}_{{R}}^{2}$. We stress that the expression in Equation (27) assumes Gaussian perturbations, whose statistics are fully described by the power spectrum ${{ \mathcal P }}_{{R}}$. We do not consider any non-Gaussian contributions to the primordial scalar power spectrum in our analysis. The impact of non-Gaussianities on the SIGW signal was studied in Nakama et al. (2017), Cai et al. (2019), Unal (2019), Atal et al. (2019), Yuan & Huang (2021b), Atal & Domènech (2021), Adshead et al. (2021), and Ferrante et al. (2023).

To remain as model-independent as possible, we refrain from choosing a particular inflation model capable of generating an enhanced spectrum ${{ \mathcal P }}_{{R}}$. Instead, we ignore the microphysics of inflation and work with three characteristic templates for ${{ \mathcal P }}_{{R}}$ that reflect the range of possibilities that one may expect in realistic models. Specifically, we consider the following templates:

sigw-delta: Sharp spectral feature in ${{ \mathcal P }}_{{R}}$ modeled by a Dirac delta function in logarithmic k space,

$$\begin{eqnarray}&&{{ \mathcal P }}_{{R}}\left(k\right)=A\,\delta \left(\mathrm{ln}k-\mathrm{ln}{k}_{* }\right).\end{eqnarray} \tag{ 29 }$$

sigw-Gauss: Broad spectral feature in ${{ \mathcal P }}_{{R}}$ modeled by a Gaussian peak in logarithmic k space,

$$\begin{eqnarray}&&{{ \mathcal P }}_{{R}}\left(k\right)=\displaystyle \frac{A}{\sqrt{2\pi }\,{\rm{\Delta }}}\,\exp \left[-\displaystyle \frac{1}{2}{\left(\displaystyle \frac{\mathrm{ln}k-\mathrm{ln}{k}_{* }}{{\rm{\Delta }}}\right)}^{2}\right].\end{eqnarray} \tag{ 30 }$$

sigw-box: Flat and continuous spectral feature in ${{ \mathcal P }}_{{R}}$ modeled by a box function in logarithmic k space,

$$\begin{eqnarray}&&{{ \mathcal P }}_{{R}}\left(k\right)=A\,{\rm{\Theta }}\left(\mathrm{ln}{k}_{\max }-\mathrm{ln}k\right){\rm{\Theta }}\left(\mathrm{ln}k-\mathrm{ln}{k}_{\min }\right).\end{eqnarray} \tag{ 31 }$$

Note that the Gaussian power spectrum in logarithmic k space that we assume in the sigw-Gauss model corresponds to a lognormal power spectrum in linear k space. As evident from the above expressions, sigw-delta represents a two-parameter model, while sigw-Gauss and sigw-box are three-parameter models. Our prior choices for the respective parameters are listed in Table 3 in Appendix A, where we use again $f=k/\left(2\pi {a}_{0}\right)$ to convert from wavenumber to frequency. For a given set of parameter values, we are then able to use the scalar power spectrum in Equation (29), Equation (30), or Equation (31) to evaluate the integrals in Equation (27) 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 Wang et al. (2019) and Yuan & Huang (2021a).

We now turn to the outcome of our Bayesian fit analysis. Compared to the igw model discussed in Section 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 ${ \mathcal B }=44\pm 2$, ${ \mathcal B }=57\pm 3$, and ${ \mathcal B }=21\pm 1$ for sigw-delta, sigw-Gauss, and sigw-box, respectively, and ${ \mathcal B }=38\pm 2$, ${ \mathcal B }=59\pm 3$, and ${ \mathcal B }=23\pm 1$ for sigw-delta+smbhb, sigw-Gauss+smbhb, and sigw-box+smbhb, respectively (see 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 uncertaintity of our bootstrapping analysis.

The reconstructed posterior distributions for all SIGW models under consideration are shown in Figures 6 and 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. We can quantify this statement in terms of the lower limits of the 95% Bayesian credible intervals for A that we obtain from integrating the corresponding 1D marginalized posteriors over the regions of highest posterior density: for sigw-delta, sigw-Gauss, and sigw-box, we respectively find ${\mathrm{log}}_{10}A\gtrsim -\,1.55$, ${\mathrm{log}}_{10}A\gtrsim -\,1.43$, and ${\mathrm{log}}_{10}A\gtrsim -\,1.90$. At the same time, the enhanced power in scalar perturbations must be localized at the right scales, so that the resulting SIGW signal falls into the PTA band. This requirement leads to bounds on the frequencies ${f}_{\min }$, ${f}_{\max }$, and f_* that can be read off from Figures 6 and 7 and that are summarized in Table 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A notable feature in this context is that the posterior distributions for ${f}_{\min }$, ${f}_{\max }$, and f_* all extend to large frequencies, much beyond the upper limit of the NANOGrav band. The reason for this is that the NG15 data are best fit by the low-frequency tail of the SIGW spectrum (see Figure 19). Beyond the NANOGrav band, the SIGW spectrum keeps increasing until it reaches its maximal value at f ≫ 1 nHz. This observation also explains the flat directions in the 2D posterior distribution for A and ${f}_{\min }$ in Figure 6 and the 2D posterior distributions for A and f_* in Figure 7. A simultaneous increase in A and ${f}_{\min }$ or f_* along these flat directions moves the peak in the GW spectrum to higher frequencies, but it preserves the shape of the low-frequency tail in the PTA band and hence does not affect the quality of the fit.

We also observe that the 2D posterior distributions for A and f_* for sigw-delta and sigw-Gauss are roughly similar to each other, with the distribution for the sigw-Gauss model being slightly broader. This result is consistent with the fact that sigw-delta and sigw-Gauss are nested models; sigw-delta can be obtained from sigw-Gauss in the limit Δ → 0. The slightly broader posterior distribution for A and f_* in the right panel of 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 ratio in Equation (10) falls below 1/10. In these regions, which are shaded in gray in 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 an SIGW signal. For sigw-box, we are not able to place a K-ratio bound on the 2D parameter space spanned by A and ${f}_{\min }$, because we additionally marginalize over ${f}_{\max }$. That is, for any pair of values of A and ${f}_{\min }$, the SIGW signal can be made arbitrarily small if we choose ${f}_{\max }$ close enough to ${f}_{\min }$.

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,

$$\begin{eqnarray}&&{f}_{\mathrm{PBH}}\leqslant 1,\end{eqnarray} \tag{ 32 }$$

where f_PBH = Ω_PBH/Ω_DM denotes the PBH DM fraction integrated over the entire PBH mass spectrum. We provide more details on how we compute f_PBH in Appendix C.2; here we simply report our final results in terms of the teal contour lines labeled f_PBH = 1 in Figures 6 and 7. For sigw-box, the PBH bound in the A–${f}_{\max }$ plane is computed for ${f}_{\min }={10}^{-11}\,{\rm{Hz}}$ (the lower end of our prior range), and the PBH bound in the A–${f}_{\min }$ plane is computed for ${f}_{\max }={10}^{-5}\,{\rm{Hz}}$ (the upper end of our prior range). For sigw-Gauss, we show the PBH bound in the A–f_* plane for Δ = 1 (solid contour line) and Δ = 2 (dashed contour line).

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_PBH is very sensitive to various assumptions and numerical steps in the analysis. Slight changes in the computational strategy may result in very different results for f_PBH, which is why the results reported in Figures 6 and 7 need to be treated with caution. In view of the large conceptional uncertainty in the computation of f_PBH, one needs to be careful not to draw any premature conclusions. At the same time, the PBH bounds in Figures 6 and 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_PBH = 1 inside the 68% credible regions. That is, along the f_PBH = 1 contour lines and inside the 68% credible regions, we find scenarios where SIGWs manage to provide a good fit to the NG15 data and PBHs account for the entire DM in our universe.

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 2022, 2023). Kristiano & Yokoyama (2022, 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 must be small, ${\mathrm{log}}_{10}A\ll -\,2$ otherwise, the loop corrections to the scalar power spectrum will exceed the measured amplitude at CMB scales, A_s ≃ 2.10 × 10⁻⁹. At present, this claim is subject to an ongoing debate; it is notably challenged by Riotto (2023a, 2023b) and Firouzjahi & Riotto (2023). However, if it should prove to be valid, the requirement of ${\mathrm{log}}_{10}A\ll -2$ would clash with the lower bounds on A listed above. In this case, one would then have to either give up on the SIGW interpretation of the NANOGrav signal or seek to construct inflation models that evade the arguments of Kristiano & Yokoyama (2022, 2023) and still lead to a large SIGW signal.

In the cosmological context, first-order phase transitions take place when a potential barrier separates the true and false minima of a scalar field potential. ⁸² The phase transition is triggered by quantum or thermal fluctuations that cause the scalar field to tunnel through or fluctuate over the barrier, which results in the nucleation of bubbles within which the scalar field is in the true vacuum configuration. If sufficiently large, these bubbles expand in the surrounding plasma where the scalar field still resides in the false vacuum. The expansion and collision of these bubbles (Kosowsky et al. 1992a, 1992b; Kosowsky & Turner 1993; Kamionkowski et al. 1994; Caprini et al. 2008; Huber & Konstandin 2008), together with sound waves generated in the plasma (Giblin & Mertens 2013, 2014; Hindmarsh et al. 2014, 2015), can source a primordial GWB (see Witten 1984; Hogan 1986 for seminal work). ⁸³ For earlier Bayesian searches for a phase transition signal in PTA data, see Arzoumanian et al. (2021) and Xue et al. (2021).

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., Ai et al. 2023; Krajewski et al. 2023; Li et al. 2023, for recent work on this topic). Therefore, in this work, we perform two separate analyses: a sound-wave-only analysis (pt-sound), where we assume that the runaway regime is not reached 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:

• 1.  
  T_*, the percolation temperature, i.e., the temperature of the universe when ∼34% of its volume has been converted to the true vacuum (Ellis et al. 2019a). For weak transitions, this temperature coincides with the temperature at the time of bubble nucleation, T_n ∼ T_*. Conversely, for supercooled transitions, we typically have T_n ≪ T_*. Barring the case of extremely strong transitions, α_* ⋙ 1 (see below), which we do not consider in this work, T_* also determines the reheating temperature after percolation, T_rh ∼ T_* (Ellis et al. 2019a).
• 2.  
  α_*, the strength of the transition, i.e., the ratio of the change in the trace of the energy–momentum tensor across the transition and the radiation energy density at percolation (Ellis et al. 2019b; Caprini et al. 2020).
• 3.  
  ${H}_{* }{R}_{* }={R}_{* }/{H}_{* }^{-1}$, the average bubble separation in units of the Hubble radius at percolation, ${H}_{* }^{-1}$. For relativistic bubble velocities, which is what we consider in the following, R_* is related to the bubble nucleation rate, β, by the relation H_* R_* = (8π)^1/3 H_*/β.

In addition to the parameters T_*, α_*, and H_* R_*, the GWB produced by a phase transition also depends on the velocity of the expanding bubble walls, v_w . However, deriving the precise value of this quantity is an open theoretical problem, which depends on model-dependent quantities, such as the strength of the interactions between the bubble walls and the SM plasma. Therefore, in our analysis, we fix the bubble velocity to unity (i.e., the speed of light in natural units). This assumption is well justified for strong phase transitions (Bodeker & Moore 2017), which, realistically, are the only ones that could lead to a detectable signal in our current data. In particular, we fix v_w = 1 for both phase transition scenarios that we are interested in, pt-sound and pt-bubble. In the latter case, v_w → 1 is automatically implied by the runaway behavior of the phase transition; in the former case, one actually expects a subluminal terminal velocity, v_w < 1. In this sense, our decision to fix v_w = 1 amounts to the optimistic assumption that this terminal velocity is numerically close to v_w = 1. A similar approach is followed by the authors of the LISA review paper (Caprini et al. 2020), who work with v_w = 0.95 throughout most of their analysis in the absence of more detailed microphysical calculations. Finally, we point out that the parameterization of the GWB signal in terms of H_* R_* = (8π)^1/3 v_w H_*/β already absorbs a large part of the dependence on the bubble wall velocity. The remaining v_w dependence is mostly contained in the efficiency factor κ_s (see below). However, in the regime of large α_* values, α_* ∼ 0.3 ⋯ 10, which turn out to be preferred by the NG15 data (see Figure 8), this dependence is rather weak (see Figure 13 in Espinosa et al. 2010), which justifies again keeping v_w fixed.

Zoom In Zoom Out Reset image size

The GWB spectrum sourced by bubbles and sound waves can be written in terms of these parameters as

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{b}(f)={ \mathcal D }\,{\tilde{{\rm{\Omega }}}}_{b}{\left(\displaystyle \frac{{\alpha }_{* }}{1+{\alpha }_{* }}\right)}^{2}{\left({H}_{* }{R}_{* }\right)}^{2}{ \mathcal S }(f/{f}_{b})\end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{s}(f)={ \mathcal D }\,{\tilde{{\rm{\Omega }}}}_{s}{\rm{\Upsilon }}({\tau }_{\mathrm{sw}}){\left(\displaystyle \frac{{\kappa }_{s}\,{\alpha }_{* }}{1+{\alpha }_{* }}\right)}^{2}\left({H}_{* }{R}_{* }\right){ \mathcal S }(f/{f}_{s}).\end{eqnarray} \tag{ 34 }$$

Here ${\tilde{{\rm{\Omega }}}}_{b}=0.0049$ (Jinno & Takimoto 2017) and ${\tilde{{\rm{\Omega }}}}_{s}=0.036$ (Hindmarsh et al. 2017), the efficiency factor ${\kappa }_{s}={\alpha }_{* }/(0.73+0.083\sqrt{{\alpha }_{* }}+{\alpha }_{* })$ (Espinosa et al. 2010) gives the fraction of the released energy that is transferred to plasma motion in the form of sound waves, and ${ \mathcal D }$ accounts for the redshift of the GW energy density,

$$\begin{eqnarray}&&{ \mathcal D }=\displaystyle \frac{{\pi }^{2}}{90}\displaystyle \frac{{T}_{0}^{4}}{{M}_{\mathrm{Pl}}^{2}{H}_{0}^{2}}\,{g}_{* }{\left(\displaystyle \frac{{g}_{* ,s}^{\mathrm{eq}}}{{g}_{* ,s}}\right)}^{4/3}\simeq 1.67\times {10}^{-5}.\end{eqnarray} \tag{ 35 }$$

We recall that T₀ and H₀ denote the photon temperature and Hubble rate today. The degrees of freedom g_* and g_*,s in Equation (35) are evaluated at T = T_*, and ${g}_{* ,s}^{\mathrm{eq}}$ is the number of degrees of freedom contributing to the radiation entropy at the time of matter–radiation equality. The production of GWs from sound waves stops after a period τ_sw, when the plasma motion turns turbulent (Ellis et al. 2019a, 2019b, 2020; Guo et al. 2021). In Equation (34), this effect is taken into account by the suppression factor

$$\begin{eqnarray}&&{\rm{\Upsilon }}({\tau }_{\mathrm{sw}})=1-{\left(1+2{\tau }_{\mathrm{sw}}{H}_{* }\right)}^{-1/2},\end{eqnarray} \tag{ 36 }$$

where the shock formation timescale, τ_sw, can be written in terms of the phase transition parameters as ${\tau }_{\mathrm{sw}}\approx {R}_{* }/{\bar{U}}_{f}$, with ${\bar{U}}_{f}^{2}\approx 3{\kappa }_{s}{\alpha }_{* }/[4(1+{\alpha }_{* })]$ (Weir 2018).

The functions ${{ \mathcal S }}_{b,s}$ describe the spectral shape of the GWB and are expected to peak at the frequencies

$$\begin{eqnarray}&&{f}_{b,s}\simeq 48.5\,\mathrm{nHz}\ {g}_{* }^{1/2}{\left(\displaystyle \frac{{g}_{* ,s}^{\mathrm{eq}}}{{g}_{* ,s}}\right)}^{1/3}\left(\displaystyle \frac{{T}_{* }}{1\,\mathrm{GeV}}\right)\ \displaystyle \frac{{f}_{b,s}^{* }{R}_{* }}{{H}_{* }{R}_{* }},\end{eqnarray} \tag{ 37 }$$

where the values of the peak frequencies at the time of GW emission are given by ${f}_{b}^{* }=0.58/{R}_{* }$ (Jinno & Takimoto 2017) and ${f}_{s}^{* }=1.58/{R}_{* }$ (Hindmarsh et al. 2017). In passing, we mention that the numerical factors in these estimates may still change in the future, as our understanding of cosmological phase transitions improves. However, at the level of our Bayesian fit analysis, changes in these prefactors can be absorbed in the temperature scale T_*, which in its role as an independent fit parameter only controls the peak frequencies in Equation (37). A similar argument applies to the numerical factors in Equations (33) and (34): changes in these prefactors can always be absorbed in a rescaled version of the fit parameter α_*, which only appears in the expressions for the peak amplitudes of the GWB signal.

The shape of the spectral functions can be usually approximated with a broken power law of the form

$$\begin{eqnarray}&&{ \mathcal S }(x)=\displaystyle \frac{1}{{ \mathcal N }}\displaystyle \frac{{\left(a+b\right)}^{c}}{{\left({{bx}}^{-a/c}+{{ax}}^{b/c}\right)}^{c}}.\end{eqnarray} \tag{ 38 }$$

Here a and b are two real and positive numbers that give the slope of the spectrum at low and high frequencies, respectively; c parameterizes the width of the peak. The normalization constant, ${ \mathcal N }$, ensures that the logarithmic integral of ${ \mathcal S }$ is normalized to unity and is given by

$$\begin{eqnarray}&&{ \mathcal N }={\left(\displaystyle \frac{b}{a}\right)}^{a/n}{\left(\displaystyle \frac{{nc}}{b}\right)}^{c}\ \displaystyle \frac{{\rm{\Gamma }}\left(a/n\right){\rm{\Gamma }}\left(b/n\right)}{n\,{\rm{\Gamma }}(c)},\end{eqnarray} \tag{ 39 }$$

where n = (a + b)/c and Γ denotes the gamma function. While the values of the coefficients a, b, and c can in principle be estimated from numerical simulations, we allow them to float within the prior ranges given in Table 3. These prior ranges were chosen to roughly reflect the typical uncertainties of numerical simulations and any possible model dependency of these coefficients (see, e.g., Hindmarsh et al. 2017, 2021; Cutting et al. 2018, 2021; Lewicki & Vaskonen 2020, 2021). ⁸⁴

The reconstructed posterior distributions for the parameters α_*, T_*, and H_* R_* of the pt-sound and pt-bubble models are reported in 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). ⁸⁵ Corner plots including the posterior distributions for the spectral shape parameters a, b, c and SMBHB parameters A_BHB and γ_BHB are reported in Figures 26 and 27 in Appendix C.3.

In all analyses, the data prefer a relatively strong and slow phase transition. Specifically, for pt-bubble, we find α_* > 1.1 (0.29) and H_* R_* > 0.28 (0.14) at the 68% (95%) credible level. When the SMBHB signal is added on top of the GWB predicted by pt-bubble, we find α_* > 1.0 (0.23) and H_* R_* > 0.26 (0.11) at the 68% (95%) credible level. For the pt-sound model, we find α_* > 0.42 (0.37) and H_* R_* ∈ [0.053, 0.27] ([0.046, 0.89]) at the 68% (95%) credible level. Including the SMBHB signal on top of the one predicted by pt-sound, we find α_* ∈ [0.46, 5.4] (>0.16) and H_* R_* ∈ [0.054, 0.35] (>0.0015) at the 68% (95%) credible level.

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] ([0.023, 1.75]) GeV at the 68% (95%) credible level, whereas for the pt-sound model we get T_* ∈ [4.7, 33] ([2.7, 93]) MeV at the 68% (95%) credible level. The shift between these T_* intervals is partially explained by the different numerical factors in the frequencies ${f}_{s}^{* }$ and ${f}_{b}^{* }$ (see Equation (37)). As explained below Equation (37), any change in these numerical factors can be reabsorbed in a redefinition of the fit parameter T_*.

The inclusion of the SMBHB signal, by adding power to the lowest frequency bins, allows the T_* posterior for the pt-sound model to extend to higher values. In this case, we find that T_* ∈ [4.9, 50] ([0.8, 2 × 10⁶]) MeV at the 68% (95%) credible level. Here the increase in the 68% upper limit is reflected in the slight shift between the red and blue dashed vertical lines in the 1D marginalized posterior distribution for T_* in the right panel of 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] ([0.017, 3.27]) GeV at the 68% (95%) credible level.

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}^{* }$, but also of the lower prior range for the low-frequency spectral index adopted for the pt-bubble model. Indeed, a shallower low-frequency tail allows spectra with a higher peak frequency to still produce a sizable signal in the lowest frequency bins. In Appendix C.3, we report the results of the analysis in which the low-frequency slope is set to the value predicted by causality (a = 3). In this case, as expected, the reconstructed phase transition temperatures for the two phase transition models are closer to each other.

The corner plots in Figure 8 also illustrate that, as expected from the expression for the peak frequency in Equation (37), there is an approximately linear correlation between ${\mathrm{log}}_{10}{T}_{* }$ and ${\mathrm{log}}_{10}{H}_{* }{R}_{* }$. For α_* ≲ 1, we instead find α_* ∼ 1/(H_* R_*) for the pt-bubble model and ${\alpha }_{* }^{2}\sim 1/({H}_{* }{R}_{* })$ for the pt-sound model as expected from the factors in Equations (33) and (34).

We also notice that, for both models, the posterior distribution for T_* is peaked at significantly larger values compared to what was derived in the 12.5 yr analysis (Arzoumanian et al. 2021). This shift results from the reconstructed h²Ω_GW spectrum being bluer than the one derived for the common process observed in the 12.5 yr data set. As a result, the lowest frequency bins, which were fit by the high-frequency tail of the phase transition spectrum in the 12.5 yr analysis, are now fit by the low-frequency tail of the spectrum. This then translates into a higher peak frequency and therefore a higher transition temperature. Incidentally, the larger reconstructed value for the transition temperature allows the phase transition signal to safely evade bounds from BBN and CMB observations (Bai & Korwar 2022; Deng & Bian 2023) for both the pt-bubble and pt-sound models, which constrain the phase transition parameter space at temperatures around T_* ∼ 1 MeV.

Instead, we conclude that the reconstructed posterior distribution of T_* is compatible with phase transition scenarios that have been discussed in the literature as a possible source of GWs in the PTA band: (i) BSM models in which the chiral-symmetry-breaking phase transition in quantum chromodynamics (QCD) is a strong first-order phase transition (see, e.g., Li et al. 2021; Neronov 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. 2021; Ratzinger & Schwaller 2021). 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. 2017).

Finally, we report that, like the models studied in Sections 5.1 and 5.2, 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 ${ \mathcal B }=18.1\pm 0.6$ and ${ \mathcal B }=3.7\pm 0.1$, respectively, while the Bayes factors for pt-bubble+smbhb and pt-sound+smbhb versus smbhb are ${ \mathcal B }=12.6\pm 0.5$ and ${ \mathcal B }=6.5\pm 0.3$, respectively. 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 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, which can as be as low as a = 1 for pt-bubble, but which we require to be at least a = 3 for pt-sound (see 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 ${\mathrm{log}}_{10}{A}_{\mathrm{BHB}}\sim -(15\cdots 14)$ and γ_BHB ∼ 3 ⋯ 4 (see Figure 27).

Cosmic strings are effectively 1D topological defects that can form in the early universe as a consequence of a cosmological phase transition (Kibble 1976). 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) symmetry. The breaking of global U(1) symmetries results in the formation of global-string networks, which happens, e.g., in axion models. The GW signal in this case is suppressed because global strings lose most of their energy by emitting light pseudo-Nambu–Goldstone bosons (i.e., "axions"). In the following, we therefore focus on local strings from the spontaneous breaking of a local U(1) symmetry as occurs in many particle physics models of the early universe, such as grand unified theories (GUTs), where intermediate symmetry breaking stages can be realized in the form of string-producing phase transitions. Cosmic strings are thus well motivated by ideas about particle physics at very high energies, which prompts us to consider them as yet another possible source of exotic GWs (Vilenkin 1981a; Hogan & Rees 1984; Vachaspati & Vilenkin 1985; Accetta & Krauss 1989; 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. (2021), Ellis & Lewicki (2021), and Blanco-Pillado et al. (2021).

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, μ. In the case of metastable strings, there is a second parameter, κ, which controls the lifetime of the network. For superstrings, the relevant parameters are μ and the intercommutation probability P (see below). Before we move on, we note that, as an alternative approach to the Nambu–Goto approximation, it is possible to describe cosmic strings in a field-theoretical language, e.g., in terms of the Abelian-Higgs model. The comparison of these two approaches, i.e., the relation between Nambu–Goto strings and Abelian–Higgs strings, is the subject of ongoing research (Blanco-Pillado et al. 2023), and we refer the reader to Section 3.4 of Auclair et al. (2020) 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. 2020).

Shortly after their formation, cosmic strings enter a scaling regime where the total energy stored in the network, ρ_cs, remains a constant fraction of the critical energy density, ${{\rm{\Omega }}}_{\mathrm{cs}}={\rho }_{\mathrm{cs}}/{\rho }_{c}\approx {\rm{const}}$ (Hindmarsh & Kibble 1995). 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,

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{\mathrm{GW}}^{\mathrm{cs}}\left(f\right)=\displaystyle \frac{8\pi }{3{H}_{0}^{2}}{\left(G\mu \right)}^{2}\displaystyle \sum _{k=1}^{{k}_{\max }}{P}_{k}\,{{ \mathcal I }}_{k}\left(f\right).\end{eqnarray} \tag{ 40 }$$

Here the dimensionless factor G μ denotes the cosmic-string tension in units of Newton's constant, for which we choose a log-uniform prior in our numerical analysis, ${\mathrm{log}}_{10}G\mu \in [-14,-6]$ for stable strings and superstrings and ${\mathrm{log}}_{10}G\mu \in [-14,-1.5]$ for metastable strings. The sum in Equation (40) runs over the harmonic excitations of the closed string loops that, given a loop of length ℓ, correspond to oscillations at frequency f = 2k/ℓ. We evaluate the sum starting at the fundamental oscillation mode, k = 1, and including terms up to ${k}_{\max }={10}^{5}$, which ensures good convergence of the GW spectrum.

The dimensionless factor P_k inside the sum describes the GW power, in units of G μ², that is emitted by a loop in its kth excitation. For large ${k}_{\max }$, we can write

$$\begin{eqnarray}&&{P}_{k}=\displaystyle \frac{{\rm{\Gamma }}}{\zeta \left(q\right)}\displaystyle \frac{1}{{k}^{q}},\end{eqnarray} \tag{ 41 }$$

where the prefactor ${\rm{\Gamma }}/\zeta \left(q\right)$ ensures that the total power emitted in GWs amounts to ∑_k P_k = Γ and where the power-law index q depends on the predominant source of GWs on the loops. In our analysis, we set Γ = 50, as suggested by numerical simulations (Blanco-Pillado & Olum 2017), and discuss different possibilities for q. The actual average GW spectrum from non-self-intersecting loops is still uncertain. We therefore choose several different models that give an idea of the range of possibilities. Specifically, we consider four different models of stable cosmic strings:

stable-c: Stable strings emitting GWs only in the form of GW bursts from cusps on closed loops, q = 4/3 (Vachaspati & Vilenkin 1985).

stable-k: Stable strings emitting GWs only in the form of GW bursts from kinks on closed loops, q = 5/3 (Damour & Vilenkin 2001).

stable-m: Stable strings emitting monochromatic GWs at the fundamental frequency f = 2/ℓ of closed loops.

stable-n: Stable strings described by numerical simulations including GWs from cusps and kinks (Blanco-Pillado et al. 2011, 2015).

For stable-n, P_k is dominated by cusp emission at large k, i.e., q ≈ 4/3 for k ≳ 100, while at smaller k it deviates from the pure cusp spectrum, reaching q values of up to q ∼ 1.5 around k ∼ 10. Meanwhile, Equation (41) is irrelevant for stable-m; for this model, we simply set P_k = Γ if k = 1 and P_k = 0 otherwise. More details on our four stable-string models can be found in Blanco-Pillado et al. 2021.

Finally, the frequency-dependent factor ${{ \mathcal I }}_{k}$ in Equation (40) represents an integral of the number density of closed string loops, n_l (ℓ, t), over all possible GW emission times,

$$\begin{eqnarray}&&{{ \mathcal I }}_{k}\left(f\right)=\displaystyle \frac{2k}{f}{\int }_{{t}_{\mathrm{ini}}}^{{t}_{0}}{dt}{\left(\displaystyle \frac{a\left(t\right)}{a({t}_{0})}\right)}^{5}{n}_{l}\left(\displaystyle \frac{2k}{f}\displaystyle \frac{a\left(t\right)}{a({t}_{0})},t\right).\end{eqnarray} \tag{ 42 }$$

Here t₀ is the present time and t_ini is the time when the network reaches the scaling attractor solution. The precise value of t_ini only affects the high-frequency part of the GW spectrum and plays no role in our analysis. The loop number density n_l can be estimated based on the velocity-dependent one-scale (VOS) model,

$$\begin{eqnarray}&&{n}_{l}\left({\ell },t\right)={ \mathcal F }\ \displaystyle \frac{{C}_{* }\,{\rm{\Theta }}\left(t-{t}_{* }\right){\rm{\Theta }}\left({t}_{* }-{t}_{\mathrm{ini}}\right)}{{\alpha }_{* }\left({\alpha }_{* }+{\rm{\Gamma }}G\mu +{\dot{\alpha }}_{* }{t}_{* }\right){t}_{* }^{4}}{\left(\displaystyle \frac{{a}_{* }}{a\left(t\right)}\right)}^{3},\end{eqnarray} \tag{ 43 }$$

where ${ \mathcal F }=0.1$ is an efficiency factor (Blanco-Pillado et al. 2014; Auclair et al. 2020) 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_*:

$$\begin{eqnarray}&&{t}_{* }=\displaystyle \frac{{\ell }+{\rm{\Gamma }}G\mu \,t}{{\alpha }_{* }+{\rm{\Gamma }}G\mu },\qquad {\alpha }_{* }=\alpha \left({t}_{* }\right).\end{eqnarray} \tag{ 44 }$$

The time-dependent functions C and α characterize the efficiency of loop formation from the network and the typical loop size at the time of production, respectively,

$$\begin{eqnarray}&&C\left(t\right)=\displaystyle \frac{\tilde{c}}{\sqrt{2}}\displaystyle \frac{\bar{v}\left(t\right)}{{\xi }^{3}\left(t\right)},\qquad \alpha \left(t\right)={\alpha }_{L}\,\xi \left(t\right).\end{eqnarray} \tag{ 45 }$$

Here $\tilde{c}$ is the loop-chopping parameter and α_L controls the loop length at the time of production in units of the correlation length of the long-string network, L = ξ t. We set $\tilde{c}\simeq 0.23$ and α_L ≃ 0.37, in agreement with numerical simulations (Martins & Shellard 2002; Blanco-Pillado et al. 2011, 2014). With $\tilde{c}$ fixed, we are able to solve the VOS equations for the dimensionless correlation length ξ and the rms velocity of the long strings, $\bar{v}$, and hence determine the time dependence of C and α. In doing so, we account for the exact evolution of the scale factor, a, in ΛCDM, including the temperature-dependent variation in the number of relativistic degrees of freedom. At very high temperatures, this analysis returns ξ_r ≃ 0.27 and ${\bar{v}}_{r}\simeq 0.66$, such that α ≃ 0.10 deep in the radiation-dominated era. The loop number density n_l obtained in this way agrees very well with the result of numerical simulations (Blanco-Pillado et al. 2014) in the limit of constant g_* and g_*,s .

In the case of metastable strings, the loop number density in Equation (43) receives two correction factors,

$$\begin{eqnarray}&&{n}_{l}^{\mathrm{meta}}\left({\ell },t\right)={\rm{\Theta }}\left({t}_{s}-{t}_{* }\right)E\left({\ell },t\right){n}_{l}\left({\ell },t\right).\end{eqnarray} \tag{ 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

$$\begin{eqnarray}&&{t}_{s}=\displaystyle \frac{1}{{{\rm{\Gamma }}}_{d}^{1/2}},\end{eqnarray} \tag{ 47 }$$

with the decay rate, Γ_d , counting the number of monopole nucleation events per time and per unit string length,

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{d}=\displaystyle \frac{\mu }{2\pi }\,{e}^{-\pi \kappa }.\end{eqnarray} \tag{ 48 }$$

Here κ is a measure for the ratio of the GUT and U(1) symmetry breaking scales in the underlying GUT model. Specifically, $\sqrt{\kappa }$ describes the ratio of the GUT monopole mass and the square root of the U(1) string tension,

$$\begin{eqnarray}&&\sqrt{\kappa }=\displaystyle \frac{{m}_{\mathrm{GUT}}}{{\mu }^{1/2}}\sim \displaystyle \frac{{{\rm{\Lambda }}}_{\mathrm{GUT}}}{{{\rm{\Lambda }}}_{U(1)}}.\end{eqnarray} \tag{ 49 }$$

Meanwhile, the second new factor in Equation (46) represents an exponential suppression term, reflecting the depletion of the loop number density in the decaying network,

$$\begin{eqnarray}&&E\left({\ell },t\right)={e}^{-{{\rm{\Gamma }}}_{d}\left[{\ell }\left(t-{t}_{* }\right)+\displaystyle \frac{1}{2}\,{\rm{\Gamma }}G\mu {\left(t-{t}_{* }\right)}^{2}\right]}.\end{eqnarray} \tag{ 50 }$$

As evident from Equations (48) and (50), the loop number density of a metastable network is exponentially sensitive to the decay parameter κ. For $\sqrt{\kappa }\gg 10$, the lifetime of the network exceeds the age of the universe, such that the resulting GW signal is indistinguishable from the signal of a stable-string network. For $\sqrt{\kappa }\sim 1$, on the other hand, the network decays very fast in the early universe, such that no GW signal at PTA frequencies is generated. For these reasons, we choose a uniform prior on $\sqrt{\kappa }$ in a rather narrow range, $\sqrt{\kappa }\in [7.0,9.5]$, which is enough to capture the nontrivial aspects of the metastable scenario.

Monopole nucleation in a metastable network results in string segments, dumbbell-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

$$\begin{eqnarray}&&{SU}{(2)}_{1}\times U{(1)}_{2}\to U{(1)}_{1}\times U{(1)}_{2}\to U(1).\end{eqnarray} \tag{ 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 C.4. Moreover, we need to specify the power spectrum of GWs emitted by string segments. Following Buchmuller et al. (2021), we use again Equation (41) and set q = 1, which can be analytically derived in the straight-string approximation (Martin & Vilenkin 1997). At the same time, we assume a pure cusp spectrum (q = 4/3) for the GW emission from loops in both meta-l and meta-ls.

Finally, we also consider cosmic superstrings:

super: Cosmic superstrings; suppressed intercommutation probability P; GWs from cusps on loops, q = 4/3.

Cosmic superstrings are characterized by a reduced intercommutation probability P as a consequence of their quantum-mechanical nature. In superstring networks, intercommutation events are, therefore, less frequent, which demands longer numerical simulations (potentially up to a factor 10³ longer). Such simulations have not yet been carried out, which makes the prediction of the GW signal from cosmic superstrings rather uncertain at present. In the absence of a more sophisticated treatment, we follow the standard practice in the literature and simply estimate the GW spectrum from cosmic superstrings by rescaling the spectrum in Equation (40),

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{\mathrm{GW}}^{\mathrm{cs}}\left(f\right)\quad \to \quad \displaystyle \frac{1}{P}\,{{\rm{\Omega }}}_{\mathrm{GW}}^{\mathrm{cs}}\left(f\right),\end{eqnarray} \tag{ 52 }$$

where we allow for P values drawn from a log-uniform prior, ${\mathrm{log}}_{10}P\in \left[-3,0\right]$, which covers the theoretically expected range for cosmic F- and D-superstrings (Jackson et al. 2005).

We now turn to the outcome of our Bayesian fit analysis, discussing first our results for stable strings. As evident from Figures 2 and 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 ${ \mathcal B }=0.277\pm 0.006$, ${ \mathcal B }=0.364\pm 0.008$, ${ \mathcal B }=0.379\,\pm 0.008$, and ${ \mathcal B }=0.307\pm 0.006$, respectively; comparing the stable-c+smbhb, stable-k+smbhb, stable-m+smbhb, and stable-n+smbhb models to the smbhb model, we obtain ${ \mathcal B }=0.76\pm 0.01$, ${ \mathcal B }=0.89\pm 0.02$, ${ \mathcal B }=0.84\pm 0.02$, and ${ \mathcal B }=0.83\pm 0.01$, respectively. These Bayes factors are close to unity, which means that adding a cosmic-string contribution to the GWB signal on top of the SMBHB signal does not improve the fit. Meanwhile, the larger prior volume pushes the Bayes factors to values slightly smaller than unity.

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 μ is necessary. Large G μ values, however, tend to result in a rather flat GW spectrum in the PTA band as seen in Figure 20 (see also, e.g., Figure 1 in Blanco-Pillado et al. 2011 or Figure 1 in Blasi et al. 2021), which clashes with the fact that the data seem to prefer a blue-tilted h²Ω_GW spectrum. If we approximate the GW signal from stable strings by a simple power law, this observation can be reformulated in the language of the γ–A plot in Figure 1. In terms of γ and A, we then conclude that stable strings allow one to obtain γ ≲ 4 only for ${\mathrm{log}}_{10}A\lesssim -15.0$, while larger values in the range ${\mathrm{log}}_{10}A\sim -\left(14.5\cdots 14.0\right)$ can only be achieved for γ ∼ 5. The GW spectrum from stable strings is therefore always either too weak or too flat.

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-SMBHBs models. In the former case, we find peaked distributions centered around values of the order of ${\mathrm{log}}_{10}G\mu \sim -\left(10.5\cdots 10.0\right)$. Values of the cosmic-string tension in this range represent a compromise between the two extremes discussed above: at smaller ${\mathrm{log}}_{10}G\mu$ the GW signal becomes too weak, while at larger ${\mathrm{log}}_{10}G\mu$ it becomes too flat. Moreover, we note that the order of the peaks in the posterior distributions—stable-m, stable-k, stable-n, stable-c from left to right—agrees with the ordering found in Blanco-Pillado et al. (2021) at γ ≲ 4.6. The underlying reason is that, for fixed G μ, stable-m predicts the strongest GW signal in the nanohertz frequency band, followed by stable-k and stable-n, while stable-c predicts the weakest signal.

Zoom In Zoom Out Reset image size

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 ${\mathrm{log}}_{10}G\mu$ range but also an extended plateau at small ${\mathrm{log}}_{10}G\mu$ values. This plateau corresponds precisely to the plateau that we alluded to in the definition of the K ratio in Equation (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: ${\mathrm{log}}_{10}G\mu \lt -9.67$, −9.87, −10.10, and −9.71. We recall that, for cosmic-string tensions exceeding these values, the likelihood of the strings-plus-SMBHBs model drops below 1/10 of the likelihood of the smbhb benchmark model. In this region of parameter space, adding GWs from stable strings to the signal does more harm than good and is strongly disfavored by the data. For comparison, we also quote the upper limits of the 95% Bayesian credible intervals for G μ, specifically of the 95% highest posterior density intervals (HPDIs). For each model, the HPDI is the narrowest possible range that includes 95% of the posterior probability. By construction, the boundaries of an HPDI are at points with the same posterior probability density, and each parameter value inside the HPDI has higher posterior probability density than any parameter value outside the HPDI. ⁸⁶ We find ${\mathrm{log}}_{10}G\mu \lt -9.88$, −10.04, −10.26, and −9.90 for stable-c+smbhb, stable-k+smbhb, stable-m+smbhb, and stable-n+smbhb, respectively.

Next, we turn to metastable strings, which, as can be read off from 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 ${ \mathcal B }=13.4\pm 0.4$ and ${ \mathcal B }=21.3\pm 0.8$, respectively. Adding an SMBHB contribution to the GW signal, the Bayes factors for the meta-l and meta-ls models are instead ${ \mathcal B }=11.1\pm 0.3$ and ${ \mathcal B }=18.9\pm 0.7$, respectively. Again, the fit does not become better, but the larger prior volume results in a decrease of the Bayes factors.

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 }$ of around $\sqrt{\kappa }\sim 8$ in combination with a large cosmic-string tension, ${\mathrm{log}}_{10}G\mu \sim -\left(7\cdots 4\right)$. Here $\sqrt{\kappa }\sim 8$ causes a suppression of the GW signal in the nanohertz band compared to stable strings because of the exponential factor in Equation (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.

Zoom In Zoom Out Reset image size

Let us comment on a few characteristic features of the posterior distributions shown in Figure 10. First, we find that the ${\mathrm{log}}_{10}G\mu$ posterior in the meta-ls model extends to slightly larger values than in the meta-l model, which demonstrates the effect of the additional contribution to the GW signal from string segments. In fact, in the region of highest posterior density, the segment contribution dominates over the loop contribution in the meta-ls model. Second, we point out that the 1D marginalized posterior distributions for ${\mathrm{log}}_{10}G\mu$ exhibit small local maxima at ${\mathrm{log}}_{10}G\mu \sim -\left(11\cdots 10\right)$, which correspond to the stable-string limit within the metastable-string models. This limit is realized for large values of the decay parameter, $\sqrt{\kappa }\gtrsim 9$, which pushes the effect of the network decay to frequencies below the PTA band. Next, we observe that for meta-l the ${\mathrm{log}}_{10}G\mu$ posterior experiences a sharp drop-off at ${\mathrm{log}}_{10}G\mu \sim -5$, whereas for meta-ls there is a small dip in the ${\mathrm{log}}_{10}G\mu$ posterior at ${\mathrm{log}}_{10}G\mu \sim -5$. Both features can be traced back to the Heaviside theta function in Equation (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 ${\mathrm{log}}_{10}G\mu$ above ${\mathrm{log}}_{10}G\mu \sim -5$. In this sense, the drop-off and the dip in the ${\mathrm{log}}_{10}G\mu$ posteriors should be regarded as being due to our simplified theoretical modeling of the GW spectrum. More work is necessary to improve on the description in terms of a simple Heaviside theta function and obtain a better understanding of the evolution of the decaying network. We expect that a more accurate description of the transition from the scaling to the decay regime would result in smoother ${\mathrm{log}}_{10}G\mu$ posteriors. Finally, we note that some of the fluctuations in the posteriors in 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 Equation (40).

We also use the corner plots in Figure 10 to highlight the relevant bounds on the parameter space spanned by ${\mathrm{log}}_{10}G\mu$ and $\sqrt{\kappa }$. The gray shaded areas notably indicate the K-ratio bound, which marks the parameter regions that are ruled out by the NANOGrav data. In these regions, the GW signal from metastable strings exceeds the observed signal and hence is unacceptably large. We find that the NANOGrav bound is stronger than the well-known CMB bound, which demands that a cosmic-string network that has not yet decayed by the time of recombination must not have a cosmic-string tension larger than ${\mathrm{log}}_{10}G\mu \simeq -7$ (Ade et al. 2016; Charnock et al. 2016; Lizarraga et al. 2016). In order to derive the CMB bound, we estimate that a decaying network completely disappears because of GW emission at times around ${t}_{e}\simeq {\left(2/\left({\rm{\Gamma }}G\mu \right)\right)}^{1/2}{t}_{s}$, which is the time when the second term in the exponent in Equation (50) becomes large. The teal shaded areas in Figure 10 indicate where the conditions ${\mathrm{log}}_{10}G\mu \gt -7$ and t_e > t_rec ≃ 370,000 yr are satisfied simultaneously. Last but not least, the yellow solid lines represent the LVK bound in Equation (24), which translates to the upper limit G μ ≲ 2 × 10⁻⁷. ⁸⁷

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 ${\mathrm{log}}_{10}G\mu$–$\sqrt{\kappa }$ parameter plane are mostly ruled out by this bound. However, we stress that the LVK bound in Figure 10 relies on an extrapolation of the GW signal across 10 orders of magnitude in frequency, from PTA frequencies, f ∼ few × 10⁻⁹ nHz, to LVK frequencies, f ∼ few × 10 Hz. In order to perform this extrapolation, we have to assume a cosmological expansion history across 10 orders of magnitude in temperature, even though not much is known about the equation of state of the universe prior to BBN. Any deviation from the expansion history of standard big bang cosmology can therefore affect the extrapolation of the GW spectrum to higher frequencies and potentially render the LVK bound harmless. A simple example is an early stage of matter domination (Allahverdi et al. 2021), 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 $\sqrt{\kappa }\sim 8$ and ${\mathrm{log}}_{10}G\mu \sim -7$ are still able to provide a good fit of the data, which represents a particularly interesting scenario for two reasons: a value of ${\mathrm{log}}_{10}G\mu \sim -7$ would point to a GUT origin of the string network, and ground-based interferometers would be poised to detect the stochastic GWB signal from such a network in the near future.

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: ${ \mathcal B }=46\pm 2$ for GWs from superstrings alone and ${ \mathcal B }=37\pm 2$ for an SMBHB contribution added to the superstring GW signal, both compared to the smbhb model.

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, ${\mathrm{log}}_{10}G\mu \lesssim -11$, yield a blue-tilted h²Ω_GW spectrum at nanohertz frequencies (see our discussion above). In the case of ordinary strings, the amplitude of this spectrum would be too weak; however, for superstrings the suppression of the GW signal at ${\mathrm{log}}_{10}G\mu \lesssim -11$ can be compensated by the 1/P enhancement factor in Equation (52). By choosing P appropriately, the amplitude of the GW signal can then be raised until a good fit of the data is reached. This interplay between the two parameters of the super model can also be observed in Figure 11. The 2D posterior distribution for ${\mathrm{log}}_{10}G\mu$ and P in this corner plot displays a strong covariance in line with our heuristic understanding.

Zoom In Zoom Out Reset image size

The highest posterior density is achieved at small intercommutation probabilities and cosmic-string tensions, ${\mathrm{log}}_{10}P\sim -3$ and ${\mathrm{log}}_{10}G\mu \sim -12$, where ${\mathrm{log}}_{10}P=-3$ corresponds in fact to the edge of the prior range for P. This parameter region is in tension with the LVK bound in Equation (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 ratio of less than 1/10. As expected, this region is located to the lower right of the 68% and 95% credible regions, where small ${\mathrm{log}}_{10}P$ values cause the amplitude of the GW spectrum to exceed the measured strength of the signal. The tension of cosmic superstrings can also be constrained by CMB observations (Charnock et al. 2016). Existing upper limits are, however, located at larger G μ values and thus not relevant for our results in 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 are 2D topological defects that form when a cosmological phase transition results in the spontaneous breaking of a discrete symmetry, filling the universe with patches in different degenerate vacua (Kibble 1976). In the absence of significant interactions with relativistic particles, domain wall networks are expected to reach a scaling regime in which each Hubble volume, H⁻³, contains ${ \mathcal A }\sim { \mathcal O }(1)$ domain walls (see, e.g., Press et al. 1989; Hindmarsh 1996, 2003; Garagounis & Hindmarsh 2003). In this regime, the energy density stored in domain walls is given by

$$\begin{eqnarray}&&{\rho }_{\mathrm{DW}}\sim { \mathcal A }\,\sigma H,\end{eqnarray} \tag{ 53 }$$

where ${ \mathcal A }$ is nearly constant and σ is the domain wall tension, which gives the domain walls their surface energy density. During the scaling regime, the energy density stored in domain walls dilutes more slowly than that of relativistic radiation and nonrelativistic matter. Indeed, the total energy density of the background always scales like ρ_c ∝ H², which means that Ω_DW = ρ_DW/ρ_c grows like the Hubble radius when the domain wall network is in the scaling regime, Ω_DW ∝ H⁻¹. Therefore, domain walls eventually overclose the universe and alter the cosmological evolution in a way that is incompatible with CMB observations (Zeldovich et al. 1974). A possible solution to this problem is to have domain walls decay at some temperature T_* by assuming, e.g., that the global symmetry responsible for domain wall formation is explicitly broken. Indeed, an explicit breaking of the symmetry introduces a bias among the possible low-energy vacua, thus lifting their degeneracy, which eventually leads to the collapse of the domain wall network (Vilenkin 1981b; Gelmini et al. 1989; Larsson et al. 1997).

During the scaling regime, domain walls continuously change their configuration and shrink in order to maintain the scaling relation in Equation (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 (Vilenkin 1981b; Preskill et al. 1991; Gleiser & Roberts 1998; Chang et al. 1999; Hiramatsu et al. 2010; Kawasaki & Saikawa 2011)

$$\begin{eqnarray}&&{h}^{2}{{\rm{\Omega }}}_{\mathrm{GW}}(f)=\displaystyle \frac{3}{32\pi }{ \mathcal D }\,\tilde{\epsilon }\,{\alpha }_{* }^{2}\,{ \mathcal S }(f/{f}_{p}),\end{eqnarray} \tag{ 54 }$$

where ${ \mathcal D }$ is the dilution factor defined in Equation (35), $\tilde{\epsilon }=0.7$ is an efficiency coefficient derived from numerical simulations (Hiramatsu et al. 2014), and α_* is the fraction of the total energy density stored in domain walls at T_*,

$$\begin{eqnarray}&&{\alpha }_{* }\equiv {{\rm{\Omega }}}_{\mathrm{DW}}\left({T}_{* }\right)={\left.\displaystyle \frac{{\rho }_{\mathrm{DW}}}{3{H}^{2}{M}_{\mathrm{Pl}}^{2}}\right|}_{T={T}_{* }}.\end{eqnarray} \tag{ 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_*, i.e., the Hubble scale at T_*. Therefore, after redshifting until today, the GWB is expected to exhibit a peak frequency

$$\begin{eqnarray}&&{f}_{p}=1.14\,\mathrm{nHz}{\left(\displaystyle \frac{10.75}{{g}_{* ,s}}\right)}^{1/3}{\left(\displaystyle \frac{{g}_{* }}{10.75}\right)}^{1/2}\left(\displaystyle \frac{{T}_{* }}{10\,\mathrm{MeV}}\right),\end{eqnarray} \tag{ 56 }$$

with both g_* and g_*,s evaluated at T_*. For the spectral shape, ${ \mathcal S }(x)$, we use a parameterization similar to the one used for the phase transition spectrum in Section 5.3: ⁸⁸

$$\begin{eqnarray}&&{ \mathcal S }(x)=\displaystyle \frac{{\left(a+b\right)}^{c}}{{\left({{bx}}^{-a/c}+{{ax}}^{b/c}\right)}^{c}}.\end{eqnarray} \tag{ 57 }$$

Causality fixes a = 3, while numerical simulations (Hiramatsu et al. 2014) suggest b ≃ c ≃ 1. In our analysis, we therefore fix the low-frequency slope of the domain wall signal to the value predicted by causality, and, following Ferreira et al. (2022), we allow b and c to float within the uniform prior ranges b ∈ [0.5, 1] and c ∈ [0.3, 3].

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 α_*, we parameterize the strength of the domain wall signal in terms of the dark radiation produced in the decay, ρ_DR. This quantity is usually expressed in terms of the effective number of extra neutrino species, ΔN_eff, which we defined in Section 5.1 and which is bounded from above by BBN and CMB observations, ${\rm{\Delta }}{N}_{\mathrm{eff}}^{\max }\sim {\rm{few}}\times 0.1$ (Pisanti et al. 2021; Yeh et al. 2021). In our MCMC analysis of the dw-dr model, we follow Ferreira et al. (2022) and use ${\rm{\Delta }}{N}_{\mathrm{eff}}^{\max }=0.39$ as the upper bound on our prior for the parameter ΔN_eff. Assuming that all the domain wall energy is converted to dark radiation at T_*, ΔN_eff is related to α_* by

$$\begin{eqnarray}&&{\rm{\Delta }}{N}_{\mathrm{eff}}\simeq 0.62{\left(\displaystyle \frac{10.75}{{g}_{* }}\right)}^{1/3}\left(\displaystyle \frac{{g}_{* ,s}}{{g}_{* }}\right)\left(\displaystyle \frac{{\alpha }_{* }}{0.1}\right),\end{eqnarray} \tag{ 58 }$$

where, as before, both g_* and g_*,s are evaluated at T_*.

In the dw-sm model, BBN restricts the possible values of the decay temperature to T_* ≳ 2.7 MeV (Jedamzik 2006; Bai & Korwar 2022) for any detectable value of α_*. Following Ferreira et al. (2022), we also impose α_* < 0.3 to avoid any possible deviation from radiation domination and to evade bounds from ΔN_eff.

The reconstructed posterior distributions for the parameters T_* and α_* (T_* and ΔN_eff) of the dw-sm (dw-dr) model are reported in 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 and c are reported in Figure 31 in Appendix C.5.

Zoom In Zoom Out Reset image size

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⁻⁸ Hz such that most of the low frequency bins are fit by the low-frequency tail of the spectrum (see Figures 3 and 19). Specifically, for the dw-sm (dw-dr) model, we find that T_* ∈ [110, 275] ([79, 153]) MeV at the 68% credible level and T_* ∈ [76, 505] ([54, 198]) MeV at the 95% credible level. When including the SMBHB contribution to the GWB, we find T_* ∈ [108, 309] ([54, 216]) MeV at the 68% credible level and T_* ∈ [67, 843] (no bound on T_*) MeV at the 95% credible level for the dw-sm (dw-dr) model. We notice that, with and without the inclusion of SMBHBs, the recovered transition temperature for the dw-sm model is high enough to evade BBN constraints.

For the dw-dr model, the posterior distribution for ΔN_eff is peaked near the upper prior boundary, signaling that larger values fit the observed signal better. Specifically, we find ΔN_eff ≳ 0.32 at the 68% credible level and ΔN_eff ≳ 0.25 at the 95% credible level. Including the contribution from SMBHBs allows the distribution for ΔN_eff to extend to lower values, and we find ΔN_eff ≳ 0.23 at the 68% credible level and no bound at the 95% credible level. We thus conclude that the dw-dr model prefers large ΔN_eff values in the vicinity of existing bounds. This means that the most promising parameter regions, i.e., regions that are not yet ruled by ΔN_eff but still manage to fit the NANOGrav signal, point to ΔN_eff values within the reach of upcoming experiments, including CMB-S4 (Abazajian et al. 2022), which promises to be sensitive to ΔN_eff values as small as ΔN_eff ≃ 0.06 at the 95% confidence level.

For the dw-sm model, we find α_* ∈ [0.080, 0.19] at the 68% credible level and α_* ∈ [0.053, 0.35] at the 95% credible level. With the inclusion of the SMBHB contribution, smaller values of α_* become allowed, and we find α_* ∈ [0.079, 0.22] at the 68% credible level and α_* ∈ [0.047, 0.61] at the 95% credible level. We also notice a partial degeneracy between T_* and α_*, due to low frequency bins (which are the ones contributing the most to the evidence for a GWB) being fit by the low-frequency tail of the domain wall spectrum. Therefore, a shift to higher transition temperatures (and therefore to higher peak frequencies) can be partially reabsorbed with a shift to higher α_* values.

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 ${ \mathcal B }=14.8\pm 0.5$ and ${ \mathcal B }=1.62\pm 0.05$, respectively, while the Bayes factors for dw-sm+smbhb and dw-dr+smbhb versus smbhb are ${ \mathcal B }=21.1\pm 0.9$ and ${ \mathcal B }=2.53\pm 0.10$, respectively. In both cases, the extra SMBHB contribution helps to improve the fit of the NG15 data. In the case of dw-dr, we moreover observe the same effect as for pt-sound in Section 5.3: 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.

While DM constitutes roughly 27% (Zyla et al. 2020) of the energy density of the universe, very little is known about its fundamental properties. Consequently, a wide range of DM models remain consistent with 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 (Tremaine & Gunn 1979; Di Paolo et al. 2018; Savchenko & Rudakovskyi 2019; Alvey et al. 2020; Davoudiasl et al. 2021). These ULDM models also generically suppress structure on small scales, allowing them to be potential solutions to the small-scale structure problems of the standard ΛCDM paradigm (Bullock & Boylan-Kolchin 2017). However, too much suppression on large scales would be in conflict with CMB measurements (Hlozek et al. 2015; Hložek et al. 2017, 2018), which sets a lower bound on the DM mass, 10⁻²⁴ eV < m_ϕ . PTAs can probe these miniscule masses, since, as we describe below, the frequency of the ULDM signal, f, is generally proportional to the ULDM mass, 2π f ∼ m_ϕ . Therefore, the sensitivity window of NANOGrav is expected to be 10⁻²³ eV ≲ m_ϕ ≲ 10⁻²⁰ eV.

Other astrophysical constraints, such as measurements of the Lyα forest (Armengaud et al. 2017; Iršič et al. 2017; Kobayashi et al. 2017; Rogers & Peiris 2021), galactic subhalo mass functions (Schutz 2020; Banik et al. 2021; Nadler et al. 2021), and stellar kinematics (Dalal & Kravtsov 2022), can push this bound further up, ranging from 10⁻²¹ to 10⁻¹⁹ eV, and the nonobservation of superradiance at SMBHs (Arvanitaki et al. 2015; Stott 2020; Ünal et al. 2021) pushes the bound up to 10⁻²¹ ⋯ 10⁻¹⁷ eV. PTA searches for ULDM in the 10⁻²³ eV ≲ m_ϕ ≲ 10⁻²⁰ eV window can then be viewed as complementary to these astrophysical searches. Signals in PTAs do not depend on the same astrophysical uncertainties, e.g., modeling the nonlinear small-scale matter power spectrum with analytic or numeric methods (Zhang et al. 2018, 2019) or the evolution of specific density profiles. Moreover, the power of these methods quickly degrades when considering subcomponents of DM, or populations of DM that do not compose the whole DM density. The searches discussed here have a weaker dependence on the subcomponent fraction and are therefore still useful in the hunt for DM.

While there is a large phenomenology of ULDM signals that can be produced in PTAs, the timing residuals for the Ith pulsar, h_I , are deterministic and can be written in the form

$$\begin{eqnarray}\displaystyle \begin{array}{rcl}{h}_{I}(t) & = & \displaystyle \sum _{i}{A}_{{\rm{E}},I}^{i}({{\boldsymbol{x}}}_{{\rm{E}}})\sin (\omega t+{\gamma }_{{\rm{E}}}^{i})\\ & & +\,{A}_{P,I}^{i}({{\boldsymbol{x}}}_{P,I})\sin (\omega t+{\gamma }_{P,I}^{i}),\end{array}\end{eqnarray} \tag{ 59 }$$

where i indexes the DM field polarization, ⁸⁹ ω is the signal frequency, and we have split the signal into two terms: an "Earth term" with amplitude ${A}_{{\rm{E}},I}^{i}({{\boldsymbol{x}}}_{{\rm{E}}})$ and time-independent phase ${\gamma }_{{\rm{E}}}^{i}$, and a "pulsar term" with amplitude ${A}_{P,I}^{i}({{\boldsymbol{x}}}_{P,I})$ and time-independent phase ${\gamma }_{P,I}^{i}$. The phases can be written as

$$\begin{eqnarray}&&{\gamma }_{{\rm{E}}}^{i}={\alpha }^{i}({{\boldsymbol{x}}}_{E}),\end{eqnarray} \tag{ 60 }$$

$$\begin{eqnarray}&&{\gamma }_{P,I}^{i}=\omega | {{\boldsymbol{x}}}_{{\rm{E}}}-{{\boldsymbol{x}}}_{P,I}| +{\alpha }^{i}({{\boldsymbol{x}}}_{P,I}),\end{eqnarray} \tag{ 61 }$$

where αⁱ ( x ) is dependent on the underlying ULDM field phase at point x , and the additional term in ${\gamma }_{P,I}^{i}$ is due to the light-travel time between Earth and the pulsar.

An important scale in understanding these signals is given by the correlation length of the ULDM field, ℓ _c ,

$$\begin{eqnarray}&&{{\ell }}_{c}\simeq \displaystyle \frac{2\pi }{{m}_{\varphi }{v}_{\varphi }}\sim 0.4\,\mathrm{kpc}\left(\displaystyle \frac{{10}^{-22}\,\mathrm{eV}}{{m}_{\varphi }}\right),\end{eqnarray} \tag{ 62 }$$

where v_ϕ ∼ 10⁻³ is the ULDM velocity. If the correlation length is much larger than the distance between Earth and the pulsar, ℓ_c ≫ ∣ x _E − x _P,I ∣, both experience the same DM field. In this "correlated" limit of the signals, the amplitudes and phases can be taken to be position independent:

$$\begin{eqnarray}&&{A}_{{\rm{E}},I}^{i}({{\boldsymbol{x}}}_{{\rm{E}}})\to {A}_{{\rm{E}}}^{i},\quad {A}_{P}^{i}({{\boldsymbol{x}}}_{P,I})\to {A}_{P,I}^{i},\quad {\alpha }^{i}({\boldsymbol{x}})\to {\alpha }^{i}.\end{eqnarray} \tag{ 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 , instead of one for each pulsar, to account for fluctuations in the ULDM field. However, since m_ϕ ∣ x _E − x _P,I ∣ ≫ 1, γ_P,I and γ_E must still be taken to be independent. In the uncorrelated limit, ℓ_c ≪ ∣ x _E − x _P,I ∣, the DM field is no longer correlated, and its amplitude fluctuations in different patches need to be accounted for.

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 Equation (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:

• 1.  
  Metric fluctuations (Khmelnitsky & Rubakov 2014; Porayko & Postnov 2014; Porayko et al. 2018; Kato & Soda 2020; Nomura et al. 2020; Unal et al. 2022; Wu et al. 2022): Oscillations in the ULDM field generate fluctuations in the local stress-energy tensor, T_{μ ν} , which are independent of any direct couplings the ULDM may have with SM fields. These fluctuations in T_{μ ν} generate fluctuations in metric perturbations by Einstein's equations. These metric perturbations then affect the photon geodesic on its path from the pulsar and generate a timing residual with pulsar and Earth amplitudes given in Table 1. While the amplitudes shown in Table 1 are 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. 2020), as off-diagonal components of the metric fluctuations can generate nontrivial angular correlations between signals seen in different pulsars. In our analysis, we ignore these spatial correlations and set approximate bounds by treating the vector as three scalar components, such that the vector signal is three times larger than the scalar one (Nomura et al. 2020).
• 2.  
  Doppler–U(1) forces (Graham et al. 2016; Xue et al. 2022): A background of ultralight vector DM can create dark "electric" fields, which can accelerate pulsars, or Earth, and Doppler-shift light arrival times. For example, if the DM is the gauge field of a local baryon symmetry, U(1)_B , then it couples to the Earth and pulsar baryon current density. This coupling generates a force on both Earth and pulsars proportional to their baryon number, Q_B = N_B × q_B , where N_B is the number of baryons and q_B is the gauge charge. This scenario is straightforwardly generalized to other U(1) models, e.g., models where the difference between baryon and lepton numbers, B − L, is a local symmetry. These forces cause periodic displacements between Earth and pulsars and generate timing residuals with amplitudes given in Table 1. These forces also exist for scalar DM coupling to SM operators (e.g., $\varphi \,\bar{n}n$, where ϕ is the DM and n is the neutron field). However, in the scalar case, these forces originate from the field gradient ∇ϕ (Graham et al. 2016), which causes them to be velocity suppressed compared to the vector DM scenario.

Table 1. Summary of the Different Effects Induced by ULDM That Generate a Deterministic Signal of the Form Given in Equation (59)

Effect	${A}_{E}^{i}({\boldsymbol{x}})$	${A}_{P,I}^{i}({\boldsymbol{x}})$	ω = 2π f	Spin (Sϕ )
Metric fluctuations	$(2{S}_{\varphi }+1)\tfrac{\pi G{\rho }_{\varphi }}{2{m}_{\varphi }^{3}}{\hat{\varphi }}^{2}({\boldsymbol{x}})$	$(2{S}_{\varphi }+1)\tfrac{\pi G{\rho }_{\varphi }}{2{m}_{\varphi }^{3}}{\hat{\varphi }}^{2}({\boldsymbol{x}})$	2mϕ	0, 1
Doppler–U(1) forces	${g}_{{\rm{E}}}\tfrac{{Q}_{{\rm{E}}}}{{M}_{{\rm{E}}}}\tfrac{\sqrt{2{\rho }_{\varphi }}}{\sqrt{3}{m}_{\varphi }^{2}}\left({\hat{{\boldsymbol{n}}}}_{I}\cdot {{\boldsymbol{\epsilon }}}_{i}\right){\hat{\varphi }}_{i}({\boldsymbol{x}})$	${g}_{P}\tfrac{{Q}_{P}}{{M}_{P}}\tfrac{\sqrt{2{\rho }_{\varphi }}}{\sqrt{3}{m}_{\varphi }^{2}}\left({\hat{{\boldsymbol{n}}}}_{I}\cdot {{\boldsymbol{\epsilon }}}_{i}\right){\hat{\varphi }}_{i}({\boldsymbol{x}})$	mϕ	1
Pulsar spin fluctuations	0	${y}_{P}{d}_{P}\tfrac{\sqrt{2{\rho }_{\varphi }}}{{m}_{\varphi }^{2}{\rm{\Lambda }}}\hat{\varphi }({\boldsymbol{x}})$	mϕ	0
Reference clock shift	${y}_{{\rm{E}}}{d}_{{\rm{E}}}\tfrac{\sqrt{2{\rho }_{\varphi }}}{{m}_{\varphi }^{2}{\rm{\Lambda }}}\hat{\varphi }({\boldsymbol{x}})$	0	mϕ	0

Note. "E" and "P" subscripts denote Earth and pulsar term contributions, respectively. Aⁱ is the signal amplitude for the ith DM polarization, ω is the signal frequency, and the "Spin" column refers to whether the effect occurs for scalar (spin-0) or vector (spin-1) ULDM candidates. ρ_ϕ is the local DM density, taken to be 0.4 GeV cm⁻³, m_ϕ is the ULDM mass, g are gauge couplings, Q is the charge under the corresponding gauge group, y parameterizes the sensitivity to different fundamental constant fluctuations, d are defined in Equation (64), and ${\rm{\Lambda }}={M}_{\mathrm{Pl}}/\sqrt{4\pi }$. $\hat{\varphi }({\boldsymbol{x}})$ is a random variable representing the fluctuations in the ULDM field in different correlation patches; the correlated limit corresponds to $\hat{\varphi }({\boldsymbol{x}})\to \hat{\varphi }$. n_I is a vector pointing from Earth to the Ith pulsar, and _i are the ULDM polarization vectors.

Download table as:  ASCIITypeset image

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,

$$\begin{eqnarray}&&\displaystyle \begin{array}{l}{ \mathcal L }\supset \displaystyle \frac{\varphi }{{\rm{\Lambda }}}\left[\displaystyle \frac{{d}_{\gamma }}{4{e}^{2}}{F}_{\mu \nu }{F}^{\mu \nu }+\displaystyle \frac{{d}_{g}{\beta }_{3}}{2{g}_{3}}{G}_{\mu \nu }^{A}{G}_{A}^{\mu \nu }\right.\\ \quad -\,\left.\displaystyle \sum _{f=e,\mu }{d}_{f}{m}_{f}\bar{f}f-\displaystyle \sum _{q=u,d}({d}_{q}+{\gamma }_{q}{d}_{g}){m}_{q}\bar{q}q\right],\end{array}\end{eqnarray} \tag{ 64 }$$

where F_{μ ν} and G_{μ ν} are the photon and gluon field strengths, respectively, ${\rm{\Lambda }}={M}_{\mathrm{Pl}}/\sqrt{4\pi }$, β₃ is the QCD beta function, γ_q are the light quark anomalous dimensions, m_f are the fermion masses, and the d values are dimensionless couplings. These couplings induce periodic oscillations in the values of fundamental constants, i.e., particle masses and couplings, which can affect timing residuals in two ways:

• 1.  
  Pulsar spin fluctuations (Kaplan et al. 2022): 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 y parameters are pulsar specific and denote the pulsar sensitivity to a specific coupling, e.g., y_e parameterizes the pulsar sensitivity to changes in the electron mass. ${y}_{\hat{m}}$ parameterizes the pulsar sensitivity to the mass-weighted combination of quark mass couplings, ${d}_{\hat{m}}=({m}_{u}{d}_{u}+{m}_{d}{d}_{d})/({m}_{u}+{m}_{d})$. To be explicit, following Kaplan et al. (2022), we employ in our analysis

  $$\begin{eqnarray}\begin{array}{rcl}{y}_{g} & = & -5,\qquad \qquad {y}_{\hat{m}}=-2.4\,\times \,{10}^{-1},\\ {y}_{\mu } & = & 2\times {10}^{-3},\qquad {y}_{e}=1.7\,\times \,{10}^{-5},\end{array}\end{eqnarray} \tag{ 65 }$$

  for each pulsar, which assumes the simplest possible model for the pulsars' moment of inertia.
• 2.  
  Reference clock shifts (Graham et al. 2016; Kaplan et al. 2022): PTA timing residuals are measured with respect to a collection of mostly cesium (McCarthy 2009) atomic clocks. Therefore, changes to these atomic clock frequencies can appear in PTA data as apparent shifts in timing residuals. These shifts are perfectly correlated among pulsars and have the amplitude reported in Table 1. Similar to the pulsar spin fluctuation signal, the y parameters denote the atomic clock sensitivity to specific couplings. To be explicit, following Kaplan et al. (2022), we take these parameters to be

  $$\begin{eqnarray}\begin{array}{rcl}{y}_{g} & = & 1\qquad \qquad \qquad {y}_{\hat{m}}=0.158\\ {y}_{e} & = & 2\quad \qquad \qquad \quad {y}_{\gamma }=4.83,\end{array}\end{eqnarray} \tag{ 66 }$$

  assuming cesium clocks to be used as a reference.

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: ${\mathrm{log}}_{10}{A}_{\mathrm{GWB}}\,\in [-18,-11]$ and γ_GWB ∈ [0, 7].

Over most of the ULDM mass range, 10⁻²⁴ eV ≲ m_ϕ ≲ 10⁻²⁰ eV, we find no significant evidence for any of the ULDM signals described in the previous section (with ${ \mathcal B }\sim 1$ in favor of models including a ULDM signal on top of the GWB). However, for ULDM models that give rise to an "Earth term" signal, we find mild evidence for an ULDM signal with frequency f ≃ 4 nHz. Specifically, restricting the prior range to f ∈ [3.05, 6.09] nHz, we find a Bayes factor of ${ \mathcal B }\sim 2$ (${ \mathcal B }\sim 1.5$) in favor of the model including the ULDM signal in the correlated (uncorrelated) regime. This result is expected given the excess of monopolar-correlated power observed in the second frequency bin of the common red-noise process (see the discussion in NG15gwb for more details). Unfortunately, an ULDM interpretation of this monopolar signal is difficult, since the corresponding ULDM masses, m_ϕ ∼ 2 × 10⁻²³ eV, needed to explain this excess are in tension with other astrophysical bounds: Lyα forest (Armengaud et al. 2017; Iršič et al. 2017; Kobayashi et al. 2017; Rogers & Peiris 2021), galactic subhalo mass functions (Schutz 2020; Banik et al. 2021; Nadler et al. 2021), and stellar kinematics (Dalal & Kravtsov 2022).

Without convincing evidence for a signal, we compute constraints on the ULDM model parameters, shown in Figures 13–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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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. 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⁻²³ eV, where we nearly constrain ULDM to be a subcomponent of the total DM abundance. While we show constraints down to m_ϕ = 10⁻²⁴ eV, it is easy to extrapolate them to lower masses, where we expect them to remain flat down to m ∼ 10⁻²⁶ eV, where 1/m_ϕ becomes of the same order of the interpulsar separation and the ULDM signal is additionally suppressed (Khmelnitsky & Rubakov 2014; Unal et al. 2022). While future PTA analyses will be able to improve on these constraints, constraining ULDM with m_ϕ ≳ 10⁻²³ eV to have a local abundance smaller than ρ_ϕ = 0.4 GeV cm⁻³ will be challenging, given that the constraints scale as ${\rho }_{\varphi }^{\mathrm{lim}}\propto {m}_{\varphi }^{3}$ for m_ϕ ≳ 1/T_obs.

In Figure 14, we show the constraints for all the d parameters describing the scalar ULDM Lagrangian in Equation (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. (2022) and competitive with laboratory constraints (Boulder Atomic Clock Optical Network et al. 2005; Hees et al. 2016; Bergé et al. 2018; Kumar Poddar et al. 2019; Dror et al. 2020). The strongest constraints, relative to laboratory bounds, are for ULDM models coupled to the electron, d_e , and muon, d_μ , mass terms. Indeed, relative shifts of the energy levels utilized in atomic clock experiments are insensitive to the d_e coupling, since atomic energy levels in different atoms scale identically with electron mass, leading to no relative energy level shifts (Van Tilburg et al. 2015; Kaplan et al. 2022). 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_μ is simply because there are not a large number of muons to study on Earth, whereas pulsars host a large number of them. As for the gravitational signal, we can extrapolate the constraints to lower masses, where we expect them to scale as d ∝ 1/m_ϕ down to m_ϕ ∼ 10⁻²⁶ eV, where 1/m_ϕ becomes comparable to the interpulsar separation.

In Figure 15, we show constraints on the gauge coupling of models where the ULDM is the gauge boson of either U(1)_B or U(1)_B−L . Our constraints are roughly consistent with those published by the PPTA Collaboration (Porayko et al. 2018). 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 Figures 14 and 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\sqrt{{\rho }_{\varphi }}$ for the scalar case and $Q\sqrt{{\rho }_{\varphi }}$ for the vector case. Therefore, if ULDM is only a faction f_ϕ of the total DM abundance, the constraints are weakened by a factor $\sqrt{{f}_{\varphi }}$.

Lastly, we comment on the prior choice for $\hat{\varphi }$. Previous studies (Porayko & Postnov 2014; Kato & Soda 2020) assumed that $\hat{\varphi }$ is not a random variable in the problem and placed constraints assuming that this parameter is simply 1. However, as pointed out by Centers et al. (2021), this assumption is not appropriate. Since the observation time of PTAs is much smaller than the coherence time of the DM, $\tau \sim {({m}_{\varphi }{v}^{2})}^{-1}$, only one instance of the field is being sampled within each correlation patch. The DM density in any correlation patch is then a random number, which follows a Rayleigh distribution with mean ρ_ϕ , and the priors should be chosen to reflect this. We also note that while Porayko et al. (2018) find the distribution of $\hat{\varphi }$ through numerical simulation, these results are consistent with the analytic predictions by Foster et al. (2018).

In the ΛCDM model, the structures we observe in the universe are seeded by primordial curvature fluctuations generated during inflation and then imprinted onto the DM density field. CMB observations indicate that these fluctuations have a nearly scale-invariant power spectrum on large scales (i.e., for comoving wavenumbers k ≃ Mpc⁻¹). However, on smaller scales, various theories of DM leave unique fingerprints on primordial perturbations or their evolution, resulting in different predictions for the amount of subgalactic DM substructures. Consequently, measuring local DM substructures could be crucial in determining the correct model of DM.

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 Section 5.2). Several studies investigated the possibility of identifying a galactic PBH population by analyzing the Doppler and Shapiro signals they can leave in PTA data (Seto & Cooray 2007; Siegel et al. 2007; Dror et al. 2019; Ramani et al. 2020; Lee et al. 2021a, 2021b). In this analysis, we will closely follow the method outlined by Lee et al. (2021b) to constrain the local PBH abundance. ⁹⁰

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, τ, the signal can be further classified as dynamic (static) if τ is much smaller (larger) than the observation time, T_obs. In the static regime, the leading-order term of the Doppler signal that is not degenerate with the timing model is given by (Ramani et al. 2020; Lee et al. 2021a, 2021b)

$$\begin{eqnarray}&&{h}_{D,\mathrm{sta}}(t)=\displaystyle \frac{{A}_{D,\mathrm{sta}}}{{\mathrm{yr}}^{2}}\,{t}^{3},\end{eqnarray} \tag{ 67 }$$

where A_D,sta is a dimensionless parameter that can be related to the kinematic parameters of the transiting event (see Appendix C.6 for more details). In the dynamic limit, and assuming that the signal is dominated by the closest transiting PBH, we get

$$\begin{eqnarray}&&{h}_{D,\mathrm{dyn}}(t)={A}_{D,\mathrm{dyn}}(t-{t}_{0}){\rm{\Theta }}(t-{t}_{0}),\end{eqnarray} \tag{ 68 }$$

where Θ is the Heaviside step function, A_D,dyn is a dimensionless amplitude that can be related to kinematic parameters of the transiting event, and t₀ is the time of closest approach (see Appendix C.6 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

$$\begin{eqnarray}&&{h}_{S,\mathrm{sta}}(t)=\displaystyle \frac{{A}_{S,\mathrm{sta}}}{{\mathrm{yr}}^{2}}\,{t}^{3},\end{eqnarray} \tag{ 69 }$$

where, as for the Doppler case, A_S,sta is a dimensionless parameter that can be related to the kinematic parameters of the transiting event (see Appendix C.6 for more details). In the dynamic limit, there is no simple parameterization 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_PBH ≡ Ω_PBH/Ω_DM, as a function of the PBH mass, M: p(f_PBH∣ δ t , M). We do this as follows:

• 1.  
  For each given value of f_PBH and M, we use the Monte Carlo code developed by Lee et al. (2021a) to generate a PBH population surrounding each of the pulsars in our array. From this distribution, we derive the amplitude of the static Doppler and Shapiro signals generated by the entire PBH population and the amplitude for the dynamic Doppler signal generated by the closest transiting PBH. Finally, we repeat this procedure for 2.5 × 10³ realizations to obtain the conditional probability distributions $p({\mathrm{log}}_{10}{A}_{I}| {f}_{\mathrm{PBH}})$, where I indexes pulsars in the array and A refers to any of the PBH signal amplitudes introduced in Equations (67), (68), and (69). In Figure 16 we report some of the distributions derived in this way. ⁹¹
• 2.  
  One at a time, we include the PBH signals given in Equations (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({\mathrm{log}}_{10}{\boldsymbol{A}}| {\boldsymbol{\delta }}{\boldsymbol{t}})$. Since the PBH signal in different pulsars is assumed to be independent, these distributions can be factorized as

  $$\begin{eqnarray}&&p({\mathrm{log}}_{10}{\boldsymbol{A}}| {\boldsymbol{\delta }}{\boldsymbol{t}})=\displaystyle \prod _{I=1}^{{N}_{P}}p({\mathrm{log}}_{10}{A}_{I}| {\boldsymbol{\delta }}{\boldsymbol{t}}).\end{eqnarray} \tag{ 70 }$$

  Some of the $p({\mathrm{log}}_{10}{A}_{I}| {\boldsymbol{\delta }}{\boldsymbol{t}})$ are reported in Figure 16.
• 3.  
  Finally, we can write

  $$\begin{eqnarray}&&\begin{array}{l}p({f}_{\mathrm{PBH}}| {\boldsymbol{\delta }}{\boldsymbol{t}})=\displaystyle \prod _{I=1}^{{N}_{P}}\int d{\mathrm{log}}_{10}{A}_{I}\ p({f}_{\mathrm{PBH}}| {\mathrm{log}}_{10}{A}_{I})p({\mathrm{log}}_{10}{A}_{I}| {\boldsymbol{\delta }}{\boldsymbol{t}})\\ \ \propto \,\displaystyle \prod _{I=1}^{{N}_{P}}\int d{\mathrm{log}}_{10}{A}_{I}\ p({\mathrm{log}}_{10}{A}_{I}| {f}_{\mathrm{PBH}})p({\mathrm{log}}_{10}{A}_{I}| {\boldsymbol{\delta }}{\boldsymbol{t}}),\end{array}\end{eqnarray} \tag{ 71 }$$

  where, in the second step, we used Bayes's theorem and assumed uniform priors on ${\mathrm{log}}_{10}{A}_{I}$ and f_PBH. More details on each of these three steps can be found in Appendix C.6 or in Lee et al. (2021b).

Zoom In Zoom Out Reset image size

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

$$\begin{eqnarray}&&{V}_{\mathrm{fifth}}(r)=-\tilde{\alpha }\displaystyle \frac{{{GMM}}_{P}}{r}{e}^{-r/\lambda },\end{eqnarray} \tag{ 72 }$$

where M and M_P are the masses of the DM and pulsar, respectively, λ is the range of the interaction, and $\tilde{\alpha }$ is the effective DM−baryon coupling, normalized against the gravitational coupling (also known as the Yukawa parameter). Here the DM can be either a particle or a macroscopic object such as a nugget of asymmetric DM (Detmold et al. 2014; Wise & Zhang 2014; Hardy et al. 2015; Krnjaic & Sigurdson 2015; Gresham et al. 2017, 2018). These Yukawa interactions can arise from an effective Lagrangian of the form ${ \mathcal L }\supset {g}_{X}\varphi \bar{X}X+{g}_{n}\varphi \bar{n}n$, where X and n are the effective DM and neutron fields, respectively, and ϕ is a massive (but potentially light) scalar or vector field. The effective coupling is related to the coupling constants by $\tilde{\alpha }\approx \tfrac{{g}_{n}{g}_{X}}{4\pi {{Gm}}_{X}{m}_{n}}$, where m_n is the neutron mass. These interactions are constrained to be weaker than gravity for the mass range M ≲ 100 GeV by the CMB, Lyα forest (Xu et al. 2018), and direct detection experiments such as X-ray Quantum Calorimeter (XQC; Mahdawi & Farrar 2018) and Cryogenic Rare Event Search with Superconducting Thermometers (CRESST; Angloher et al. 2017) (for a review on these constraints see Xu & Farrar 2020). However, stronger-than-gravity fifth forces are allowed if M ≫ 100 GeV, even when X accounts for the entirety of the DM population.

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 C.6). Therefore, as shown by Gresham et al. (2023), following a procedure similar to the one used to constrain the abundance of PBHs, we can constrain the value of the Yukawa parameter, $\tilde{\alpha }$ . Specifically, for each given value of $\tilde{\alpha }$ and M, we use the Monte Carlo code developed by Lee et al. (2021a) to generate a population of DM substructure surrounding each of the pulsars in our array. From this distribution, we derive the amplitude of the static Doppler signal generated by the closest transiting substructure by considering the acceleration induced by both the gravitational and Yukawa interaction. By repeating this procedure for multiple populations of DM substructure, we derive the distribution $p({\mathrm{log}}_{10}{A}_{I}| \tilde{\alpha })$. By plugging this quantity into an expression similar to the one given in Equation (71), we can derive $p(\tilde{\alpha }| {\boldsymbol{\delta }}{\boldsymbol{t}})$ and use this quantity to constrain $\tilde{\alpha }$.

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: ${\mathrm{log}}_{10}{A}_{\mathrm{GWB}}\in [-11,-18]$ and γ_GWB ∈ [0, 7]. We find no statistically significant evidence for any of the PBH signals described in the previous section. Therefore, we proceed to set constraints on the local PBH abundance. The prior distributions used for the PBH signal parameters are reported in Table 3.

The 95% upper limits on f_PBH derived from the static Doppler and Shapiro signals are reported in Figure 17. The dynamic Doppler signal is too weak to produce any detectable signal for any of the f_PBH values considered. These are the first constraints on f_PBH derived using real PTA data. As expected, our constraints are much weaker compared to the projections that were derived by Lee et al. (2021b) using mock data and including only white noise. Indeed, as already discussed by Lee et al. (2021b), the presence of a common red-noise process significantly reduces the sensitivity to PBH signals.

Zoom In Zoom Out Reset image size

Finally, in Figure 18 we show the constraints on $\tilde{\alpha }$ set by NG15 data. These constraints are compared with several other constraints that can be placed on $\tilde{\alpha }$. Specifically, in teal we show weak equivalence principle (WEP) constraints (Wagner et al. 2012; Shao et al. 2018; Sun et al. 2019; properly rescaled to take into account the finite range of the interaction, Gresham et al. 2023) derived by considering differential acceleration of baryonic test bodies toward the galactic center. In blue we report constraints from neutron star (NS) heating (assuming additional short-range DM−baryon interaction) induced by DM capture (Gresham et al. 2023), derived from the coldest known NS to date—PSR J2144−3933 (Guillot et al. 2019). In gray we report the indirect constraints that can be derived by combining the fifth-force constraints on baryon−baryon interactions (Bergé et al. 2018; Fayet 2019) and Bullet Cluster constraints on DM−DM interactions (Spergel & Steinhardt 2000; Kahlhoefer et al. 2014; see also Coskuner et al. 2019; Gresham et al. 2023). We find that the NG15 constraints are competitive with WEP and NS kinetic heating, especially in the M ≳ 10⁻¹² M_⊙ regime. However, the combined constraint from the Bullet Cluster and MICROSCOPE dominates over all other constraints in the entire mass range of interest. We note, however, that the combined constraint is an indirect probe of the Yukawa parameter and depends on the specific form of the light mediator coupling to DM and baryons. Specifically, in deriving the combined constraints, we have assumed that the Yukawa DM−baryon interaction arises from a Lagrangian of the form ${ \mathcal L }\supset {g}_{X}\varphi \bar{X}X+{g}_{n}\varphi \bar{n}n$. In addition, if only a subcomponent of DM interacts through the long-range fifth force, then the Bullet Cluster constraint will quickly lift, while the other constraints only deteriorate linearly with the DM fraction.

Zoom In Zoom Out Reset image size