Finite-distance effects and scalar suppression in massive gravity

Introduction

Pulsar timing arrays (PTAs) search for gravitational waves (GWs) by monitoring the arrival times of electromagnetic pulses from millisecond pulsars over a period of many years [1, 2]. A passing GW perturbs the observed frequency of the pulse, inducing a redshift. The corresponding timing residuals are obtained from integrating this redshift, and can be correlated between different pulsars to search for a stochastic GW background (SGWB). For an isotropic background, the timing correlations between a pair of pulsars are described by polarization dependent GW energy density spectra and overlap reduction functions (ORFs), which generally depend on frequency, angular separation, and pulsar distances.

In general relativity (GR), the tensor ORF for an SGWB reduces to the Hellings–Downs (HD) correlation [3] in the infinite-distance approximation, making it dependent only on the separation angle between pulsars. Theories beyond GR can modify this correlation pattern in various ways, including its shape and dependence on other parameters, such as frequency or group velocity of the GW, $v_g$. PTA observations have also been used to constrain deviations from GR, including the presence of non-Einsteinian polarization modes [4, 5]. PTA data sets are thus a rich avenue for testing alternative theories of gravity.

One theory of particular interest is massive gravity, in which the graviton has a nonzero mass, $m_g \ne0$. In such a theory, a massive spin-2 field generically carries five polarization degrees of freedom: in addition to the two helicity-$\pm2$ tensor modes of general relativity, it also includes two helicity-$\pm1$ vector modes and one helicity-0 scalar mode. The effects of massive gravity are most naturally relevant on scales comparable to the graviton Compton wavelength. A graviton mass of cosmological order, $m_g \sim H_0$, where $H_0$ is the present-day Hubble parameter, corresponds to a Compton wavelength comparable to the Hubble scale and can therefore affect late-time cosmology [6]. Certain branches of massive gravity admit self-accelerating solutions, potentially providing an alternative to a cosmological constant [7, 8]. However, constructing a consistent theory of massive gravity is highly nontrivial: nonlinear extensions of Fierz–Pauli massive gravity suffer from a pathological sixth degree of freedom, the Boulware–Deser ghost [9]. More recently, the construction of ghost-free nonlinear massive gravity by de Rham, Gabadadze, and Tolley (dRGT) [10, 11], together with its extension to bigravity [12], have established a consistent framework for studying massive spin-2 fields (see Refs. [13, 14, 15, 16] for reviews). This framework has since motivated a broad literature on cosmological backgrounds and self-accelerating solutions in massive gravity and bigravity [6, 7, 8, 17, 18, 19, 20, 21]. Other work has studied perturbative stability, extensions to cosmological settings, and the viability of different branches of solutions [22, 23, 24, 25, 26, 27, 28, 29, 30], while observational analyses have placed constraints on bimetric and massive-gravity cosmologies [31, 32, 33, 34, 35, 36, 37].

Pulsar timing arrays provide a way to test a nonzero graviton mass through its effects on GW propagation and polarization content. Previous PTA studies have explored subluminal GW propagation [38, 39, 40, 41, 42], modifications to the HD correlation [43, 44, 45, 46], and non-Einsteinian polarizations [4, 5, 47, 48, 49, 50]. These effects enter the observed timing correlations through both the pulsar response to each polarization and the relative power carried by the different sectors.

At the level of the pulsar response, an important consideration is whether the pulsar-dependent terms can be neglected. These terms are often averaged over or neglected in ORF calculations in PTA analyses [3, 39, 43, 51, 52, 53]. Although this approximation is generally appropriate for the GR tensor modes, it does not necessarily hold across the propagation regimes and polarization sectors relevant to massive gravity. Recent studies have examined finite-distance effects and the limitations of the infinite-distance approximation [45, 46, 50, 54, 55]. These considerations motivate retaining the full finite-distance response and adopting a normalization that remains consistent when the pulsar terms are included.

The relative sector energy densities are particularly important near the massless limit. As we demonstrate below, the scalar ORF approaches a finite and nonzero correlation as $m_g \to0$, whereas the vector ORF vanishes. A model intended to recover GR must also specify how the scalar energy density behaves in this limit. We introduce a phenomenological scalar energy-density prescription. We establish sufficient conditions for uniform scalar decoupling and recovery of the finite-distance GR tensor correlation. This work develops the theoretical framework and studies its phenomenology; parameter estimation and observational constraints using PTA data are left for future work.

This paper is organized as follows. In § 1.1, we introduce the formalism for the analysis of the SGWB in massive gravity. In Sec. III, we study the finite-distance effect and the massive and massless limits. In § 1.3, we introduce a prescription for suppressing the scalar contribution. In § 1.4, we discuss our findings and the possibility of future exploration. In what follows, we use natural units, setting $c = \hbar= k_B = 1$, and we use the $+---$ metric convention.

Formalism

Dispersion and polarization conventions

When the region of observation is small compared with the background curvature scale, we may express the metric perturbation $h_{\mu\nu}$ as a plane wave [43, 56]

where $\omega= 2\pi f$ is the angular frequency and $\mathbf{k}$ is the spatial wave vector. We adopt the Fourier normalization and integration measure of Ref. [43]. In this convention, the factor $2/|\mathbf{k}|$ is chosen so that after performing the integral, the measure reduces to $df\,d^2\hat{\Omega}$. We have also used the notation $kx \equiv k_\mu x^\mu$, where $x^\mu= (t,\mathbf{x})$ is the spacetime coordinate at which the perturbation is evaluated and $k^\mu= (\omega,\mathbf{k})$ is the wave four-vector of the metric perturbation. The measure $d^4 k$ denotes integration over all components of $k^\mu$, while the delta function restricts the modes to the mass shell via the dispersion relation,

