Constraining gravitational wave propagation using pulsar timing array correlations

TL;DR

This paper uses pulsar timing array data to test for subluminal gravitational wave modes, constraining the graviton mass and finding no evidence for certain GW correlation types.

Contribution

It introduces a method to constrain graviton mass using PTA correlations and applies it to real data, providing new upper bounds.

Findings

No evidence for tensor- or vector-induced GW correlations.

Vector correlations are disfavored by the data.

Upper bound on graviton mass is approximately 10^{-22} eV.

Abstract

Pulsar timing arrays (PTA) are a promising probe to the cosmologically novel nanohertz gravitational wave (GW) regime through the stochastic GW background. In this work, we consider subluminal GW modes as a possible source of correlations in a PTA, utilizing the public code PTAfast and the 12.5 years correlations data by NANOGrav, which we hypothesize are sourced by GWs. Our results show no evidence in support of tensor- or vector-induced GW correlations in the data, and that vector correlations are disfavored. This places an upper bound to the graviton mass, $m_{\rm g} \lesssim 10^{-22}$ eV, characteristic of the PTA GW energy scale.

Tables

Table 1. Table 1: Marginalized statistics for subluminal GW cross correlations constrained by PTA Arzoumanian  et al. ( 2020 ) . Results for the HD and the GW-like monopole (GW mon) are presented for comparison. The performance statistics (chi-squared, AIC, and BIC Trotta ( 2008 ); Liddle ( 2007 ) ) are relative to the GW monopole, i.e., a positive value means statistical preference over the GW monopole.

$\frac{1}{1}$ mode	111 $v$ 111	$A^2$ [$\times {10}^{-30}$]	$\Delta {\overline{\chi }}^2$	$\Delta$AIC	$\Delta$BIC
$\frac{1}{1}$ Tensor	${0.45}_{-0.42}^{+0.17}$	${7.4}_{-2.4}^{+2.1}$	$0.06$	$-1.01$	$-1.79$
$\frac{1}{1}$ Vector	${0.39}_{-0.37}^{+0.14}$	${9.2}_{-6.0}^{+2.4}$	$-0.11$	$-3.76$	$-4.53$
$\frac{1}{1}$ HD	$v=1$	$3.9\pm 1.1$	$-0.10$	$-1.66$	$-1.66$
$\frac{1}{1}$ GW mon	——-	$2.0\pm 0.5$	$0$	$0$	$0$

Equations

$$1-v=1-\sqrt{1-{\left(\frac{m_{\mathrm{g}}}{\omega }\right)}^2}.$$

$${\left(\frac{m_{\mathrm{g}}}{\omega }\right)}^2\lesssim 2\times 10^{-15}.$$

$$m_{\mathrm{g}}\lesssim \left(\frac{f}{100\mathrm{hertz}}\right)\times 10^{-19}\mathrm{eV}.$$

$${\gamma }_{ab}\left(\zeta \right)=\sum _l\frac{2l+1}{4\pi }{C}_l{P}_l\left(\cos \zeta \right),$$

$$C_l=\frac{J_l\left(f{D}_a\right)J_l^{\ast }\left(f{D}_b\right)}{\sqrt{\pi }},$$

$$J_l\left(y\right)=\sqrt{2}\pi i^l\sqrt{\frac{(l+2)!}{(l-2)!}}{\int }_0^{2\pi yv}\frac{dx}{v}{e}^{ix/v}\frac{j_l\left(x\right)}{x^2},$$

$$J_l\left(y\right)=2\sqrt{2}\pi i^l\sqrt{l\left(l+1\right)}{\int }_0^{2\pi yv}\frac{dx}{v}{e}^{ix/v}\frac{d}{dx}\left(\frac{j_l(x)}{x}\right).$$

$${\gamma }_{aa}={\int }_0^{\pi }\frac{d\theta }{\sqrt{4\pi }}\left(\frac{2\pi {\sin }^5\theta {\sin }^2(\pi fD(1+v\cos \theta ))}{(1+v\cos \theta {)}^2}\right),$$

$${\gamma }_{aa}={\int }_0^{\pi }\frac{d\theta }{\sqrt{4\pi }}\left(\frac{8\pi {\sin }^3\theta {\cos }^2\theta {\sin }^2(\pi fD(1+v\cos \theta ))}{(1+v\cos \theta {)}^2}\right).$$

