Novel tests of gravity using nano-Hertz stochastic gravitational-wave background signals

Gravity theories that modify General Relativity in the slow-motion regime can introduce nonperturbative corrections to the stochastic gravitational-wave background (SGWB) from supermassive black-hole binaries in the nano-Hertz band, while not affecting the quadrupolar nature of the gravitational-wave radiation and remaining perturbative in the highly-relativistic regime, as to satisfy current post-Newtonian (PN) constraints. We present a model-agnostic formalism to map such theories into a modified tilt for the SGWB spectrum, showing that negative PN corrections (in particular -2PN) can alleviate the tension in the recent pulsar-timing-array data if the detected SGWB is interpreted as arising from supermassive binaries. Despite being preliminary, current data have already strong constraining power, for example they set a novel (conservative) upper bound on theories with time-varying Newton's constant (a -4PN correction) at least at the level of Ġ/G ≲ 10^-5 yr^-1 for redshift z=[0.1÷1]. We also show that NANOGrav data are best fitted by a broken power-law interpolating between a dominant -2PN or -3PN modification at low frequency, and the standard general-relativity scaling at high frequency. Nonetheless, a modified gravity explanation should be confronted with binary eccentricity, environmental effects, nonastrophysical origins of the signal, and scrutinized against statistical uncertainties. These novel tests of gravity will soon become more stringent when combining all pulsar-timing-array facilities and when collecting more data.


Introduction
Pulsar timing arrays (PTAs) offer a unique way to probe gravitational-wave (GW) astrophysics at the nano-Hertz (nHz) scale.In 2020, the NANOGrav collaboration first reported evidence in their 12.5 year dataset [1] for a common spectrum of a stochastic nature, which provided the first hint of a stochastic gravitational wave background (SGWB) signature.However, in these data, there was no statistical evidence for Hellings-Down (HD) correlation pattern, necessary to interpret the signal as a quadrupolar GW background.Remarkably, the more recent PTA data released in 2023 by the NANOGrav [2,3], EPTA (in combination with InPTA) [4][5][6], PPTA [7][8][9] and CPTA [10] collaborations, found evidence for a HD angular correlation, typical of an homogeneous spin-2 GW background and consistent with the quadrupolar nature of GWs in General Relativity (GR) [11].
In principle, also modifications to GR (see, e.g., [83][84][85] for some reviews) may play a crucial role and lead to a different prediction than the standard scaling expected for SMBH binaries (including correlation patterns deviating from HD [86][87][88][89]).Indeed, many known theories beyond GR induce new effects at negative post-Newtonian (PN) orders.Examples include dipole radiation in scalar-tensor theories at −1PN, or the effects of extra-dimensions or a time-varying Newton constant, both at −4PN (see e.g.[85,90,91]).Given that BH binaries in the PTA band have extremely large orbital separation during the early-inspiral phase, one could expect negative PN modifications to GR in this regime to play a much more relevant role than in the coalescence phase typically explored by ground-and space-based detectors [92][93][94][95][96].
In this work, we focus on modified GR theories preserving the quadrupolar nature of the gravitational signal, i.e that would be compatible with the HD correlation pattern and with LIGO-Virgo-KAGRA constraints on the GW polarizations [92].For example, widely considered scalar-tensor theories introduce new dissipative degrees of freedom (e.g., scalar fields) modifying the dynamics of the binary (through extra energy fluxes) but without coupling to matter nor to the pulsar signal [85,90].These theories propagate GWs with only the two standard GR polarizations [97,98] but introduce negative PN corrections in the GW signal.We therefore analyze the impact of generic effects at negative PN orders in light of the recent results of the PTA collaborations.We show that the recently detected SGWB allows for novel tests of GR.On the one hand, slow-motion modifications to GR can alleviate the current tension in PTA data and, on the other hand, we can use current data to place stringent upper bounds on putative negative-PN modifications, e.g.theories predicting a time-varying Newton's constant, G(t), in an unconstrained region of their parameter space [99].