In general massive-gravity cosmologies, the dispersion relation does not necessarily reduce to this form. In a Friedmann–Lemaître–Robertson–Walker background, the equation of motion can contain a time-dependent propagation speed and an effective mass associated with GWs that need not equal the graviton mass parameter $m_g$, such as in extensions of quasi-dilaton theories in massive gravity [57]. As an example, after neglecting terms related to cosmological background evolution, the dispersion relation for the propagating tensor modes takes the form [58]

where $M_{\mathrm{GW}}(t)$ is the effective mass of the tensor modes, $a(t)$ is the scale factor, $K$ is the spatial curvature, and $c_g(t)$ is the tensor sound speed, whose square multiplies the spatial gradient term. The distinction between $M_{\mathrm{GW}}$ and $m_g$ is explicit in quasidilaton and extended quasidilaton massive gravity, where the effective GW mass can differ from $m_g$ [59]. In the present analysis, we assume $a = 1$, $c_g = 1$, a flat background ($K = 0$), and $M_{\mathrm{GW}} = m_g$, leading to the dispersion relation in Eq. (2).

Unlike in GR, the frequency cannot be arbitrarily small: the dispersion relation requires $\omega\geq m_g$, with the minimum value $\omega_{\min} = m_g$ reached when $|\mathbf{k}| = 0$. We define the group velocity of the GW

The physical group velocity is obtained from the full dispersion relation and is not in general equal to $c_g$; the two coincide only when the $c_g^2|\mathbf{k}|^2$ term dominates the massive term.

We express $h_{\mu\nu}$ as

where we have grouped the helicity states into their sectors, $T = \{+2,-2\}$, $V = \{+1,-1\}$, and $S = \{0\}$, following standard polarization decompositions [43, 49, 56]. Polarization tensors $e^{(i)}_{\mu\nu}$ can be constructed using a basis of polarization vectors adapted from $\hat{\Omega}$, the propagation direction of the GW, which in general is expressed as $\hat{\Omega} = (\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta)$.

We use the massive spin-1 polarization vectors [43]

These vectors are transverse to $k_\mu$, satisfying $k_\mu\epsilon^\mu_n=0$. We construct the massive spin-2 polarization tensors from the spin-1 polarization vectors via Clebsch–Gordan coefficients [60, 61],

These tensors are transverse $k^\mu\epsilon^{(i)}_{\mu\nu}=0$ and traceless $\eta^{\mu\nu}\epsilon^{(i)}_{\mu\nu}=0$ for each helicity $i$. Transversality follows from that of the polarization vectors. In our phase convention, $\eta^{\mu\nu}\epsilon^+_\mu\epsilon^-_\nu=\eta^{\mu\nu}\epsilon^0_\mu\epsilon^0_\nu=-1$, so the relative minus sign in the helicity-0 tensor ensures tracelessness. Additionally, they are unit normalized, $\epsilon^{(i)*}_{\mu\nu}\epsilon^{(j)\mu\nu}=\delta_{ij}$.

Pulsar timing array response

We now consider the response of a pulsar signal to a massive GW propagating along the direction $\hat{\Omega}$. The GW perturbs the spacetime along the line of sight of the pulsar signals, originally emitted with a constant frequency $\nu_0$, as they travel from their source pulsars to our telescopes. The wave-vector of the photon pulse is explicitly $\hat{p}^{\mu}_r=(1,-\hat{p}_r)$, where $\hat{p}_r$ points from Earth to the pulsar $r$. In the barycenter reference frame of the solar system, the frequency shift is characterized by the redshift $z(t)\equiv(\nu_0-\nu(t))/\nu_0$, following the PTA convention of Refs. [43, 51]. PTA collaborations directly measure the times of arrival of pulsar signals, which correspond to the integral of $z(t)$ [51, 62], otherwise known as the timing residual. In Fourier space, the redshift for signals emitted by pulsar $r$, located a distance $L_r$ away, can be written as

where the tilde on the redshift indicates that it is a quantity in Fourier space, and the bracketed factor is the difference between the pulsar term and the Earth term, which comes from the retarded waveforms. $F^{(i)}_r(\hat{\Omega})$ is the geometrical projection factor, which in our decomposition is equal to [43]

The second and third terms come from the time component of the metric perturbation in the redshift response. These terms vanish in the synchronous gauge, as seen in Ref. [49], but they are nonzero in the unitary gauge that we adopt.

To detect an SGWB, we look at the cross-correlation of the residuals across all pulsars, which we can obtain via the two-point correlation of the total redshift. We assume that there are two pulsars $a$ and $b$, located at distances $L_a$ and $L_b$ from the solar-system barycenter in directions $\hat{p}_a$ and $\hat{p}_b$. We assume that pulsar $a$ is situated along the $z$-axis and pulsar $b$ is located in the $x$–$z$ plane at an angle $\xi_{ab}$ from pulsar $a$. Then, the pulsar directions are explicitly

Assuming that the SGWB is stationary and isotropic, the corresponding two-point correlation function of the sky-integrated redshift, $\tilde{z}_a(f)\equiv\int d^2\hat{\Omega}\,z_a(f,\hat{\Omega})$, is [43, 51]

Here $\Omega_i$ is the dimensionless energy density per logarithmic frequency for each helicity $i$, $\tilde{\Gamma}^{i}_{ab}$ is the corresponding normalized ORF, and $\beta_i$ is the normalization factor. The raw signal $\langle\tilde{z}^{*}_a(f)\tilde{z}_b(f')\rangle$ itself is independent of normalization convention, since the normalized ORF is being divided by its normalization. The factor $\delta(f-f')$ follows from stationarity, and isotropy removes dependence on the absolute orientation of the pulsar pair via rotational symmetry, leaving the angular dependence only through $\xi_{ab}$.