$$\log \mathcal{L}=-\frac{1}{2}\sum _{\zeta }{\left(\frac{A^2{\Gamma }_{ab}\left(\zeta \right)-\mathrm{CCP}\left(\zeta \right)}{\Delta \mathrm{CCP}\left(\zeta \right)}\right)}^2,$$

$${\Gamma }_{ab}\left(\zeta \right)=\frac{3}{2}\left[\frac{1}{3}+\left(\frac{1-\cos \zeta }{2}\right)\left(\ln \left(\frac{1-\cos \zeta }{2}\right)-\frac{1}{6}\right)\right]$$

$${\Gamma }_{ab}\left(\zeta \right)=\frac{{\delta }_{ab}}{2}+\frac{1}{2},$$

$$C_2(\mathrm{quadrupole})\gtrsim C_1(\mathrm{dipole}),C_3(\mathrm{octupole})>\cdots ,$$

Full text

I Introduction

The next gravitational wave (GW) astronomy breakthrough is the detection of the stochastic gravitational wave background (SGWB) – a galaxy size superposition of GWs from various cosmological sources – that correlates astronomical observations Burke-Spolaor et al. (2019); Pol et al. (2021). Along this theme, it was advocated decades ago that the precise timing of the emission of radio pulses by nearby millisecond pulsars is equivalent to a galactic size GW detector Detweiler (1979). This proposal reached fruition in so–called ‘pulsar timing array’ (PTA) missions that aim to capture the wavelike features of the SGWB, and work as probes of the cosmologically novel nanohertz GW regime Romano and Cornish (2017); Romano (2019).

As the CMB is host to physics during last scattering, the SGWB is host to a variety of sources tied to the early Universe, such as supermassive black hole binary coalescences Phinney (2001); Wyithe and Loeb (2003); Shannon et al. (2015), and more exotic phenomena like cosmic strings Ellis and Lewicki (2021); Buchmuller et al. (2021); Hindmarsh and Kume (2022), phase transitions Arzoumanian et al. (2021a); Xue et al. (2021), and primordial black holes Nakama et al. (2017); Chen et al. (2020), among others Caprini and Figueroa (2018); Sharma (2022). PTA astronomy caters to such exciting physics that come with the SGWB.

The present PTAs, spearheaded by the North American Nanohertz Observatory for Gravitational Waves (NANOGrav) Arzoumanian et al. (2020), European PTA Chen et al. (2021a), and Parkes PTA Goncharov et al. (2021), working together under the banner International PTA (IPTA) Hobbs et al. (2010); Chen et al. (2021b), finds strong evidence for a stochastic common process across millisecond pulsars. However, the present data have yet to reach maturity to be able to tell with confidence whether this common process is due to the SGWB. This is where PTA cross correlations play a significant role, since the SGWB gives a distinct quadrupolar spatial correlation Hellings and Downs (1983); Ng (2022); Allen (2023); Allen and Romano (2022) that cannot be associated with systematic effects such as clock error and ephemeris uncertainty Tiburzi et al. (2016); Roebber (2019); Vallisneri et al. (2020). As the PTAs continue to monitor more pulsars, SGWB correlations should soon manifest in the data Pol et al. (2022), marking PTA science precision era.

Adding to the science prospects of PTAs, in this work, we take a conservative stand on SGWB correlations and consider subluminal GW modes as a source of the observed cross correlations in a PTA Qin et al. (2021); Bernardo and Ng (2023, 2022). We use the public code PTAfast Bernardo and Ng (2022) and the correlations data set by NANOGrav to constrain the GW speed, and henceforth give its take on the graviton mass.

The rest of this paper proceeds as follows. We first address the feasibility of subluminal GWs given the known tight constraint on GW propagation at sub-kilohertz frequencies Abbott et al. (2017a) (Section II). We then briefly go over the power spectrum method that is implemented to calculate the SGWB spatial correlation signals (Section III.1), and describe the data (Section III.2) to be considered in upcoming data analysis (Section IV.1). We present our main results in two parts, beginning with only subluminal GW modes (Section IV.2), and following it up with an assessment of its significance (Section IV.3). We conclude with a summary of our results and draw future directions on subluminal GW propagation (Section V).

II Subluminal GWs