SGWB spectrum in modified gravity theories
The Hamiltonian of a binary, corresponding to the centre of mass binding energy reads, at the Newtonian level, where m 1 , m 2 are the masses of the objects, r is the relative orbital separation, v is the relative orbital velocity, and µ is the reduced mass.Here we adopt a theory agnostic approach and consider model-independent modifications to GR.In particular, we consider dissipative corrections to the rate of change of the binding energy, at generic PN order 1 : where m = m 1 + m 2 is the total mass and ĖGR = 32 5 Gµ 2 ω 6 r 4 is the standard GR energy flux.This approach was widely considered in the literature in the contest of the parametrized post-Einsteinian (ppE) framework, wherein the parameter q is theorydependent and accounts for the dominant correction, while the parameter B can be mapped to the coupling constants of a given theory [91,[100][101][102].Nevertheless, in the ppE formalism these corrections are always assumed to be much smaller than the GR term, and hence treated perturbatively.Here we remove this assumption, as we are interested in a seldom explored regime where these terms might dominate over the GR ones at the nHz scale.To compute the impact of this modification in the energy spectrum of GWs, we can simply express the latter as As we focused on purely dissipative corrections to GR, the binding energy is not modified, and one can therefore assume the standard Kepler's law relating the orbital radius to the frequency at the leading order, r = (Gm/ω 2 ) 1/3 .Using the latter, and combining Eqs.(2.1) and (2.2), we can compute the rate of change of the orbital frequency, Finally, the modified GW spectrum reads In the B = 0 limit this equation coincides with the standard GR one, where the only dissipation channel is ordinary GW emission, and leads to the standard scaling , where f = ω/π is the GW frequency.On the other hand, if the modified scaling reads Modifications of the energy flux for negative PN orders (q < 0) lead to a steeper scaling, as expected.In particular, a dominant −2PN correction (q = −2) would give Ω GW ∝ f 2 , compatible with the peak of NANOGrav15 posteriors for the tilt [2], as we shall also discuss later on.Taking q = −2 as an example, from Eq. (2.6) one has (2.8) On the other hand, the requirement that the same correction remains perturbative during the inspiral of a coalescing binary imposes B ≪ (Gmω ISCO ) −2q/3 ∼ 0.03 (for q = −2 and where ω ISCO is the reference frequency of the innermost stable circular orbit), regardless of the binary mass.Thus, there is a wide range in which the coupling B satisfies the inequality (2.8) while introducing a small correction in the coalescence phase.This range becomes even wider for more negative PN corrections.Note that this general formalism includes several classes of gravity theories [90], but also several environmental effects that provide dissipation mechanisms at negative PN orders [103][104][105], leading to a steeper scaling [16,20,106].Examples include accretion and dynamical friction (q = −5.5),stellar scattering (q = −5), and interaction with circumbinary gas (q = −3.5).Finally, large eccentricity (e ≳ 0.6) also provides a steeper scaling [107,108] which can be fitted by multiple power-law phases with q ≈ −7.1 and q ≈ −2.8.

The case of varying G
Before delving into an actual confrontation with PTA data in the next section, let us discuss the strong constraining power of current data through a heuristic argument and considering a specific example.Alternatives to GR can violate the strong equivalence principle and in some case break local invariance [83], while remaining (mostly) quadrupolar in nature.Some of them predict a spacetime variation of the effective Newton constant, for example mediated by a scalar field on cosmological scales [109] or in the presence of energy leakage into small extra dimensions [110,111].Here we apply our formalism to the simplest and most studied case, promoting G to a function of time, G(t), with G 0 being its present value [112,113].Several constraints to the first derivative of the Newton constant, Ġ, have been considered at different scales.At cosmological scales, bounds come from Big Bang nucleosynthesis [114][115][116] and cosmic microwave background [117], which estimated Ġ/G 0 ≲ 10 −12 yr −1 , while in the solar system the most stringent constraint comes from a detailed analysis of Mercury's orbits, yielding Ġ/G 0 ≲ 10 −14 yr −1 [118].Nevertheless, the former constraint assumes a linear scaling of G(t) throughout the entire cosmic history, while the latter is obtained in the solar system at zero redshift.On the other hand, GWs offer a unique probe of local variation Ġ at intermediate epochs [113,119], being thus complementary to the aforementioned constraints.Depending on the specific GW source, one can place constraints at different redshift z.In particular, PTAs are sensible to systems at z ≈ [0.1 ÷ 1] [13], and hence they can provide complementary constraints to the lowredshift ones which were placed using binary pulsars [112,120,121] and LIGO [92].On the other hand, the future space mission LISA will be able to provide bounds at similar redshifts [113,119,122,123].Here we show that using the recent PTA detection of the SGWB it is possible to place competitive bounds already with current data.The GW emission power in these theories depends on the time-varying Newton constant, where M c is the binary chirp mass.Using Eq. (2.1) (with G → G(t)) and the balance law, one can compute an expression for the secular variation of the GW frequency valid at all orders in Ġ [119]: The first term above is the standard GR one, while the second term is the modification due to the variation of G. Using the above relations, one can straightforwardly compute Eq. (2.3).If the first term dominates over the second one in Eq. (3.2), one obtains the standard scaling Ω GW ∝ f 2/3 , while in the opposite regime, i.e. if Ġ is large enough, one obtains the scaling which corresponds to q = −4, i.e. to a −4PN order effect, in agreement with the ppE mapping of this correction in the perturbative regime [90,113].Hence, theories with a time-varying G can be directly mapped into our generic framework.The obtained scaling (3.3) deviates more than 3σ from the recent NANOGrav15 measurement, which can in turn be used to place constraints on Ġ.Assuming for simplicity that all binaries have same redshift and equal masses, m 1 = m 2 = m/2, one finds that the spectrum (3.3) is obtained if which is therefore excluded by PTA measurements.

Confrontation with NANOGrav15 dataset
After the heuristic analysis above, we now focus on the most constraining PTA dataset, which is the one reported by the NANOGrav collaboration [2].Assuming the signal to be a power-law (PL), this corresponds to a fraction energy density in GW today as [124] where H 0 is the current Hubble rate, f yr ≡ (1 yr) −1 ≃ 32 nHz, whereas A and n T are the amplitude and spectral tilt, respectively.In Fig. 1, we show the posterior distribution for the SGWB abundance Ω GW (f ) (gray violins), obtained by the NANOGrav collaboration fitting their data with the inclusion of HD correlations, in the first 14 bins.We perform a maximum likelihood analysis to derive bounds on the tilt of the spectrum as a function of frequency, and compare this with the circular inspiral prediction Ω GW ∝ f 2/3 .The full log-likelihood is computed as  For f > f b we always assume n T = 2/3 as predicted by GW-driven, circular SMBH binaries in GR.Left: The low-frequency tilt n T is left free to vary.As f b increases, the constraint becomes more stringent, and converges towards the single power-law result.The colored vertical lines indicate the 2σ range for n T (f < f b ) while the dashed vertical lines indicate n T = 2/3 and 10/3 corresponding to GR and −4PN corrections, respectively.Right: The low frequency tilt is fixed assuming a given negative qPN correction, n T = 2/3(1 − q) for q = −1, . . ., −6.The posterior distribution for f b indicates the preference for negative PN corrections (i.e.steeper spectra than ∝ f 2/3 ).As a consequence, we always find lower bounds on f b , while for q ≲ −3 upper bounds are also obtained.The light vertical lines indicate the location of the NANOGrav frequency bins.
where θ is the parameter vector characterising the models (e.g., amplitude and tilt for the power-law model), while p k is the likelihood for each frequency bin k.
First, we assume that the SGWB spectrum is described by a power law in a certain number of the lowest frequency bins,hereafter denoted as N bin .We show the resulting 2σ range for the spectral tilt in Table 1 (PL case).
We can also assume the SGWB spectra to be described by a broken power-law (BPL) where the low-frequency tilt (f < f b ) is fixed assuming a given negative PN correction, n T = 2/3(1 − q), while Ω GW ∝ f 2/3 for f > f b .In our model (2.5), this occurs naturally: 1.9  1.We report the range of n T at 2σ obtained fitting the first N bin frequency bins with a power-law (PL), or with a broken power-law (BPL), where we set n T = 2/3 only at f i > f b .The first frequency bin is at f 1 = 1.98 nHz.
for fixed q and B, f b corresponds to saturating Eq. (2.6), namely 2q /(πGm).The corresponding 2-σ range for n T , assuming different values of f b (corresponding to the location of each NANOGrav15 bin) are also shown in Table 1.This second version of the bounds are slightly more stringent, as additional information is brought by considering all NANOGrav15 posteriors, while assuming the high frequency range is described by the GW prediction.In Fig. 2 (left panel), we show the corresponding posterior distribution for n T (f < f b ), assuming the BPL model.
One can also perform a different analysis and focus on a single PN correction (i.e.fixing q) and derive the best fit parameters for the location of the break frequency f b .The posterior distribution for f b shown in the right panel of Fig. 2 has always significant support in the NANOGrav band, thus indicating the preference for negative PN corrections (i.e.steeper spectra than f2/3 ).In all cases, we find lower bounds on f b , while upper bounds are also obtained when q ≲ −3.
Constraints on f b can be directly translated into constraints on the parameter B using Eq.(2.6).We show this in Fig. 3, assuming for simplicity that all binaries have the same total mass.The absence of support for B = 0 in all cases with q < 0 is due to the existing tension in PTA data with the standard scenario Ω GW ∝ f 2/3 .
For q = −4 the posterior on B can be mapped into a bound on Ġ/G 0 , assuming the latter to be approximately constant at the source's redshift z ∼ 0.5 dominating the signal [13], see Fig. 4.Even in the most conservative scenario (m = 10 9 M ⊙ ), the bound on Ġ/G 0 is comparable to the projected bounds that LISA is expected to place using quasi-monochromatic sources [119].Due to the Ġ ∝ m 5/3 dependence, in more optimistic scenarios (m < 10 9 M ⊙ ) the bound is only slightly less stringent than what will be achievable in the LISA era by detecting SMBH coalescences and extreme mass-ratio inspirals [113,122,123] at similar redshift (z ≈ [0.1 ÷ 1]).Overall, the bound on Ġ/G 0 is also several orders of magnitude better than the current bounds with LIGO/Virgo blackhole 2 binaries [90,92], consistently with the fact that a small-velocity nonperturbative modification in the nHz band can become perturbatively small in the relativistic regime.

Conclusions and outlook
We argued that the recent groundbreaking PTA detection of a SGWB provides novel avenues to test gravity in the slow-motion regime.The novel constraints derived here are already competitive with those that the future LISA mission will place using nonrelativistic binaries.More accurate and stringent constraints should be derived by taking into account the mass distribution of SMBH binaries and its uncertainties.Furthermore, we expect that the upcoming joint analysis involving all collaborations within the International PTA framework will soon strengthen the constraints derived in this work.Likewise, the PTA posteriors on the tilt will shrink in time as longer duration datasets will be analyzed.Following the estimates provided in Ref. [126], one may expect the uncertainty on the spectral tilt to scale with the observation time as σ n T ∝ T −0.9(5−n T ) obs (∝ T −0.4 obs ) in the intermediate (weak) signal regime, with additional improvements given by the growing number of observed pulsars.
Depending on the specific modified gravity theory, binary pulsars [85,127] can provide more stringent constraints on negative PN terms than PTAs.However, a direct comparison is challenging, since the coupling B can be source-dependent, for example being zero for pulsars and nonzero for SMBHs, as it happens in several theories for black-hole dipolar emission [91,127].As expected, PTA constraints are more stringent for more negative PN terms.For example, even if a −1PN modification were excluded at each frequency bin, the bounds derived from Eq. (2.6) at the nHz scale for q = −1 would be at most the level of B ≲ 6 × 10 −4 .This is more stringent than those that can be obtained with LIGO at design sensitivity, but much less stringent than the projected bounds on B with LISA sources [91,122,123].
In addition to placing competitive constraints on modified gravity, a more ambitious possibility is to invoke slow-motion, beyond-GR effects to alleviate or solve the tension in the latest PTA data if the detected SGWB is interpreted as arising from SMBH binaries.While we showed that this is in principle possible without violating current constraints, one should be careful with degeneracies with astrophysical effects [13][14][15][16][17][18][19]21], including large eccentricity [107,108], and environmental modifications such as stellar scattering, interaction with circumbinary gas, and dynamical friction [16,20,106].In this context, since our framework can accommodate both beyond-GR and environmental effects, it can be used to disentangle these modifications (which often enter at different negative PN order) through dedicated Bayesian inferences with more constraining future datasets.Finally, future work can be devoted to confront our general model with the hypothesis that the signal is due to, or contaminated by, sources of cosmological origin [14].