Within the tensor and vector sectors, we assume equal power in the two helicity modes. This follows from taking the stochastic background to be parity invariant, so that there is no net circular polarization in either the tensor or vector sector. We also assume that there is no cross-correlation between distinct helicities or sectors, motivated by the fact that the tensor, vector, and scalar modes decouple with an appropriate gauge-fixing term, allowing the quadratic kinetic terms to be diagonalized [13]. We can write the dimensionless energy density and ORF for the sectors themselves, with $\Omega_I=\sum_{i\in I}\Omega_i and unnormalized ORF $\Gamma^{I}_{ab} = \sum_{i \in I} \Gamma^{i}_{ab}$. Thus, with the equal power assumption, $\Omega_{+2} = \Omega_{-2} = \Omega_{T}/2$ and $\Omega_{+1} = \Omega_{-1} = \Omega_{V}/2$, while $\Omega_{0} = \Omega_{S}$. The signal is proportional to

leaving the scalar contribution with a relative factor of 2 compared to the tensor and vector contributions.

Following Ref. [43], we define the ORF for each polarization sector $I$ as

We note that the phases have opposite signs because the redshift from pulsar $a$ is complex conjugated in $\langle\tilde{z}^{*}_{a}(f)\tilde{z}_{b}(f')\rangle$, while the redshift from pulsar $b$ is not. If we cannot distinguish among the different polarization modes, we observe the angular correlation as an effective ORF normalized relative to the tensor sector [43], which has the following form

Here, we notice that the ratios $\beta_{T}/\beta_{V}$ and $\beta_{T}/\beta_{S}$ cause the normalization factors for the vector and scalar modes to cancel out in their ORFs. We may factor $\beta_{T}$ out of the expression and equivalently write the effective ORF as

Alternatively, if the combined correlation were interpreted as an effective vector or scalar signal, we could analogously factor out $\beta_{V}$ or $\beta_{S}$, with the relative energy densities expressed using $\Omega_{V}$ or $\Omega_{S}$ in the denominator, respectively.

Finite-distance effect

The finite-distance effect enters through the difference between the metric perturbation evaluated at the pulsar and at the Earth. Following Ref. [43], we define the difference in metric perturbation as

Using $x_{e}=0$, $x_{p}=L_{r}\hat{p}_{r}$, $t_{p}=t_{e}-L_{r}$, and writing the metric perturbation as a plane wave, this difference becomes

The factor in square brackets contains both the pulsar term and the Earth term, which we notate as

The ORFs therefore contain the product

The four terms in Eq. (19) are the pulsar–pulsar, pulsar–Earth, Earth–pulsar, and Earth–Earth contributions, respectively. We refer to the first three collectively as the finite-distance terms because they retain the dependence on the pulsar distances through $\mathcal{P}^{*}_{a}$ and $\mathcal{P}_{b}$. We use the term “finite-distance effect” for the combined modification to the ORF that results from the retention of these three terms. In the usual large-$fL$ limit, their rapidly oscillating phases in the angular integral cause them to average away, but this behavior is not uniform as $v_{g}\to0$ or as the two pulsar phases become coherent.

It is useful to assign these four terms their own ORFs. For each sector $I\in\{T,V,S\}$, we write

where the superscript $a^{*}$ denotes that the term involves pulsar $a$ and is complex conjugated, the superscript $b$ denotes the term involves pulsar $b$, and the superscript $0$ denotes the Earth term. Explicitly, each term is defined as

where $x\in\{0,a\}$ and $y\in\{0,b\}$ identify the Earth or pulsar contribution in the response of each member of the pair. This naturally extends the notation $\Gamma_{0,I}$ of Ref. [43] for the Earth–Earth contribution, denoted here by $\Gamma^{00,I}_{ab}$.

Finite-distance normalization

We first revisit the normalization convention for the tensor mode in GR. The usual normalization is one where $\widetilde{\Gamma}^{\mathrm{HD}}_{ab}$ is normalized to 1 in the case of the autocorrelation where pulsar $a=b$, and is normalized to $1/2$ in the case of the cross correlation, where pulsar $a \ne b$ [45, 63]. The normalized HD correlation is thus written as

The $\delta_{ab}$ term introduces a discontinuity between the coaligned cross-correlation and the autocorrelation. This distinction arises from the pulsar terms, as we will see.

We consider the finite-distance factor in the GR limit, from Eq. (19). For two distinct pulsars, the pulsar phases are rapidly oscillatory in the large-$fL$ limit. The pulsar terms, i.e. those that contain $\mathcal{P}_a^*$ or $\mathcal{P}_b$, average away in the sky average integral, leaving only the Earth–Earth term, 1. This is the contribution that gives the usual HD curve, whose $\xi_{ab} \to0^+$ limit is $1/2$.

For the autocorrelation, $a=b$, the pulsar–pulsar term becomes $\mathcal{P}_a^*\mathcal{P}_a=1$, so the finite-distance factor becomes

The remaining two terms are still rapidly oscillatory and average away, but the non-oscillatory part of the sum is now 2 rather than 1. Thus the autocorrelation receives an Earth–Earth contribution and a pulsar–pulsar contribution that is equal in magnitude. This is where the $\delta_{ab}/2$ term comes from in the normalized HD curve, as emphasized in Ref. [45].

However, this convention for normalizing the finite-distance GR tensor ORF, which we denote $\Gamma^{\mathrm{GR}}_{ab}$ to distinguish it from HD, is not able to properly account for pulsars that are very similar in distance. The Kronecker delta differentiates exactly identical pulsars from all of the distinct pulsar pairs, but the finite-distance response itself is not discontinuous in this way. The coherence of the pulsar–pulsar term is controlled by the relative phase between the two pulsar terms. Even for two distinct pulsars, this phase can be small if their distances are sufficiently similar. In that case, the pulsar–pulsar term does not necessarily behave like a rapidly oscillating factor, and the ORF may behave more like the autocorrelation.

This can be seen in Figure 1 for the tensor mode in the GR limit, where for small $f\lvert L_a-L_b\rvert$, the pulsar phases are coherent and the response approaches 1, as in an autocorrelation. For large $f\lvert L_a-L_b\rvert$, the pulsar-dependent pieces decohere and the response approaches $1/2$. For distinct pulsars with nearly equal distances, partial coherence allows the coaligned response to take intermediate values between $1/2$ and 1 instead of being fixed at $1/2$ whenever $a \ne b$. To retain this continuous dependence, we normalize $\Gamma^{\mathrm{GR}}_{ab}(\xi_{ab},f)$ by twice the tensor Earth–Earth response evaluated at $\xi_{ab}=0$, $2\Gamma^{00,\mathrm{GR}}$, which smoothly interpolates between the two limiting behaviors, as shown in Figure 1. Here, we note that the Earth–Earth contribution of the ORF no longer has the subscript $ab$ because we have evaluated the ORF at $\xi_{ab}=0$ and the pulsar distances no longer appear in the expression.

Vector (PDF)

With the tensor polarization convention adopted here, the Earth–Earth response in GR satisfies $\Gamma^{00,\mathrm{GR}}=2\pi/3$, and therefore

which is also a standard normalization used in PTA ORF calculations [4, 49, 51]. Although the interchange between this normalization and the $\delta_{ab}$-based normalization is justified when pulsar terms are ignored, it can potentially yield ORFs that disagree in their overall amplitude when pulsar terms are kept. We find that keeping the constant $4\pi/3$ factor as the normalization explicitly preserves the continuous dependence on pulsar separation and distances when pulsar terms are kept.

This is relevant for real PTA observations. For example, pulsars J2017+0603 and J2317+1439 are both located at a distance of 1.57 kpc, with uncertainties $\pm0.31$ kpc and $\pm0.29$ kpc, respectively. Similarly, pulsars J2124−3358 and J1744−1134 are located $0.413\pm0.055$ kpc and $0.4141\pm0.0093$ kpc away. For this pair, the finite-distance normalization gives $\Gamma^{\mathrm{GR}}_{ab}(\xi_{ab}=0)/(2\Gamma^{00,\mathrm{GR}})\simeq0.82$ rather than $1/2$.

Each pulsar pair fixes the distances $L_a$ and $L_b$ and defines its own ORF as a function of angular separation. The correlation entering the analysis is evaluated at the actual observed separation $\xi_{ab}$, while the coaligned limit is obtained by varying the separation with both distances held fixed. This limit is relevant to the normalization of the entire ORF, regardless of the pair’s observed separation. Thus, if we want to compare the ORF values for the different pulsar pairs, especially for a set that contains pairs with similar distances, a shared normalization that puts them on equal footing is necessary.

Our procedure for normalizing the ORF in GR elucidates an analogous normalization procedure for the effective ORF in massive gravity. As shown in Eq. (15), the effective ORF carries the normalization coefficient $\beta_T$. We therefore define $\beta_T$ by using the massive-gravity tensor Earth–Earth contribution,

This quantity can be analytically computed:

Its value in the massless limit is $2\pi/3$, recovering the standard HD normalization $\beta_T=3/(4\pi)$. In the opposite limit, as $m_g/\omega\to1$ and $v_g \to0$, $\Gamma^{00,T}$ approaches the value $4\pi/15$, yielding $\beta_T=15/(8\pi)$. Thus, the tensor Earth–Earth normalization remains finite across the full range of $v_g$.

We adopt the tensor Earth–Earth normalization to compare finite-distance responses across pulsar pairs and propagation regimes. The behavior of these responses depends on both $v_g$ and $\xi_{ab}$. In massive gravity, the pulsar phase contains the factor $1+v_g\hat{\Omega}\cdot\hat{p}_a$. As $v_g$ decreases, its angular oscillations become weaker, and the pulsar-dependent terms no longer necessarily average away over the sky. We examine the massive limit and, following the autocorrelation-first treatment of Ref. [45], the massless autocorrelation limit obtained by first setting $b=a$, $\xi_{aa}=0$, and $L_b=L_a$ at finite $fL_a$, and then taking $v_g\to1$ at fixed $f$ and $L_a$.

ORF expressions

To simplify our expressions, we define the angular projection factors for pulsar $a$ and $b$

Substituting the polarization tensors into the general receiving functions in Eq. (9) gives the following expressions

The intermediate contractions are provided in Appendix A. Factoring out the common denominators $(1+v_g\mu_a)$ and $(1+v_g\mu_b)$ and grouping them with the finite-distance factors of the ORF, since they naturally appear in the phases of the pulsar terms, gives us the angular kernels $\mathcal{N}_I$ for each polarization sector:

The sector ORFs therefore take the common form

for $I\in\{T,V,S\}$. Note that because the receiving function sums are all real, the imaginary component of the integrand comes solely from the pulsar terms. We see that the tensor kernel is independent of $v_g$ because the helicity-$\pm2$ polarization tensors are purely spatial and transverse. Nevertheless, the tensor ORF has $v_g$ dependence through the denominators and the pulsar phases. The vector kernel contains a factor of $1-v_g^2$, that converges to $0$ when $v_g\to1$, causing the vector ORF to vanish in the GR limit. Like the tensor ORF, the scalar ORF remains finite in the GR limit. These ORFs are not analytically integrable in their full general forms, but in certain limits, they permit analytical forms.

Massive limit

We first analyze the massive limit, in which $m_g/\omega\to1$, or equivalently $v_g\to0$ when keeping $\omega$ fixed. In this limit, the angular dependence of the phases of the pulsar terms disappear because the oscillatory part of the pulsar term, $1+v_g\mu_a$ and $1+v_g\mu_b$, both approach 1. The pulsar terms thus approach the following forms,

Taking the real part of the finite-distance product in the ORF integrals, we get

Here we used the identities $e^{\pm ix}-1=\pm2ie^{\pm ix/2}\sin(x/2)$ and $\operatorname{Re}(e^{ix})=\cos x$. $\mathcal{R}_{ab}^0$ is periodic under $fL_a\to fL_a+n$ and $fL_b\to fL_b+m$, for integers $n$ and m. Because the pulsar-distance uncertainties relevant to PTAs generally span many cycles at nanohertz frequencies, the phases $2\pi fL_a and $2\pi fL_b$ are effectively unconstrained modulo $2\pi$. As a result, the factor can vary between pulsar pairs and may take either sign. Its maximum value of 4 occurs when both $fL_a$ and $fL_b$ are half-integers, and its minimum value of $-1/2$ occurs when the fractional parts of $fL_a$ and $fL_b$ are $1/6$ and $5/6$ in either order. Due to the lack of angular dependence of the function, it factors out of the angular integral and as a result, the angular dependence of the ORFs is determined entirely by $\mathcal{N}_T$, $\mathcal{N}_V$, and $\mathcal{N}_S$.

The angular integrals of these kernels give the sector ORFs in the massive limit:

where $P_2(\cos\xi_{ab}) = (3\cos^2\xi_{ab}-1)/2$ is the Legendre polynomial of order 2. All three sectors share this quadrupolar dependence, with identical tensor and vector responses and a scalar response half their magnitude, as illustrated in Figure 2.

We form the effective ORF by combining the three modes, setting the energy density ratios $\Omega_V/\Omega_T = 1$ and $\Omega_S/\Omega_T = 1/2$ and applying the normalization of Eq. (25). We find the massive limit to be

In this limit, the pulsar terms no longer average away, and $R^{0}_{ab}$ modifies the normalized effective ORF, giving a possible enhancement, suppression, or a change in sign, depending on the pulsar phases. If the pulsar terms are instead neglected before computing the angular integral, this factor is replaced by 1 and the result becomes $5P_2(\cos\xi_{ab})/4$.

For the parameters in Figure 3, $R^{0}_{ab}=-0.3523$, so the pulsar terms reverse the signs of all three Earth–Earth responses and reduce their magnitudes to approximately 35%. We can see that even for large pulsar distances,

increasing $m_g$ can prevent the pulsar-dependent terms from averaging away.

Massless limit

We now examine the ORFs in the massless limit $m_g \to0$, or equivalently $v_g \to1$ keeping $\omega$ fixed. We first consider the ORFs obtained by dropping the pulsar terms to illustrate the necessity of the pulsar terms in this limit. In the massless limit, the tensor Earth–Earth contribution to the ORF remains finite and gives the usual HD correlation

The autocorrelation of the massless tensor ORF with the pulsar terms is found to be

We see that the autocorrelation value of $4\pi/3$ is only modified by small correction terms whose leading order scales as $(fL)^{-2}$ for large $fL$. If the tensor modes are the only modes considered in an analysis, then ignoring pulsar terms provides a very good approximation for cross-correlations between distant, sufficiently separated pulsars. For autocorrelations, however, the pulsar–pulsar contribution must still be retained, producing the familiar factor of two in the large-$fL$ limit.

To examine the Earth–Earth vector response near the massless limit, we factor out $1-v_g^2$, obtaining

For fixed $0<\xi_{ab}\leq\pi$, this rescaled limit is finite, and the Earth–Earth vector ORF therefore vanishes as $v_g\to1$. The additional limit $\xi_{ab}\to0$ produces a logarithmic divergence in the rescaled expression. To examine the coaligned configuration, we first set $\xi_{ab}=0$ at finite $v_g$, obtaining

This expression also vanishes when $v_g\to1$. The Cauchy–Schwarz bound $|\Gamma_{ab}^{00,V}|\leq\Gamma_{aa}^{00,V}$ extends this conclusion to all angular separations. Retaining the pulsar terms also makes the rescaled response regular at fixed finite pulsar distance, since

The finite-distance pulsar phase therefore regulates the apparent singularity. For the vector ORF, the autocorrelation is

Following Ref. [45], we set $b=a$ in the full response before taking $v_g\to1$ at fixed $f$ and finite $L_a$, obtaining

where $\gamma_E$ is the Euler–Mascheroni constant and Ci is the cosine integral function, defined by

$\operatorname{Ci}(x)$ tends to 0 as $x\to\infty$. The rescaled finite-distance autocorrelation therefore has a finite massless limit, and $\Gamma_{aa}^{V}=O(1-v_g^2)$. Its large-$fL_a$ behavior is obtained only after taking this limit, yielding the logarithmic dependence on $fL_a$. Although both the Earth–Earth term and

the full expression with pulsar terms vanish in the massless limit due to the factor of $1-v_g^2$, retaining the pulsar terms changes the rate at which the vector mode decouples. This can be observed in Figure 4, where we see that the full vector autocorrelation with pulsar terms and the Earth–Earth term of the vector autocorrelation begin to noticeably differ when $1-v_g<10^{-5}$ for the chosen parameters. Nevertheless, the vector contribution vanishes in the massless limit because of the factor $1-v_g^2$, even though the finite-distance effects can modify its behavior relative to the Earth–Earth approximation.

Vector (PDF)

The same analysis can be done for the scalar sector. In the massless limit, the scalar kernel takes the form

This numerator cancels the apparent singularities in the denominator

When we compute the scalar ORF while ignoring pulsar terms, we obtain a rather simple expression,

In general, the helicity-0 polarization contains both breathing and longitudinal contributions. Although individual longitudinal polarization components diverge as $m_g\to0$, the inverse powers of $m_g$ cancel in the complete receiving function, leaving a finite breathing contribution and a longitudinal correction that is suppressed by $1-v_g^2$, as shown in Appendix A. With pulsar terms retained, the latter vanishes uniformly over propagation directions as $v_g\to1$. The surviving response therefore yields the finite breathing-mode ORF [49, 54, 64].

Although the scalar ORF can be analytically integrated when pulsar terms are included in the $v_g=1$ limit (see Eq. (25) of Ref. [64]), the resulting expression is rather difficult to parse, so we shall analyze the autocorrelation. Setting pulsar $b=a$, the scalar ORF autocorrelation is found to be

which is identical to the tensor autocorrelation in Eq. (36) divided by a factor of three. It is clear from this expression that, like for the tensor mode, the pulsar terms also do not appreciably affect the scalar ORF in the massless limit, with a leading correction proportional to $(fL)^{-2}$ for large $fL$.

Figure 5 illustrates the massless limit ORF for each polarization sector. The tensor mode becomes the HD correlation, with small corrections due to the pulsar terms, the vector mode exhibits a pronounced peak near zero angular separation after the vanishing factor is removed, and the scalar mode approaches the breathing-mode correlation—a sum of the monopole and dipole correlations.

Vector (PDF)

Phenomenological scalar suppression

An important aspect that has not been adequately explored in earlier studies of ORFs in massive gravity is the treatment of the scalar mode as $m_g\to0$. In this limit, the dispersion relation reduces to $\omega=|\mathbf{k}|$, and the effective ORF should recover the ORF in GR. However, in linearized massive gravity, the helicity-0 mode remains coupled to matter sources and can remain non-vanishing. In nonlinear completions of massive gravity, the helicity-0 contribution can be suppressed by the Vainshtein mechanism [11, 14, 65, 66], for example. Taking Vainshtein screening as a qualitative motivation, we introduce a phenomenological prescription for suppressing the scalar contribution to the effective ORF.

We parametrize the scalar-to-tensor energy-density ratio as

where $n_S$ is the scalar suppression index that controls how quickly the scalar contribution decouples in the massless limit, and $A_S$ is the scalar amplitude parameter. We will show that $n_S>0$ is sufficient to make the scalar contribution vanish as $m_g/\omega\to0$.

Convergence

We now prove that the effective ORF, with the scalar suppression, converges uniformly in $\xi_{ab}$ to the GR tensor ORF when $m_g\to0$. We start with the scalar sector, whose contribution to the effective ORF is $2\beta_TA_S(m_g/\omega)^{n_S}\Gamma^S_{ab}$. The ORF is the real part of an inner product of the two pulsar response functions, as shown by Eq. (13). The scalar ORF therefore satisfies the Cauchy-Bunyakovsky-Schwarz (CBS) inequality,

for every $\xi_{ab}\in[0,\pi]$. The autocorrelations are nonnegative, since they are sky integrals of the squared moduli of the receiving function and finite-distance terms. At fixed $f$ and finite $L_r$, Eq. (46), together with the property $\sin x\leq x$ for $x>0$, gives

The scalar autocorrelation is therefore bounded near the massless limit and is independent of $\xi_{ab}$. Since $\beta_T\to3/(4\pi)$, for fixed finite $A_S\geq0$ and $n_S>0$ we obtain

for every $\xi_{ab}\in[0,\pi]$. Since the CBS bound is independent of $\xi_{ab}$, and the scalar autocorrelations remain bounded as $m_g\to0$, the scalar contribution converges uniformly to zero for all $\xi_{ab}\in[0,\pi]$.

More generally, the bound implies that the suppression does not need to be in the form of Eq. (47). We see that any nonnegative scalar energy density ratio satisfying

is able to suppress the scalar contribution uniformly in $\xi_{ab}$.

The same argument applies to the vector sector. At fixed $f$ and finite $L_a$ and $L_b$, the right-hand side of Eq. (41) is finite, since it is not dependent on $v_g$. The factor $1-v_g^2$ therefore gives $\lim_{m_g\to0}\Gamma^V_{rr}=0$ for $r\in\{a,b\}$. The CBS inequality gives

for every $\xi_{ab}\in[0,\pi]$. Since the right-hand side of this inequality is independent of $\xi_{ab}$ and tends to zero when $m_g\to0$, the vector ORF itself converges to zero in the massless limit. Additionally, since $\beta_T\to3/(4\pi)$, any bounded $\Omega_V/\Omega_T$ gives

For the tensor sector, we must show that the general tensor ORF uniformly converges to the GR tensor ORF as $m_g\to0$, or equivalently, that the difference $\Gamma^T_{ab}-\Gamma^{\mathrm{GR}}_{ab}$ tends to zero for all $\xi_{ab}$. We see that Eq. (30) can be written using the pulsar factors defined in Eq. (18):

We recall that the tensor kernel $\mathcal{N}_T$ is independent of $v_g$. To examine the remaining dependence on $v_g$, we expand each exponential in its power series:

This expression remains finite at $1+v_g\mu_r=0$. At fixed $f$ and finite $L_r$, the series and its first two derivatives with respect to $v_g$ converge uniformly, since $\mu_r$ is the dot product of two unit vectors, so $|\mu_r|\leq1$. We then Taylor expand the tensor ORF about $v_g=1$:

Since $\Gamma^T_{ab}(1)=\Gamma^{\mathrm{GR}}_{ab}$, subtracting the GR ORF and using the Taylor expansion of $v_g$ with respect to $m_g/\omega$,

gives us

The derivative of $\Gamma^T_{ab}$ evaluated at $v_g=1$ is bounded independently of $\xi_{ab}$, and the same is true for the coefficient bounding the $O((m_g/\omega)^4)$ terms. We see that for sufficiently small $m_g/\omega$,

where $C$ is some finite constant. For the normalized tensor ORF, we write

The CBS inequality also bounds the GR response by its finite autocorrelations from Eq. (36),

Since $\beta_T\to\beta_{\rm GR}=3/(4\pi)$, both terms vanish as $m_g\to0$. The bound is independent of $\xi_{ab}$, establishing uniform convergence of the normalized tensor ORF to the GR ORF.

Revisiting the massless limit

We now compare the rate at which the scalar contribution vanishes with the remaining tensor and vector departures from the GR tensor correlation. When the scalar mode is strictly unsuppressed, with the choice of $n_S=0$ and nonzero $A_S$, the surviving scalar contribution adds a monopole-dipole angular dependence to the effective ORF. This can be seen in Figure 6, where the unsuppressed effective ORF at $v_g=0.999$, represented by the dotted blue curve in the rightmost panel, lies above HD. The corresponding residual in the lower panel approximately follows the scalar angular dependence $3+\cos\xi_{ab}$. For $n_S>0$, the scalar contribution vanishes as $m_g/\omega\to0$. Already at $v_g=0.999$, the choices $n_S=1,2$ give effective ORFs close to HD, as shown by the dotted red and solid purple curves.

Vector (PDF)

We wish to quantify to what extent the scalar contribution is responsible for the remaining deviation from the finite-distance GR tensor correlation. Although the tensor mode approaches $\Gamma^{\rm GR}_{ab}$ as $m_g\to0$, at nonzero mass its response differs through the denominators $1+v_g\mu_a$ and $1+v_g\mu_b$ and the pulsar phases. Similarly, the vector mode can contribute at nonzero mass despite vanishing at the massless limit. We therefore quantify departures from the expected GR signal for the three sectors:

The total deviation can be expressed as $\Delta\widetilde{\Gamma}^{\rm eff}_{ab}=\widetilde{\Gamma}^{\rm eff}_{ab}-\widetilde{\Gamma}^{\rm GR}_{ab}=\Delta\widetilde{\Gamma}^{T}_{ab}+\Delta\widetilde{\Gamma}^{V}_{ab}+\Delta\widetilde{\Gamma}^{S}_{ab}$. Then, integrating these quantities’ squares over $\xi_{ab}$, we obtain the root mean squares (RMS)

This quantity is plotted in Figure 7. For $n_S=0,1$, the scalar contribution dominates remaining deviation from GR for sufficiently small $m_g/\omega$. The asymptotic slopes in Figure 7 follow from the convergence results established earlier.

Vector (PDF)

For fixed $f$ and pulsar distances, the unnormalized tensor deviation scales as $\mathcal{O}((m_g/\omega)^2)$. The normalized tensor deviation, on the other hand, also contains the change in $\beta_T$. To find the further scaling of the deviation, we expand $\beta_T=[2\Gamma^{00}_{T}-1]^{-1}$

Combining this with our results from § 1.3 A, the normalized tensor deviation scales in the massless regime as

Since the logarithmic term, proportional to $(m_g/\omega)^2\ln(\omega/m_g)$, converges more slowly than the term proportional to $(m_g/\omega)^2$ as $m_g/\omega\to0$, it dominates the asymptotic behavior of the tensor deviation and is therefore the leading-order contribution. The vector response scales as $\mathcal{O}((m_g/\omega)^2)$ due to the factor $1-v_g^2=(m_g/\omega)^2$, while the remaining factor approaches a finite function of $\xi_{ab}$ in the massless limit. The scalar ORF approaches a finite limit, leaving its suppression factor to determine the leading scaling. Taking the angular RMS near $m_g\to0$ therefore gives the following scalings:

For $0<n_S<2$, the scalar contribution eventually dominates the remaining deviation. At exactly $n_S=2$, the scalar and vector contributions scale quadratically, while the tensor deviation contains an additional logarithmic factor. For $n_S>2$ the scalar contribution decays with an additional positive power of $m_g/\omega$, as seen in the rightmost panel of Figure 7, where the scalar curve has a greater slope in $m_g/\omega$ and therefore is subdominant to the rest of the effective ORF’s deviation from GR for low $m_g/\omega$.

Conclusions and Discussion

In this paper, we analyzed the effective ORF for an SGWB in massive gravity, retaining finite-distance pulsar terms and suppressing the scalar sector. We showed that the pulsar terms remain important in the massive limit and produce pair-dependent amplitudes. We also introduced a normalization based on twice the tensor Earth–Earth response. In the massless limit, the vector ORF vanishes as $(m_g/\omega)^2$, whereas the scalar ORF approaches a finite, generally non-zero correlation, preventing recovery of the GR tensor correlation unless the scalar contribution is suppressed. We therefore introduced phenomenological suppression $\Omega_{\mathrm{S}}/\Omega_{\mathrm{T}} = A_{\mathrm{S}}(m_g/\omega)^{n_{\mathrm{S}}}$, with scalar amplitude $A_{\mathrm{S}} \ge0$ and suppression index $n_{\mathrm{S}} > 0$, and showed that it sufficiently decouples the scalar contribution to the effective ORF in order to recover the GR tensor ORF as $m_g \to0$. Together, these results provide a consistent treatment of the effective ORF for testing massive gravity with PTA observations near the massless limit.

This work may inspire several directions of research. For example, our phenomenological suppression formalism can be extended beyond PTAs to SGWB searches with ground- or space-based interferometers, where the long-wavelength regime dominates. ORFs for alternative polarizations have already been studied for ground-based detectors [67, 68], and so a consistent treatment of the scalar mode in this context may be applied. Although ground-based searches have not yet detected an SGWB [69], their upper limits already constrain non-tensorial backgrounds [70]. Extending our framework could allow

us to find constraints using existing data and forecasts for future ground-based observations and for the Laser Interferometer Space Antenna [71].

A natural next step is to directly derive the suppression from the Vainshtein mechanism itself [65]. In ghost-free massive gravity, the Vainshtein mechanism is associated with nonlinear interactions of the scalar mode [11, 16, 66]. A theoretical prediction for the scalar energy density would require an explicit calculation that takes into account a specific model exhibiting the Vainshtein mechanism, and ties it to a specific source of GWs, such as a population of SMBHBs or a cosmological source [58, 72, 73, 74].

An assumption we make is that $m_g$ used in Eq. (2) is the mass scale that controls other manifestations of massive gravity, including static Yukawa forces and cosmological evolution. This is not necessarily the case in every nonlinear completion of massive gravity. For example, in extended quasidilaton massive gravity, the effective mass of GWs can differ from the graviton-mass parameter responsible for cosmic expansion [58, 59]. In minimal quasidilaton massive gravity, massive tensor modes can instead have a high-momentum propagation speed equal to the speed of light [75]. A more complete treatment might include deriving the radiative dispersion relation within a nonlinear model.

We have also assumed the isotropy of the SGWB, which allows us to decompose Eq. (11) as above. This assumption is not guaranteed: anisotropies may arise from discrete SMBHB sources, the large-scale distribution of matter, or features of a cosmological source [76, 77, 78, 79, 80]. Anisotropic backgrounds are naturally described in harmonic space, where the multipole content of the SGWB sky distribution enters the timing-residual correlations through generalized, multipole-dependent ORFs that depend on the individual pulsar sky locations, not only on their angular separation [81, 82, 84, 85] [83]. Fisher-matrix methods provide another way to work with real pulsar distributions [86]. More recently, frequency-resolved optimal statistics have been developed to improve the characterization of SGWB anisotropies using PTA observations [87].

Acknowledgments

We thank Marcus Bosca, Neil Cornish, Murman Gurgenidze, Arthur Kosowsky, and Qiuyue Liang for useful discussions related to the paper. C.C. and T.K. acknowledge support from the NASA Astrophysics Theory Program (ATP) Award 80NSSC22K0825 and the National Science Foundation (NSF) Astronomy and Astrophysics Research Grants (AAG) Award AST2408411. T.K. and B.P. acknowledge partial support from the NSF AAG Award AST2307698. T.K. acknowledges partial support from the Shota Rustaveli National Science Foundation of Georgia, FR24-901.

Data Availability

The source code used to reproduce the figures and numerical results in this work is publicly available on Zenodo [88].

Receiving functions

We derive the receiving functions entering the ORFs. We use $\mu_r = \hat{\boldsymbol{\Omega}} \cdot\hat{\mathbf{p}}_r$ throughout. We remind ourselves that the receiving function is

For the polarization vectors given in the main text, we express them as

The transverse vectors satisfy following identities:

Since $\epsilon_\mu^-=(\epsilon_\mu^+)^*$, the sum in Eq. (69) is twice the real part of its $+$ term.

Tensor sector

For $\epsilon_{\mu\nu}^{(\pm2)}=\epsilon_\mu^\pm\epsilon_\nu^\pm$, the time components vanish. Eq. (67) therefore reduces to

The numerator of the helicity sum is

where we have used the identity $2\operatorname{Re}(z^2)=4[\operatorname{Re}(z)]^2-2|z|^2$. Using Eq. (69), together with $\hat{p}_a\cdot\hat{p}_b=\cos\xi_{ab}$, gives

Substituting these into the tensor helicity sum yields

Vector sector

For vector polarization tensors, $\epsilon_{\mu\nu}^{(\pm1)}=(\epsilon_\mu^\pm\epsilon_\nu^0+\epsilon_\mu^0\epsilon_\nu^\pm)/\sqrt{2}$, the contractions are

Substituting into Eq. (67), we obtain

We see that the terms in the numerator combine to form $-(1-v_g^2)\mu_r$, cancelling the inverse power of $\sqrt{1-v_g^2}$. Multiplying the vector receiving functions and summing over the two helicities gives

where the second equality follows from Eq. (69).

Scalar sector

For $\epsilon_{\mu\nu}^{(0)}=(\epsilon_\mu^+\epsilon_\nu^-+\epsilon_\mu^-\epsilon_\nu^+-2\epsilon_\mu^0\epsilon_\nu^0)/\sqrt{6}$, all three terms of Eq. (67) contribute. The contractions are

After substituting into Eq. (67), we get

Combining the last three terms with a common denominator of $\sqrt{6}(1-v_g^2)(1+v_g\mu_r)$ gives us

The numerator of the second term turns out to be a perfect square:

One factor of $1-v_g^2$ therefore cancels, yielding

The scalar receiving function is real, so the receiving function product is simply a square

References