Before we push through, we discuss the feasibility of subluminal tensors in the PTA/nanohertz GW band, particularly in light of the GW astronomy constraint $|1-v|\lesssim 10^{-15}$ at $f\sim 100$ hertz Abbott et al. (2017a, b).

We start by writing down the dispersion, $\omega^{2}=k^{2}+m_{\rm g}^{2}$, expressed in terms of the velocity, $v=d\omega/dk$, wave number $k$, and particle (graviton) mass $m_{\rm g}$, as

[TABLE]

Thus, providing $1-v\sim 10^{-15}$, and taking physical considerations into account, $v<1$ and $m_{\rm g}\leq\omega$, we find

[TABLE]

At $f\sim 100$ hertz (or $\omega\sim 2\pi 100$ hertz), we see that the LIGO/Virgo band, in this crude calculation, translates to a graviton mass upper bound of $m_{\rm g}\lesssim 10^{-19}$ eV. In general, the GW speed constraint converts to

[TABLE]

This expression shows that in PTAs (at much lower frequencies, $f\sim 1$ per year) the graviton mass bounds are always going to be acceptable within the LIGO/Virgo GW speed constraint. In other words, subluminal GW propagation is allowed by the present constraints de Rham and Melville (2018).

Keeping this in mind, we proceed to the SGWB correlations induced by subluminal GWs in PTA.

III PTA Correlations

We briefly describe the PTA correlations in terms of the power spectrum, and the data which we will use. But, first we discuss the feasibility of subluminal tensors in the nanohertz GW band.

III.1 SGWB correlations

The SGWB spatially correlates the time of arrival perturbations of the radio pulses emitted by nearby millisecond pulsars. This angular correlation, $\gamma_{ab}\left(\zeta\right)$, can be written in terms of the power spectrum, $C_{l}$’s, as

[TABLE]

where $\zeta$ is the angular separation of pulsars $a$ and $b$ in the sky.

For this work, we consider only tensor and vector polarizations. Subluminal scalar GWs were discussed in Bernardo and Ng (2022). In general, the power spectrum multipoles can be written as

[TABLE]

where for tensor modes the $J_{l}(y)$’s are

[TABLE]

whereas for vector modes the $J_{l}(y)$’s are instead

[TABLE]

Above, $j_{l}(x)$’s are the spherical Bessel functions, $f$ is the GW frequency, $v$ is the GW speed, and $D_{i}$’s are the distances to the pulsars.

The overlap reduction function (ORF), $\Gamma_{ab}\left(\zeta\right)$, or the PTA’s spatial correlation observable, is obtained simply by normalizing $\gamma_{ab}\left(\zeta\right)$ such that $\Gamma_{ab}\left(0\right)=0.5$ for tensors with $v=1$. We refer the reader to Bernardo and Ng (2023, 2022) for a detailed description of this procedure, and a derivation of the results put forward in this section.

The autocorrelation, $\gamma_{aa}$, embodying a pulsar’s correlation with itself, is also of relevance in this work. This physical quantity admits analytical expressions, as one dimensional integrals, which we give below for the tensor and vector modes Bernardo and Ng (2023): for tensor polarizations the autocorrelation is

[TABLE]

and for vector polarizations the expression is

[TABLE]

For our purposes, we apply PTAfast Bernardo and Ng (2022) which implements all the intervening calculations, and provides the ORF’s and autocorrelations relevant for comparison with observations.

III.2 Data

The NANOGrav Arzoumanian et al. (2020), among other PTAs Chen et al. (2021a); Goncharov et al. (2021), confirm the existence of a stochastic common spectrum process. This is not sufficient to tell unequivocally whether the SGWB is out there, but it nonetheless is able to give bounds to the SGWB amplitude. In particular, for a spectral density, characteristic of supermassive black hole binaries, the NANOGrav gives the GW characteristic strain, $A_{\rm CP}$, a median $1.92\times 10^{-15}$, and $5-95\%$ quantiles $1.37-2.67\times 10^{-15}$ at $f=1$ yr*-1*.

Spatial correlations, on the other hand, reveal the SGWB, as it perturbs the arrival of the radio pulses from the millisecond pulsars, delaying some while advancing others, depending on when and where they stand on the wave. However, quite similar with the CMB, digging up this GW correlation from the noise requires the precise monitoring of as many as possible pulsars for long periods. For the moment, the PTAs altogether keep an eye on fewer than a hundred millisecond pulsars, which remain unperceptive of the SGWB, but nonetheless present nontrivial features in the correlation.