Figure 1 .
Figure1.The posterior distribution for Ω GW in each frequency bin f i reported by NANOGrav, including HD correlations (gray violins).We also show the best fit spectra obtained by inserting a specific PN correction in the model, i.e. assuming a BPL spectrum with n T = 2/3(1 − q) at f < f b and n T = 2/3 at f > f b .The black dashed line shows the best fit assuming GR scaling Ω GW ∝ f 2/3 .

Figure 2 .
Figure 2. Fit of the NANOGrav15 dataset, assuming a BPL model with break frequency at f b .For f > f b we always assume n T = 2/3 as predicted by GW-driven, circular SMBH binaries in GR.Left: The low-frequency tilt n T is left free to vary.As f b increases, the constraint becomes more stringent, and converges towards the single power-law result.The colored vertical lines indicate the 2σ range for n T (f < f b ) while the dashed vertical lines indicate n T = 2/3 and 10/3 corresponding to GR and −4PN corrections, respectively.Right: The low frequency tilt is fixed assuming a given negative qPN correction, n T = 2/3(1 − q) for q = −1, . . ., −6.The posterior distribution for f b indicates the preference for negative PN corrections (i.e.steeper spectra than ∝ f 2/3 ).As a consequence, we always find lower bounds on f b , while for q ≲ −3 upper bounds are also obtained.The light vertical lines indicate the location of the NANOGrav frequency bins.

Figure 3 .
Figure 3. Posterior distribution of B derived assuming only one PN correction is active at the time.We left the explicit dependence on the total mass of the binaries, conservatively normalised to the value of m = 10 9 M ⊙ .These constraints are derived from the BPL posteriors on f b in Fig. 2.