For this work, we make use of the correlation data points by NANOGrav Arzoumanian et al. (2020) that measures the average cross correlated power (CCP) across the pulsars in the PTA, ${\rm CCP}\left(\zeta\right)\pm\Delta{\rm CCP}\left(\zeta\right)$. In addition, we utilize the GW amplitude of the common process as a data point, $A_{\rm CP}^{2}=A^{2}\Gamma_{aa}$, that we associate with the autocorrelation.

The 12.5 years noise–marginalized cross correlation data by NANOGrav comprises of timing residual observations of 45 millisecond pulsars, which make up 990 pulsar pairs that we distribute evenly over 15 bins in angular space Arzoumanian et al. (2020). The error $\Delta{\rm CCP}$ thus contains all statistical and systematic uncertainty leftovers after marginalizing over the intrinsic and nuisance parameters pertaining to each millisecond pulsar.

IV Results

We present our main results in this section, starting with our data analysis implementation (Section IV.1), and then moving on to the correlations (Sections IV.2 and IV.3).

IV.1 Bayesian analysis

We model the CCP in a PTA as $A^{2}\Gamma_{ab}\left(\zeta\right)$, making up a two dimensional parameter space, $v\times A^{2}$, of the GW speed $v$ and the amplitude squared $A^{2}$. We compare this with the data through the likelihood function $\mathcal{L}$, which we take to be

[TABLE]

where the sum runs over the angular separations of pulsar pairs in the PTA.

We also make use of the amplitude of the common process, which we symmetrize for our purposes to be $A_{\rm CP}^{2}=\left(3.68\pm 1.58\right)\times 10^{-30}$. This conservative value is obtained from the median and larger error bar of the asymmetric bound of the common process Arzoumanian et al. (2020). Admittedly, this is not the most faithful reconstruction from the best estimate, but it encompasses all the relevant values. Another point that we bring up once more, as alluded to previously, is that this value $A_{\rm CP}$ was obtained assuming that the SGWB comes mainly from supermassive black hole binaries, and so carries with it a mild model dependence despite the large uncertainty. We use it anyway for the analysis, but guarantee that the main results of this paper holds even with only the correlations, which dominate the bulk of the data.

Now that the model, data, and the likelihood are laid out, we obtain the posterior, $P(\theta|D)$, through Bayes theorem, $P\left(\theta|D\right)\propto P\left(\theta\right){\cal L}$, by sampling over the parameter space $\theta=\left(v,A^{2}\right)$ via a Markov chain Monte Carlo (MCMC) algorithm Trotta (2008). We consider flat priors $P\left(\theta\right)$ for both the parameters, $v\in[0.01,1.00]$ and $A^{2}\in[0.01,50]\times 10^{-30}$, and as a starting point (each time) take the references $v\in(0.48,0.52)$ and $A^{2}\in(1,2)\times 10^{-30}$. The references give a better guarantee when the run has converged since the initial points are always different. This is in addition to specifying a strict Gelman-Rubin convergence criterion of $R-1=10^{-3}$ to mark the end of the sampling.

We implement these steps using PTAfast Bernardo and Ng (2022) together with the cosmology community code Cobaya Lewis (2019), and analyze the results using GetDist Torrado and Lewis (2021). Our codes and resulting MCMC chains are publicly accessible in GitHub.

IV.2 Subluminal GWs and the graviton

Figure 1 shows the posterior resulting from PTA correlation observations.

The results show that the PTA correlations constrain the GW amplitude and the GW speed, despite the low level of precision of the present correlation measurements. Furthermore, an interesting angle here is that the data supports subluminal GW modes ($v<1$) more than luminal ones ($v=1$) at 95% confidence. We remind that subluminal GW propagation is consistent with the GW speed constraint by LIGO/Virgo Abbott et al. (2017a). Caution must however be taken care of when interpreting the nonrelativistic speed limit ($v\ll 1$) since this regime is sensitive to the pulsar distances. Granted, the large angle correlations that the current PTAs are able to see are very weakly influenced by how far the pulsars are from the observer Ng (2022). However, this picture changes when entering the nonrelativistic limit ($v\ll 1$) where the SGWB power spectrum exhibits modes, beyond the quadrupole, that manifest at large angles Bernardo and Ng (2023, 2022). Nonetheless, the takeaway of this result is clear: that the uncertainty in the present PTA data is able to support subluminal GW propagation.

We add that these hold regardless of the nature, that is, tensor or vector, of the GWs making up the SGWB. Scalar field induced GWs were considered in Bernardo and Ng (2022) where similar conclusions were drawn. At the same time, there are noticeable differences that discriminate between these tensor and vector modes. For one, the GW amplitude is more constrained for tensors than for vectors, i.e., tensors have narrower $A^{2}$ posterior compared to vectors. It is also worth adding that the correlations strongly disfavor luminal vector GW modes, which we understand is a reflection of the divergence of vector correlations at luminal speed Arzoumanian et al. (2021b). Tensor SGWB correlations are in fact the only ones that does not lose predictability at luminal speeds Bernardo and Ng (2023, 2022).

Furthermore, as teased in Section II, GW speed constraints translate to a graviton mass bound. Working with this relation, we find the sampled graviton mass posteriors using the PTA correlations. The results are shown in figure 2.

This makes clear the relation between the GW speed and the graviton mass where now the PTA constraint on the GW speed translates to a graviton mass posterior: $m_{\rm g}\simeq\left(1.3105\times 10^{-22}{\rm\ eV}\right)\times f\left[{\rm yr}^{-1}\right]\times\sqrt{1-v^{2}}$. The peak in the graviton mass at $m_{\rm g}\sim 1.3\times 10^{-22}{\rm\ eV}$ thus relates to $v\sim 0.1$ where we see the GW speed samples cluster around. This similarly explains the squeezed shape of the graviton mass posterior as the slim mass range $m_{\rm g}\sim(1.3-1.31)\times 10^{-22}{\rm\ eV}$ maps with $0<v\lesssim 0.1$ whereas the rest of $m_{\rm g}\lesssim 1.3\times 10^{-22}{\rm\ eV}$ goes to $0.1\lesssim v\leq 1$. It is also worth stressing that this relation is nonlinear, and so the shapes of the posteriors do not transform trivially. This explains why the mass upper limit $m_{\rm g}\sim 1.4\times 10^{-22}$ eV, where the probability is nearly vanishing, do not manifest as an equivalently nearly vanishing GW speed probability at $v\ll 1$. We also note that the smallest $v$ considered is $v=0.01$, which converts to a mass upper bound $m_{\rm g}\simeq 1.31\times 10^{-22}$ eV during the sampling. A related subtle point is that as far as the sampling on $m_{\rm g}$ goes $m_{\rm g}\lesssim 1.3105\times 10^{-22}{\rm\ eV}$, where the upper bound is determined by the GW frequency that PTAs are sensitive to ($f\sim 1{\rm\ yr}^{-1}$). Values that seemingly lie outside this in figure 2 are only an artefact of the kernel density estimator considered to visualize the sampled posteriors Torrado and Lewis (2021).

A subluminal GW speed in the PTA GW band has important implications for fundamental physics. Together with the LIGO/Virgo constraint Abbott et al. (2017b), this gives rise to dispersive nature of GWs, which could be associated with massive gravity or scalar-tensor theories with kinetic couplings, e.g., $S[g,\phi]\sim\left(\partial\phi\right)^{2}R$, in the gravity sector. This revives interest in the couplings that were let go in favor of trivially satisfying the LIGO/Virgo constraint at the subkilohertz band Ezquiaga and Zumalacárregui (2017); Baker et al. (2017). Since nontrivial couplings also come with predictions in other regimes, such as cosmology Kase and Tsujikawa (2019); Ferreira (2019) or strong fields Tattersall et al. (2018); Minamitsuji et al. (2022), this would tie in PTA observations with other experiments to make stronger physics constraints.

IV.3 Subluminal GWs vs HD vs GW monopole

We take a look at subluminal GWs within the broader context of the Hellings-Downs (HD) curve and the GW-like monopole, or for short GW monopole Arzoumanian et al. (2021b). These give rise to spatial correlations with well established analytical forms,

[TABLE]

for the HD correlation, and

[TABLE]

for the GW monopole. The $\delta_{ab}$ in (12) is considered to take care of the small scale power that is dropped when pulling the pulsars to infinity. In Bernardo and Ng (2022), we have proposed a physical mechanism for the GW monopole by means of a scalar field Appleby and Bernardo (2022). Similarly, for the HD, $\Gamma_{aa}=1$ is set in by hand to account for the small scale power Arzoumanian et al. (2021b); Ng (2022).

These traditional correlation signals were given their due consideration by PTAs Arzoumanian et al. (2021b) and so it is important to assess the statistical significance of subluminal GWs with respect to these correlations. We model their cross correlated power in the PTA by multiplying the relevant ORFs, $\Gamma_{ab}\left(\zeta\right)$, by the GW amplitude squared, $A^{2}$. In contrast with subluminal tensors, both HD and GW monopole carry one parameter, $A^{2}$, that needs to be fit with the correlation observations.

We assess the statistical significance of subluminal tensors, HD, and the GW monopole by means of the reduced chi–squared statistic, $\overline{\chi}$, often considered the baseline for quantifying models’ goodness of fit. We also consider the Akaike information criterion (AIC) and the Bayesian information criterion (BIC) Trotta (2008); Liddle (2007), which in addition to providing a measure of the fit, also penalizes unnecessary addition of free parameters. To put briefly, the smaller $\overline{\chi}$, AIC, and BIC are, the better the fit. The relevant statistics of the sampling for the subluminal GW modes as well as for the HD and the GW monopole are summarized in Table 1.

We referred to the GW monopole since this phenomenological model has a competitive S/N in the NANOGrav data Arzoumanian et al. (2021b). More concretely, we expressed the metrics ${\rm X}\left(=\chi^{2},{\rm AIC},{\rm BIC}\right)$ in Table 1 as $\Delta{\rm X}={\rm X}_{\rm GWmon}-{\rm X}_{i}$ such that if ${\rm X}_{i}<{\rm X}_{\rm GWmon}$, or that correlations model $i$ fits better compared with the GW monopole, then $\Delta{\rm X}>0$, and vice versa. In short, a value in the more positive is better. At a glance, Table 1 shows almost all $\Delta\overline{\chi}^{2}$, $\Delta{\rm AIC}$, and $\Delta{\rm BIC}$ to be negative, meaning that none of the tensor and vector correlation models fit the data better than the GW monopole.

The GW monopole correlations represent the current data quite well, in agreement with earlier result Bernardo and Ng (2022). We emphasize that the GW monopole is physically distinct from the systematic monopole as it brings small scale power, and can be provided by a physical mechanism Arzoumanian et al. (2021b); Bernardo and Ng (2022). In short this is a clear cut from the clock errors that were disfavored by the PTAs Arzoumanian et al. (2020); Goncharov et al. (2021).

Letting the tensors move off the light cone improve the likelihood of their correlations in the data; however, this comes at a price since when considering the AIC and BIC, tensor correlations turn out disfavorable over the GW monopole. This means that there is no evidence for tensor correlations, may it be subluminal ones or the HD. The error bars are undesirably big to support or falsify tensorial spatial correlations. On the other hand, the results show that vector correlations are disfavored, in all significance metrics we consider. This improves the overall constraints on the space of vector correlations in PTA Arzoumanian et al. (2021b).

We visualize the mean correlation curves that resulted from the sampling together with the data in figure 3.

Subluminal vector correlations in the SGWB are shaped at large angles by its dominant multipoles Bernardo and Ng (2023)

[TABLE]

an order that follows from the formulas in Section IIIA. Near luminal GW speeds for vector correlations are also more strongly disfavored than tensor ones, as the vector correlations diverge at $\zeta\rightarrow 0^{+}$ Bernardo and Ng (2023, 2022). Figure 3 shows that the extra dipolar power, $C_{1}\neq 0$, in the vector modes does not help out in making their correlations a better fit to the data, or at least not better than its tensor counterpart ($C_{1}=0$ and $C_{2}>C_{3}>\cdots$). Of course, we expect vector GW correlations to be irrelevant, for instance, through the arguments raised in Wu et al. (2023) that there are not enough mechanisms that lead to their excitation. We also add that vector modes are diluted by cosmic expansion. On the other hand, our result disfavors the vector correlations in the SGWB not by argument but by comparing directly with data.

The large uncertainty in the correlations clearly benefit the GW monopole, which so far has the largest signal-to-noise ratio among various models Arzoumanian et al. (2021b). Nonetheless, the tensor correlations should eventually manifest in the data as more pulsars are included Pol et al. (2022). The cross correlations enhance with the square of the number of pulsars in a PTA, or, in other words, more pulsars guarantee smaller uncertainties in the cross correlations.

Upper bounds on the graviton mass can be obtained through the PTA correlations regardless of the significance of tensor modes. Particularly this work puts a bound of $m_{\rm g}\lesssim 10^{-22}$ eV on the graviton mass. This is consistent with an earlier result Wu et al. (2023) restricted within the framework of massive gravity Liang and Trodden (2021). Here, they evaluate the significance of tensor correlations per fixed graviton mass, but find that the values fall short to make any strong conclusions in support of massive gravity correlations. On the other hand, this work is conservative to the physics that endows the graviton a nonzero mass Montani and Moretti (2019); Poddar et al. (2022), or rather follows a more framework–liberated direction in the sense that there is no gravity model in particular that is being considered. We examined GW propagation in this way, considering both tensor and vector correlations in the SGWB, by constraining the GW speed and the overall SGWB amplitude directly through the correlation measurements. It should be mentioned that the dispersion relation that links this work to massive gravity can also arise in various other models particularly in the weak field limit relevant for PTAs. In any case, neither present evidence in support of massive gravity or subluminal GWs in the data, and so are only able to give upper bounds to the graviton mass cast by the PTA GW energy scale.

V Conclusions

We have explored subluminal GW modes as a possible source of spatial correlations in a PTA. This adds to the building literature on their feasibility given the available data: Arzoumanian et al. (2021b) studied luminal scalar, vector, and tensor modes; Chen et al. (2021c) on luminal scalar transverse modes; Bernardo and Ng (2022) on subluminal scalar transverse and scalar longidutinal modes; Wu et al. (2023) on subluminal tensors. This work considered subluminal tensor and, for the first time, to the best of our knowledge, subluminal vector modes that may anchor as SGWB correlations within the scope of the data.

Tensor correlations in general turned out to have low significance, which is understandable given the large uncertainties in the correlations data. Vector correlations appeared strongly disfavored, if not ruled out, by the data. Of course the result with vectors may not be surprising as we expect such modes to be diluted by the cosmic expansion, and that there’s not too much known mechanisms that lead to their excitation Wu et al. (2023). Our result on the other hand disfavors vector GW modes directly from the data without relying on other assumptions.

We emphasize that we use here the actual correlation observations to constrain the properties of the SGWB, rather than through the spectral density. The main advantage of which is that the SGWB produces a very distinct correlations signature that cannot be mistaken for systematics. We see this way also forecasts PTA/SGWB science’s precision era as the cross correlations significantly enhance in the data with more pulsars Pol et al. (2022).

We acknowledge that this analysis, and all that came before, are based on the hypothesis that the correlations in PTA are mainly due to GWs. This is reasonable as GWs are real and thus by their wavelike nature it should not come as a surprise to find them in superposition, as in the SGWB. The challenge comes from setting apart the uncorrelated common process, $A_{\rm CP}$, that is confirmed by all PTAs, from the correlations, $\Gamma_{ab}$, which we associate with GWs, since PTA correlation measurements come in the combination $A^{2}_{\rm CP}\Gamma_{ab}$. We are hopeful that the observation of more pulsars Pol et al. (2022) and the recent IPTA checklist Allen et al. (2023) should be able to settle the reservations surrounding this fundamental assumption.

Subluminal GWs open up various prospects. If the evidence for them improves in future data sets, then fundamental questions about the microscopic and cosmological theory of gravity naturally follow. Alternatively, should subluminal GW propagation be ruled out in PTA, then this may strongly put to rest further questions about the GW speed, since the PTA GW band and energy scale are orders of magnitude away from LIGO/Virgo characteristic scales. Tangential to this, the sources of the SGWB, which exist at very high energies, remain to be given further study once there is enough resolution in the data. On a practical side, incorporating the theoretical uncertainties in the correlations Allen (2023); Allen and Romano (2022); Bernardo and Ng (2022) in the SGWB search pipelines of PTAs remain for future work.

Note Added. Wu et al. (2023) appeared in the arXiv while we were writing ours.

Acknowledgements.

This work was supported in part by the National Science and Technology Council (NSTC) of Taiwan, Republic of China, under Grant No. MOST 111-2112-M-001-065.