DETECTING MASSIVE GRAVITONS USING PULSAR TIMING ARRAYS

Although a complete quantum version of gravitation has not yet been achieved (Smolin 2003; Kiefer 2006), in the weak field, linearized limit quantization of gravitation can be carried out. One can therefore ask phenomenologically whether or not the graviton in a gravity theory is massive (Rubakov & Tinyakov 2008; Goldhaber & Nieto 2010). Up to now, the most successful theory of gravitation is general relativity (Will 1993, 2006; Stairs 2003; Damour 2009; Kramer & Wex 2009). The quantization of its weak field limit (see Gupta 1952 and references therein) shows that its gravitational interaction is mediated by spin-2 massless bosons, the so-called gravitons. Recently, attention has turned to this question due to both observational and theoretical advances.

On the theoretical side, there is the issue of the vDVZ discontinuity (Iwasaki 1970; van Dam & Veltman 1970; Zakharov 1970), the prediction that for a massive graviton, light deflection is greater than in general relativity by the factor 4/3, no matter how small the mass is. If this were true, despite its paradoxical nature, then classical tests of starlight deflection would rule out massive gravitons. But recently, resolutions have been proposed to the vDVZ paradox, so that massive gravitons cannot be ruled out (Vainshtein 1972; Visser 1998; Deffayet et al. 2002; Damour et al. 2003; Finn & Sutton 2002).

On the theoretical side also, advances in higher dimensional theories have led to graviton mass-like features (e.g., DGP model by Dvali et al. 2000); meanwhile, because of possible issues with standard gravitation from the solar system scale (Anderson et al. 1998, 2008) to the cosmological scale (Sanders & McGaugh 2002), interest has been growing in alternative gravity theories, some of which predict massive gravitons (Dvali et al. 2000; Alves et al. 2010; Eckhardt et al. 2010). There is, therefore, considerable motivation for experimental or observational tests of whether the graviton can be massive.

The literature contains many estimates for an upper limit on the graviton mass, estimates that differ by many orders of magnitude. An early limit of 8 × 10⁴ eV (Hare 1973) was based on considerations of a graviton decay to two photons. At about the same time, Goldhaber & Nieto (1974) used the assumption that clusters of galaxies are bound by more-or-less standard gravity to infer an upper limit of 2 × 10⁻²⁹h₀ eV, where h₀ is the Hubble constant in units of 100 km s⁻¹ Mpc⁻¹. By using the effect of graviton mass on the generation of gravitational waves (GWs) by a binary and the rate of binary inspiral inferred from the timing of binary pulsars, Finn & Sutton (2002) inferred an upper limit of 7.6 × 10⁻²⁰ eV. Choudhury et al. (2004) considered the effect of a graviton mass on the power spectrum of weak lensing; with assumptions about dark energy and other parameters, they estimated a upper limit of 7 × 10⁻³² eV for the graviton mass. Reviews about these techniques can be found, e.g., in Goldhaber & Nieto (2010), Particle Data Group (2008), and Will (2006).

These upper limits are all of value, since they are based on different assumptions about the phenomenological effects of graviton mass. Most of the estimates are based on the effects of graviton mass on the static field of a source, typically assuming a Yukawa potential, though this may not be the general case for a particular theory of gravity, such as fifth force theories or MOND-type theories (Will 1998; Maggiore 2008). In view of the lack of a theory of the graviton, it is important to have upper limits based on different phenomenological implications of graviton mass.

The mass limit of Finn & Sutton (2002) is based on the effect of graviton mass on the generation of GWs, not on their propagation, but the dispersion relation for propagation is also an important independent approach to a mass limit, as has been recently suggested by a number of groups (Will 1998; Larson & Hiscock 2000; Cutler et al. 2003; Stavridis & Will 2009). Questions about this method are timely since the detection of GWs is expected in the near future, thanks to the progress with present ground-based laser interferometers, possible future space-based interferometers (Hough & Rowan 2000; Hough et al. 2005), and pulsar timing array projects (Sallmen et al. 1993; Stappers et al. 2006; Manchester 2006; Hobbs et al. 2009b).

The pulsar timing array is a unique technique to detect nano-Hertz GWs by timing millisecond pulsars, which are very stable celestial clocks. It turns out that a stochastic GW background leaves an angular-dependent correlation in pulsar timing residuals for widely spaced pulsars (Hellings & Downs 1983; Lee et al. 2008). That is, the correlation C(θ) between timing residual of pulsar pairs is a function of angular separation θ between the pulsars. One can analyze the timing residual and test such a correlation between pulsar timing residuals to detect GWs (Jenet et al. 2005). We find in this paper that if the graviton mass is not zero, the form of C(θ) is very different from that given by general relativity. Thus, by measuring this graviton mass-dependent correlation function, we can also detect the massive graviton.

The outline of this paper is as follows. The mass of the graviton is related to the dispersion of GWs in Section 2. The pulsar timing responses to a plane GW and to a stochastic GW background in the case of a massive graviton are calculated in Section 3. The massive graviton induces effects on the shape of the pulsar timing correlation function, which is derived in Section 4, while the detectability of a massive GW background is studied in Section 5. The algorithm to detect a massive graviton using a pulsar timing array and the sensitivity of that algorithm are examined in Section 6. We discuss several related issues and conclude in Section 7.

We incorporate the massive graviton into the linearized weak field theory of general relativity (Gupta 1952; Arnowitt & Deser 1959; Weinberg 1972). For linearized GWs, specifying the graviton mass is equivalent to specifying the GW dispersion relation that follows from the special relativistic relationship:

$$\begin{equation} E^2=p^2 c^2+m^2 c^4, \end{equation} \tag{ 1 }$$

where c is the light velocity, E is energy of the particle, and p and m are the particle's momentum and rest mass, respectively. One can derive the corresponding dispersion relation from Equation (1) by replacing the momentum by p = ℏk_g and the energy by E = ℏω_g, where ℏ is the reduced Planck constant with k_g and ω_g, respectively, the GW wave vector and the angular frequency. With these replacements, the dispersion relation for a massive vacuum GW graviton propagating in the z direction reads

$$\begin{equation} {\bf {k_{\rm g}}}({\omega _{\rm g}}) = \frac{\big({\omega _{\rm g}^{2}}-{\omega _{\rm cut}^2}\big)^{\frac{1}{2}}}{c}\, {\hat{\rm \bf e}_{z}\,\!}\,, \end{equation} \tag{ 2 }$$

where ${\hat{\rm \bf e}_{z}\,\!}$ is the unit vector in the z direction. If the GW frequency ω_g is less than the cutoff frequency ω_cut ≡ m_gc²/ℏ, then the wave vector becomes imaginary, indicating that the wave attenuates and does not propagate. (The equivalent phenomena for electromagnetic waves can be found in Section 87 of Landau & Lifshitz 1960.)

At a spacetime point (t, r), the spatial metric perturbation due to a monochromatic GW is

$$\begin{equation} h_{ab}(t,{\bf r})=\Re \left[\sum _{P=+,\times } A^{P} {\epsilon }^{P}_{ab} e^{i[{\omega _{\rm g}}t-{\bf r} \cdot {\bf {k_{\rm g}}}({\omega _{\rm g}}) ] }\right], \end{equation} \tag{ 3 }$$

where ℜ indicates the real part, and where the a, b range over spacetime indices from 0 to 3. The summation is performed over the polarizations of the GW. Since we are not assuming that general relativity is the theory of gravitation, we could, in principle, have as many as six polarization states. For definiteness, however, and to most clearly show how pulsar timing probes graviton mass, we will confine ourselves in this paper to only the two standard polarization modes of general relativity, denoted + and ×, the usual "TT" gauge (see Appendix A for the details). Thus, the polarization index takes on only the values P = +, ×, with A^P and ^P standing for the amplitude and polarization tensors for the two transverse traceless modes.

The polarization tensor ^P is described in terms of an orthonormal three-dimensional frame associated with the GW propagating direction. Let the unit vector in the direction of GW propagation be ${\hat{\rm \bf e}_{z}\,\!}$; we can choose the other two mutually orthogonal unit vectors ${\hat{\rm \bf e}_{x}\,\!}, {\hat{\rm \bf e}_{y}\,\!}$ to be both perpendicular to ${\hat{\rm \bf e}_{z}\,\!}$. In terms of these three vectors, ${\hat{\rm \bf e}_{z}\,\!}, {\hat{\rm \bf e}_{x}\,\!}$, and ${\hat{\rm \bf e}_{y}\,\!}$, the polarization tensors are given as

$$\begin{eqnarray} &&{\epsilon }^{+}_{ab}={\hat{\rm \bf e}_{x}\,\!}_{a} {\hat{\rm \bf e}_{x}\,\!}_{b} -{\hat{\rm \bf e}_{y}\,\!}_{a}{\hat{\rm \bf e}_{y}\,\!}_{b}, \nonumber \\ &&{\epsilon }^{+}_{ab}={\hat{\rm \bf e}_{x}\,\!}_{a} {\hat{\rm \bf e}_{y}\,\!}_{b} +{\hat{\rm \bf e}_{y}\,\!}_{a}{\hat{\rm \bf e}_{x}\,\!}_{b}\,. \end{eqnarray} \tag{ 4 }$$

Since the polarization tensors are purely spatial, we will have only spatial components of the metric perturbations. For a stochastic GW background, these metric perturbations are a superposition of monochromatic GWs with random phase and amplitude and can be written as

$$\begin{equation} h_{ij}(t,r^i)=\sum _{P=+,\times }\!\!\int _{-\infty }^{\infty }\!\!df_{\rm g}\!\!\int d\Omega \, h^{P}(f_{\rm g},{\hat{\rm \bf e}_{z}\,\!})\, {\epsilon }^{P}_{ij}({\hat{\rm \bf e}_{z}\,\!}) e^{ i[{\omega _{\rm g}}t-{\bf k}_g({\omega _{\rm g}})\cdot {\bf r } ]}, \end{equation} \tag{ 5 }$$

where f_g = ω_g/2π is the GW frequency, Ω is solid angle, spatial indices i, j run from 1 to 3, and h^P is the amplitude of the GW propagating in the direction of ${\hat{\rm \bf e}_{z}\,\!}$ per unit solid angle, per unit frequency interval, in polarization state P. If the GW background is isotropic, stationary, and independently polarized, we can define the characteristic strain h^P_c according to Maggiore (2000) and Lee et al. (2008), and can write

$$\begin{equation} \langle h^{P}(f_{\rm g},{\hat{\rm \bf e}_{z}\,\!}) h^{\star P^{\prime }}(f_{\rm g}^{\prime }, {\hat{\rm \bf e}_{z}\,\!}^{\prime })\rangle = \frac{\big|h_{c}^{P}\big|^2}{16 \pi f_{\rm g}} \delta _{PP^{\prime }} \delta (f_{\rm g}-f_{\rm g}^{\prime })\delta ({\hat{\rm \bf e}_{z}\,\!}-{\hat{\rm \bf e}_{z}^{\prime }}), \end{equation} \tag{ 6 }$$

where the ⋆ stands for the complex conjugate and 〈〉 is the statistical ensemble average. The symbol $\delta _{PP^{\prime }}$ is the Kronecker delta for polarization states; $\delta _{PP^{\prime }}=0$ when P and P' are different, and $\delta _{PP^{\prime }}=1$ when P and P' are the same. With the relationships above, one can show that

$$\begin{equation} \langle h_{ab}(t)h_{ab}(t)\rangle =\sum _{P=+,\times }\int _{0}^{\infty }\frac{\big|h_{c}^{P}\big|^{2}}{f_{\rm g}} df_{\rm g}. \end{equation} \tag{ 7 }$$

We now turn to the calculation of the pulsar timing effects due to the stochastic GW background prescribed by Equation (5). That background is composed of a monochromatic plane wave with random phase and with amplitude determined from some chosen spectrum. We can then first calculate the pulsar timing effect due to a monochromatic plane wave and then add together all contributions to find the effects of a stochastic background.

In the weak field case, the time component of the null geodesic equations for photons is (Liang 2000; Hobbs et al. 2009a)

$$\begin{equation} \frac{d \omega }{d \lambda }+\frac{\omega ^2}{2} \frac{\partial {h_{ij}} }{\partial t} \hat{\rm \bf n}^{i}\hat{\rm \bf n}^{j}=0\,. \end{equation} \tag{ 8 }$$

Here, ω is the angular frequency for the photon, $\hat{\rm \bf n}^{i}$ is the direction from the observer (Earth) to the photon source (pulsar), and λ is an affine parameter along the photon trajectory, normalized so that in the Minkowski background dt/dλ = ω. From this geodesic equation, the frequency shift of a pulsar timing signal induced by a monochromatic, plane, massive graviton GW is

$$\begin{eqnarray} \frac{\Delta \omega (t)}{\omega }&=& -\,\frac{ \hat{\rm \bf n}^{i}\hat{\rm \bf n}^{j} }{2(1+({c}/{{\omega _{\rm g}}}){\rm \bf {k_{\rm g}}} {\,\bm \cdot\,} \hat{\rm \bf n})}\nonumber\\ &&\times[ h_{ij}(t,0)-h_{ij}(t- |{\rm \bf D}|/c,{\rm \bf D})], \end{eqnarray} \tag{ 9 }$$

where it is assumed that the observer is at the coordinate origin and that D is the displacement vector, in the background, from the observer to the pulsar. Equation (9) is a minor generalization of the relationship for zero-mass gravitons given in many references (Estabrook & Wahlquist 1975; Sazhin 1978; Detweiler 1979; Lee et al. 2008). It is clear that the frequency shift of the pulsar timing signal only involves the metric perturbations at the observer, i.e., the term h_ij(t, 0), and the metric perturbation at the pulsar, i.e., the term h_ij(t, D). The dispersion relation of the GWs enters, through the denominator, in the geometric factor $(1+c {\rm \bf {k_{\rm g}}}{\,\bm \cdot\,} \hat{\rm \bf n}/{\omega _{\rm g}})$ as well as the phase difference between the GW at the earth h_ij(t, 0) and the GW at the pulsar h_ij(t − |D|/c, D). The induced pulsar timing residuals R(t) are given by the temporal integration of the above frequency shift at Earth, thus

$$\begin{equation} R(t)=\int _{0}^{t} \frac{\Delta \omega (\tau)}{\omega } d\tau. \end{equation} \tag{ 10 }$$

Equations (5), (9), and (10) determine the response of pulsar timing to a monochromatic plane GW. One can show the pulsar timing residual R(t) induced by the stochastic GW background is

$$\begin{eqnarray} R(t) &=&\frac{i}{2}\sum _{P=+,\times }\int _{-\infty }^{\infty } df_{\rm g} \int d\Omega \frac{ \hat{\rm \bf n}^{i}\hat{\rm \bf n}^{j} }{({\omega _{\rm g}}+c {\rm \bf {k_{\rm g}}} {\,\bm \cdot\,} \hat{\rm \bf n})} h^{P}(f_{\rm g},{\hat{\rm \bf e}_{z}\,\!}){\epsilon }^{P}_{ij}({\hat{\rm \bf e}_{z}\,\!})\nonumber\\ &&\times [1- e^{-i {\rm \bf {k_{\rm g}}({\omega _{\rm g}}) \cdot D}}][1-e^{i{\omega _{\rm g}}t}]. \end{eqnarray} \tag{ 11 }$$

This is, for example, a minor modification of Equation (A9) given by Lee et al. (2008). As we pointed out shortly before, the massive graviton dispersion relation enters via the term ${\omega _{\rm g}}+c {\rm \bf {k_{\rm g}}} {\,\bm \cdot\,} \hat{\rm \bf n}$ as in Equation (9) and the term $1- e^{-i {\rm \bf {k_{\rm g}}({\omega _{\rm g}}) \cdot D}}$, which comes from the phase difference between pulsar term and earth term.

After the nonobservable zero frequency component is removed, the autocorrelation function 〈R(t)R(t + τ)〉 can be calculated by replacing the statistical ensemble average with the time average,

$$\begin{eqnarray} \langle R(t)R(t+\tau)\rangle &=&\!\!\!\! \sum _{P=+,\times } \int _{-\infty }^{\infty } df_{\rm g} \int d\Omega \, \frac{\big|h_{c}^{P}\big|^2}{32 \pi f_{\rm g} } \left[\frac{{\epsilon }^{P}_{ij}\hat{\rm \bf n}^i \hat{\rm \bf n}^j }{{\omega _{\rm g}}+c {\rm \bf {k_{\rm g}}{\,\bm \cdot\,} \hat{\rm \bf n}}}\right]^2\nonumber\\ &\times &\left[1-\cos \left(\frac{D}{c} ({\omega _{\rm g}}+c{\rm \bf {k_{\rm g}}{\,\bm \cdot\,} \hat{\rm \bf n}})\right) \right] e^{i{\omega _{\rm g}}\tau }{,} \end{eqnarray} \tag{ 12 }$$

where D = |D| is the pulsar distance. The folded timing residual power spectra S_R(f_g) are defined to be the Fourier transform of the autocorrelation function of the timing residual,

$$\begin{equation} S_{\rm R}(f_{\rm g})=2\int _{-\infty }^{\infty } \langle R(t)R(t+\tau)\rangle e^{-2\pi i f_{\rm g} \tau } d\tau, \end{equation} \tag{ 13 }$$

for which the explicit result is

$$\begin{eqnarray} S_{\rm R}(f_{\rm g}) &=& \sum _{P=+,\times } \int d\Omega \, \frac{\big|h_{c}^{P}\big|^2}{16 \pi f_{\rm g} } \left[\frac{{\epsilon }^{P}_{ij}\hat{\rm \bf n}^i \hat{\rm \bf n}^j }{{\omega _{\rm g}}+c {\rm \bf {k_{\rm g}}{\,\bm \cdot\,} \hat{\rm \bf n}}}\right]^2\nonumber\\ &&\times\left[1-\cos \left(\frac{D}{c} ({\omega _{\rm g}}+c{\rm \bf {k_{\rm g}}{\,\bm \cdot\,} \hat{\rm \bf n}})\right) \right]. \end{eqnarray} \tag{ 14 }$$

For a power-law stochastic GW background, the power spectra of the induced pulsar timing residual are (see Appendix A for the details)

$$\begin{equation} S_{\rm R}(f_{\rm g})=\sum _{P=+,\times }\frac{\big|h_{c}^{P}(f_{\rm g})\big|^2}{24 \pi ^2 f_{\rm g}^3} \eta (f_{\rm g}), \end{equation} \tag{ 15 }$$

where the η(f_g) is

$$\begin{equation*} \eta (f_{\rm g})=\left\lbrace \begin{array}{@{}c c}\dsty\frac{-4 \zeta ^3+6 \zeta +3 \left(\zeta ^2-1\right) \log \left(\frac{\zeta +1}{1-\zeta }\right)}{2 \zeta ^5} & \hbox{if $f_{\rm g}> f_{\rm cut}$} \\ 0 &\hbox{if $f_{\rm g}\le f_{\rm cut}$} \end{array} \right. \end{equation*}$$

and

$$\begin{equation} \zeta =\sqrt{1-\frac{\omega _{\rm cut}^2}{{\omega _{\rm g}}^2}}\quad \omega _{\rm cut}\equiv {m_{\rm g}}c^2/\hbar {.} \end{equation} \tag{ 16 }$$

One can see that η(f_g) is the spectral correction factor for massive GW backgrounds, if compared with the case of general relativistic GW background (Lee et al. 2008), where η is equal to unity. To show the frequency dependence of the graviton mass effect, η(f_g) is plotted in Figure 1.

Zoom In Zoom Out Reset image size

When f_g > f_cut, a convenient polynomial approximation of η(f) with maximal error of 1.5% is

$$\begin{eqnarray} \eta (f_{\rm g}) &=& -\!0.7764 \left(\frac{f_{\rm cut}}{f_{\rm g}}\right)^4+1.748 \left(\frac{f_{\rm cut}}{f_{\rm g}}\right)^3-1.001 \left(\frac{f_{\rm cut}}{f_{\rm g}}\right)^2\nonumber\\ &&-\!0.5868\, \frac{f_{\rm cut}}{f_{\rm g}}+1.012\,. \end{eqnarray} \tag{ 17 }$$

A least-squares polynomial fitting technique was used to calculate the coefficients in the above equation. The rms level σ_R of the timing residual power is defined to be σ²_R = ∫^∞₀S_R(f)df. The spectra of GW backgrounds generated by various astrophysical processes are usually summarized as power-law spectra with power index α, i.e., the characteristic strain of GWs is h_c = A_c(f/f₀)^α. For such power-law spectra, the rms level of corresponding pulsar timing residuals is

$$\begin{equation} \sigma _{\rm R}^2\simeq \frac{A_{\rm c}^2 f_{\rm L}^{2\alpha }}{f_{0}^{2 \alpha }}\left(\beta _1 \frac{f_{\rm cut}^4}{f_{\rm L}^{6}} +\beta _2 \frac{f_{\rm cut}^3}{f_{\rm L}^{5}} +\beta _3 \frac{f_{\rm cut}^2}{f_{\rm L}^{4}} +\beta _4 \frac{f_{\rm cut}}{f_{\rm L}^{3}} +\beta _5 \frac{1}{f_{\rm L}^{2}} \right), \end{equation} \tag{ 18 }$$

where the constants β₁...β₅ take the following values:

$$\begin{eqnarray} \beta _1 &=&3.278 \times 10^{-3} (\alpha -3)^{-1}, \nonumber \\ \beta _2 &=&1.476 \times 10^{-2} (5-2\alpha)^{-1}, \nonumber \\ \beta _3 &=&4.226 \times 10^{-3} (\alpha -2)^{-1}, \\ \beta _4 &=&4.955 \times 10^{-3} (2\alpha -3)^{-1}, \nonumber \\ \beta _5 &=&4.272 \times 10^{-3} (1-\alpha)^{-1},\nonumber \end{eqnarray} \tag{ 19 }$$

and f_L = Max[T⁻¹, f_cut] is the larger of the following two frequencies: (1) the frequency cutoff (T⁻¹) due to the finite length time span T of observation and (2) the intrinsic frequency cutoff f_cut = ω_cut/(2π) due to graviton mass.

One can derive an upper limit for the GW velocity using single pulsar timing data, because of the surfing effect (Baskaran et al. 2008b), but it is unlikely that one can use single pulsar timing data to constrain the graviton mass (Baskaran et al. 2008a). Because of the correction factor η(f_g) (see Figure 1), the graviton mass reduces the GW-induced pulsar timing residuals. This prevents us from constraining the graviton mass using the amplitude of single pulsar timing residuals. However, as explained in the next section, the cross-correlation between pulsar timing residuals from different directions will help us in detecting the graviton mass.

A stochastic GW background leaves a correlation between timing residuals of pulsar pairs (Hellings & Downs 1983; Lee et al. 2008). Such a correlation, C(θ), depends on the angular distance θ between two pulsars. It turns out that the graviton mass changes the shape of this correlation function. One can therefore detect a massive graviton by examining the shapes of pulsar timing correlation functions.

As shown in Appendix B, the pulsar timing cross-correlation function for a massive GW background depends on the graviton mass, specific power spectra of GW background, and observation schedule. In this way, an analytical expression for the cross-correlation function would not be possible. We use Monte Carlo simulations in this paper to determine the shape of the correlation function for GW backgrounds with a power-law spectra. In the Monte Carlo simulations for C(θ), we randomly choose pulsars from an isotropic distribution over sky positions. We then hold constant these pulsar positions and calculate the angular separation θ between every pair of pulsars. Next, to simulate the power-law GW background, we generate 10⁴ monochromatic waves, choosing random phase and choosing the amplitude from the power law h_c = A_c(f/f₀)^α, where we take α = −2/3 (Phinney 2001; Jaffe & Backer 2003; Wyithe & Loeb 2003; Enoki et al. 2004; Sesana et al. 2004; Wen et al. 2009). The timing residuals are calculated using Equations (9) and (10). Then, the cross-correlation function, C(θ), between pulsar pairs is calculated. We repeat such processes and average over the angular-dependent correlation function C(θ) until the change in C(θ) is less than 0.1%. The averaged correlation function is then smoothed by fitting an eighth-order Legendre polynomial (see Lee et al. 2008 for the details). This smoothed C(θ) is the correlation function we need. C(θ) is plotted for various parameters in Figure 2. We also check the results by choosing different sets of pulsars to make sure that C(θ) is not sensitive to the details of the random pulsar samples. As C(θ) is the statistical expectation of the correlation function, we call it the theoretical correlation function in contrast with the observed correlation function defined in the next section.

Zoom In Zoom Out Reset image size

As one may expect, the massive graviton has stronger effects for data from long observing periods than from short periods. One can see this by comparing the 5 year and 10 year correlation functions given in Figure 2, where curves corresponding to the same range of graviton mass show considerably greater deviations from the massless case in the 10 year correlation function than in the 5 year one.

As we have explained, the massive graviton reduces the pulsar timing response to GWs through the correction factor η(f_g). In this way, a non-zero graviton mass reduces the pulsar timing array sensitivity for detecting a GW background. It is interesting to know how the sensitivity changes, if the graviton is massive.

The stochastic GW background is detected by comparing the measured cross-correlation function c(θ_l) with the theoretical correlation C(θ) calculated in the previous section. Here, the measured cross-correlation function c(θ_l) is defined by

$$\begin{equation} c(\theta _{m})=\frac{ \sum _{l=0}^{N-1} (R_{a}(t_{l})-\overline{R_{a}(t_{l})})(R_{b}(t_{l})- \overline{R_{b}(t_{l})})}{\sqrt{\sum _{l=0}^{N-1}(R_{a}(t_{l})- \overline{R_{a}(t_{l})})^{2}\sum _{l=0}^{N-1}(R_{b}(t_{l})-\overline{R_{b}(t_{l})})^{2}}}, \end{equation} \tag{ 20 }$$

where R_a(t_l) and R_b(t_l) are the timing residuals of pulsars "a" and "b" at time t_l and N is the number of observations. The $\overline{R_{b}(t_{l})}=\sum _{l=0}^{N-1} R_{b}(t_{l})/N$ and $\overline{R_{a}(t_{l})}=\sum _{l=0}^{N-1} R_{a}(t_{l})/N$. The θ_m is the angle between the direction pointing to pulsar "a" and the direction pointing to pulsar "b." Given N_p pulsars, the index m runs from 1 to the number of pulsar pairs M = (N_p − 1)N_p/2, because the autocorrelations are not used.

Following Jenet et al. (2005), we define

$$\begin{equation} \rho = \frac{ \sum _{m=1}^{M}(C(\theta _{m})-\overline{C})(c(\theta _{m})-\overline{c})}{\sqrt{\sum _{m=1}^M (C(\theta _{m})-\overline{C})^2\sum _{m=1}^M(c(\theta _{m})-\overline{c})^2}}, \end{equation} \tag{ 21 }$$

where $\overline{C}=\sum _{m=1}^{M}C(\theta _{m})/M$ and $\overline{c}=\sum _{m=1}^{M}c(\theta _{m})/M$. Then the statistic S, describing the significance of the detection, is $S=\sqrt{M\;}\rho$. In particular, when there is no GW present, c(θ_m) will be Gaussian-like white noise, the probability of getting a detection significance larger than S is about ${\rm erfc}(S/\sqrt{2})/2$ (Jenet et al. 2005).

Our aim is to determine the ability of a given pulsar timing array configuration to detect a GW background. To do this, we calculate the expected value for the detection significance S by using a second set of Monte Carlo simulations. These second Monte Carlo simulations are similar to the first ones, but instead of calculating the average value for C(θ), we inject white noise for each pulsar, to represent the intrinsic pulsar noise and instrumental noise, and we calculate the expected value of S. We summarize the following steps here.

• 1.  
  Generate a large number of GW sources (10⁴) to simulate the required GW background.
• 2.  
  Calculate the timing residual for each pulsar as described above and add white Gaussian noise.
• 3.  
  Calculate the measured correlation c(θ_m) using Equation (20) and calculate the detection significance S using Equation (21).
• 4.  
  Repeat steps 1–3 and average over the detection significance S. The converged S is the value needed to estimate the detection significance.

The results for the expectation value of S, as a function of GW amplitude A_c for various pulsar timing array configurations, are presented in Figure 3. We have also compared simulations from several different pulsar samples with the same number of pulsars to make sure such S is not sensitive to the detailed configuration of the pulsar samples.

Zoom In Zoom Out Reset image size

Two features of the curves in Figure 3 are worth noting. First, the minimal detection amplitude of a GW background becomes larger, when a massive graviton is present, i.e., the leading edge of the S–A_c curve shifts rightwards as m_g is made larger. This tells us that in order to detect a massive GW background, one needs a stronger GW background signal or a smaller pulsar intrinsic noise than in the case of a massless GW background. As previously noted, this effect is mainly due to the reduction of the pulsar timing response and the reduction of the GW amplitude at lower frequencies. Figure 3 also tells us when we can neglect the effect of a massive graviton. It is clear from Figure 3 that if m_g ⩽ 2 × 10⁻²³ eV for a 5 year observation, the minimal detection amplitude is not reduced by more than 5%. For 10 years of observation, a 5% reduction corresponds to m_g = 10⁻²³ eV.

The second noteworthy feature of the S–A_c curves in Figure 3 is that of the saturation level of detection significance. Due to the pulsar distance term of Equation (11) (the term involving the D), the detection significance achieves a saturation level when the GW-induced timing residuals are much stronger than the intrinsic pulsar timing noise (Jenet et al. 2005). From Figure 3, we note that the saturation level of detection significance is large, when the graviton is massive; i.e., the plateau at the right part of the S–A_c curve becomes higher for a more massive graviton. This is rather similar to the whitening filter discussed by Jenet et al. (2005). The graviton mass introduces a low frequency spectral cutoff, which is equivalent to applying a whitening filter to the timing signals. In this way, the saturation level of S starts to grow for a massive GW background. Because the cutoff in the frequency domain coherently removes low frequency GW components, the envelope of the these S–A_c curves is similar to the curves with the whitening filter as described by Jenet et al. (2005).

In Section  5, we have discussed how the detection significance for the stochastic GW background is affected by the graviton mass. Here, we formulate the question as a detection problem rather than a parameter estimation problem. In this section, we will describe a detector that accepts pulsar timing array data and then determines if the GW background is massless or massive. Then, we simulate pulsar timing data, passing through the detector, to determine the quality of such a detector. With the detector quality, we then discuss related technical requirements (the number of pulsars, the pulsar noise level, and so on) for performing such detections.

The question, whether the graviton mass is zero or not, is best formulated as a statistical hypothesis test composed of two hypotheses F and H given by

$$\begin{equation*} \left\lbrace \begin{array}{@{}l} \hbox{H: The graviton is massless,} \\ \hbox{F: The graviton is massive.} \end{array} \right. \end{equation*}$$

The detection algorithm determines which hypothesis is accepted. We define the detection rate P_d as the probability that the detection algorithm gives statement F, when the graviton is massive, i.e., P_d = P(F|m_g > 0); and we define the false alarm rate P_f as the probability of getting statement F when the graviton is massless, i.e., P_f = P(F|m_g = 0). The quality of the detection algorithm is evaluated by calculating the relation between the detection rate P_d and the false alarm rate P_f. Here, we take the standard approach (DiFranco & Rubin 1968) that one fixes P_f to a certain threshold level P_th and calculates the corresponding P_d. Throughout this paper, we fix P_th = 0.1%. If a detector maximizes the P_d for a prescribed value of P_f, we say that such detector is optimal (Kassam 1988).

As we have shown in Section  4, the graviton mass changes the shape of the correlation function. The best way to detect such a difference is to use the following statistics (Kassam 1988):

$$\begin{eqnarray} \gamma &=& \frac{ \sum _{m=1}^{M}(C_{\rm m}(\theta _{m})-\overline{C_{\rm m}})(c(\theta _{m})-\overline{c})}{\sqrt{\sum _{m=1}^M (C_{\rm m}(\theta _{m})-\overline{C_{\rm m}})^2\sum _{j=1}^M(c(\theta _{m})-\overline{c})^2}}\nonumber\\ &&-\frac{ \sum _{m=1}^{M}(C_0(\theta _{m})-\overline{C_0})(c(\theta _{m})-\overline{c})}{\sqrt{\sum _{j=1}^M (C_0(\theta _{m})-\overline{C_0})^2\sum _{j=1}^M(c(\theta _{m})-\overline{c})^2}},\qquad \end{eqnarray} \tag{ 22 }$$

where C₀(θ) is the correlation function for a massless GW background and C_m(θ) is the correlation function for a massive GW background, which maximizes the detection significance S among possible theoretical correlation functions for all possible values of m_g. One can show that the optimal statistical decision rule is (Kassam 1988)

$$\begin{equation*} \left\lbrace \begin{array}{@{}l} \hbox{Choose H: if } \gamma \le \gamma _{\rm th}, \\ \hbox{Choose F: if } \gamma < \gamma _{\rm th}, \end{array} \right. \end{equation*}$$

where the γ_th is the threshold of statistical decision, which is determined to guarantee the false alarm rate constraint P_f ⩽ P_th = 0.1%.

We determine the threshold γ_th by a third set of Monte Carlo simulations. These simulations take the following steps. First, a massless GW background is generated, then we search for the matched value of m_g, such that the corresponding correlation function C_m(θ) maximizes the detection significance S. Then, we use this C_m(θ) to calculate the statistic γ. We repeat this, recording the values of γ, to establish the statistical distribution of γ values. We take the value of the threshold γ_th to be that for which there is less than a 0.1% chance of getting γ>γ_th for the case of m_g = 0.

After we determined the γ_th, a fourth Monte Carlo simulation is used to calculate P_d. This simulation is very similar to the previous one, except that a massive GW background is generated. We repeat the simulation 10³ times and take as P_d the probability of getting γ>γ_th.

We summarize the results of the simulation in Figure 4, which are gray-scale contour plots for the detection rate P_d for different scenarios of using 60, 100, and 300 pulsars, respectively. The corresponding parameters are given in the legend of each panel. Intuitively, the necessary conditions for a positive detection of a graviton mass should be, first, that the GW is strong enough for the GW to be detected, and second, that the graviton mass is large enough to change the shape of the correlation function. This intuition is confirmed by our simulations, which show that the high detection rate concentrates in the upper right corner of each panel, where both graviton mass and GW amplitude are large enough.

Zoom In Zoom Out Reset image size

From Figure 4, we can also see that we need at least 60 pulsars to be able to tell the difference between a massive GW background and a massless one. For 5 year observations of 100 pulsars, we can start to detect a graviton heavier than 2.5 × 10⁻²² eV and can achieve a limit of m_g = 10⁻²² eV by using 5 year observations of 300 pulsars. We can achieve levels of 10⁻²² eV and 5 × 10⁻²³ eV in 10 year observations using 100 and 300 pulsars, respectively.

We also note that there is a positive correlation between the minimal detectable graviton mass and the GW background amplitude, as shown by the leftwards lead edge of the contours, in other words, we need a larger GW amplitude such that the GW background can be detected, if gravitons are more massive. As discussed above, this correlation is due to the reduction of pulsar timing residuals due to massive graviton.

Current efforts to detect GWs using radio pulsar observations make use of correlated arrival time fluctuations induced by the presence of a stochastic GW background. Einstein's general theory of relativity makes a specific prediction for the angular dependence of the correlation. The power spectrum of the induced fluctuations is given by the power spectrum of the GW background, determined by the physical processes of the generation process, divided by the square of the GW frequency. In this paper, we showed that the form of the expected correlation function, as well as the power spectrum of the induced timing fluctuations, will be significantly altered if the graviton had a non-zero mass. Only the transverse traceless GW modes were considered in this analysis.

A non-zero graviton mass introduces a cutoff frequency below which no GWs will propagate. If we consider a fixed amount of time over which the correlation will be measured, 5 years for example, a small increase in the graviton mass will actually increase our ability to detect a strong background. Here, we are assuming that the background is being generated by an ensemble of supermassive black hole binaries. This effect is due to the fact that the presence of the graviton mass will act to flatten out the spectrum of the induced time residuals (see Equation (15) and Figure 1), thus making the strong background easier to detect. As the mass is increased, the cutoff frequency will increase. Once this cutoff frequency is above the inverse of the observing time, the sensitivity starts to fall off dramatically. This is due to the fact that pulsar timing is most sensitive to the lowest observable frequencies and the presence of a massive graviton removes all GW power at frequencies below the cutoff frequency.

A non-zero graviton mass will also change the shape of the angular-dependent correlation function. In order to understand this, consider two GWs traveling near the direction of the line of sight between the Earth and a given pulsar. One GW is traveling toward the pulsar and the other is traveling toward the Earth. The induced timing fluctuations will be very different for each of these waves. The GW traveling toward the Earth will induce higher amplitude fluctuations then its counter part traveling toward the pulsar. Now, consider the same case but with a massive graviton. As the GW frequency approaches the cutoff frequency, the GW wavelength gets arbitrarily large. Assume that the GW wavelength is much larger than the Earth–pulsar distance. In this case, both GWs look exactly the same as far as the Earth–pulsar detector is concerned. Hence, both waves induce the same timing fluctuations. This symmetry between the two GW directions implies that the timing fluctuations from two pulsars near each other on the sky or 180° apart will have the same amount of correlation, unlike the massless graviton case. Thus, as the graviton mass increases, the correlation curve will become more symmetric about 90°.

Since the graviton mass changes the expected correlation function, we can detect the massive graviton by measuring the shape of correlation function. The optimal detection algorithm differentiating between a massive and a massless GW background was constructed in this paper. Using this algorithm, we found that at least 60 pulsars are required to discriminate between a GW background made up of massive and massless gravitons. Note that, like the case of discriminating various polarization modes of GWs (Lee et al. 2008), the minimum number of pulsars is insensitive to the intrinsic timing noise level or the total observing time.

Here, we are mainly focusing on answering the following question: what is the threshold of graviton mass such that the detector, which we described in this paper, will find out that the GW background is composed of massive gravitons rather than massless gravitons? We answered this statistical detection question in this paper. The related parameter estimation question⁴ will be investigated in the future works.

In this paper, we treat the graviton mass in the sense of a GW dispersion relation. A better graviton mass upper limit (m_g ⩽ 2 × 10⁻²⁶ eV) can be achieved by observing the GW dispersion using Laser Interferometer Space Antenna (LISA; Stavridis & Will 2009). Besides GW dispersion-based techniques, there are other good upper limits (2 × 10⁻²⁹h₀ eV) from galaxy cluster observations (Goldhaber & Nieto 1974). However, since some gravity theories (Maggiore 2008) contain mass terms for h_ab, but maintain the scalar sectors, the m_g upper limits from GW dispersion observations are independent of the upper limits from Yukawa potential experiments such as the solar system and the galaxy cluster observations.

In summary, for the task of detecting the massive graviton using a pulsar timing array, there is one critical requirement: a large sample of stable pulsars. Thus, the ongoing and upcoming projects like the Parkes Pulsar Timing Array (Hobbs et al. 2009b), European Pulsar Timing Array (Stappers et al. 2006), NANOGrav (Jenet et al. 2009), the Large European Array for Pulsars (Stappers et al. 2009), the Five-hundred-meter Aperture Spherical Radio Telescope (Nan et al. 2006; Smits et al. 2009), and the Square Kilometer Array will offer unique opportunities to find and to time these pulsars, and to detect the GW background and measure its properties.

K.J.L. gratefully acknowledges support from ERC Advanced Grant "LEAP," Grant Agreement Number 227947 (PI: M. Kramer), F.A.J. and R.H.P. gratefully acknowledge support from the NSF under grants AST0545837 and PHY0554367 and support from the UTB Center for Gravitational Wave Astronomy. We also thank M. Pshirkov and N. Straumann for very useful comments.

The power spectra of timing residuals are calculated by integrating Equation (14). Here, we give the details of the calculation.

In Figure 5, the components of ${\hat{\rm \bf e}_{z}\,\!}$ in the $\hat{\bf X}^a=(\hat{\bf X},\hat{\bf Y},\hat{\bf Z})$ frame can be seen to be (sin θ_gcos ϕ_g,sin θ_gsin ϕ_g, cos θ_g), where the θ_g and ϕ_g are, respectively, the polar angle and azimuthal angle of the GW propagation vector. To proceed, we need the components of the polarization tensors in the $\hat{\bf X}^a$ frame. The transformation from the components _0,cd given in the GW frame $\hat{\bf e}_b =(\hat{\bf e}_x,\hat{\bf e}_y,\hat{\bf e}_z)$ of Equation (4) is made with ^P_ab = T_caT_db ^P_0,cd, where $T_{ca}=\hat{\bf e}_c{\,\bm \cdot\,} \hat{\bf X}_a$ has components

$$\begin{equation} \hbox{\fontsize{5.9}{8}\selectfont $\left({ \begin{array}{@{}ccc@{}} \cos \theta _{g}\cos \phi _{g}\cos {\psi _{g}}-\sin \phi _{g}\sin {\psi _{g}} &\cos \theta _{g}\cos \psi _{g}\sin \phi _g+\cos \phi _g\sin \psi _g &-\cos \psi _{g}\sin \theta _{g} \\ -\cos \psi _{g}\sin \phi _{g}-\cos \theta _{g}\cos \phi _{g}\sin \psi _{g} &\cos \phi _g\cos \psi _g-\cos \theta _{g}\sin \phi _{g}\sin \psi _g & \sin \theta _g\sin \psi _g \\ \cos \phi _{g}\sin \theta _{g} &\sin \theta _g\sin \phi _g &\cos \theta _{g} \end{array} }\right)$}. \end{equation} \tag{ A1 }$$

Zoom In Zoom Out Reset image size

Since the GW background is isotropic, with no loss of generality one can choose $\hat{\rm \bf n}= \lbrace 0, 0, 1\rbrace$ so that

$$\begin{equation} {\omega _{\rm g}}+c\hat{\rm \bf n}{\,\bm \cdot\,} {\rm \bf {k_{\rm g}}} = (1+\zeta \cos \theta _{\rm g})\,{\omega _{\rm g}}, \end{equation} \tag{ A2 }$$

where

$$\begin{equation} \zeta =\sqrt{1-\frac{{m_{\rm g}}^2 c^4}{{\omega _{\rm g}}^2 \hbar ^2}\;}=\sqrt{1-\frac{f_{\rm cut}^2}{f_{\rm g}^2}\;}\,, \end{equation} \tag{ A3 }$$

for f_g ⩾ f_cut . The term ${\epsilon }^{P}_{ij}\hat{\rm \bf n}^i \hat{\rm \bf n}^j$ now simplifies to

$$\begin{eqnarray} {\epsilon }^{+}_{ij}\hat{\rm \bf n}^i \hat{\rm \bf n}^j &=&\sin ^2\theta _{\rm g} \cos (2\psi _{\rm g}) \nonumber \\ {\epsilon }^{\times }_{ij}\hat{\rm \bf n}^i \hat{\rm \bf n}^j &=&-\!\sin ^2\theta _{\rm g} \sin (2\psi _{\rm g}). \end{eqnarray} \tag{ A4 }$$

After the ensemble average over the polarization angle ψ_g, the integration of Equation (14) gives

$$\begin{eqnarray} S_{\rm R}(f_{\rm g})&=&\sum _{P=+,\times }\frac{\big|h_{c}^{P}(f_{\rm g})\big|^2}{128 \pi ^3 f_{\rm g}^3} \int _{0}^{\pi } d\theta _{\rm g} \int _{0}^{2\pi } d\phi _{\rm g} \frac{\sin ^5\theta _{\rm g}}{(1+\zeta \cos \theta _{\rm g})^2}\nonumber\\ &&\times\left[1-\cos \left(\frac{D{\omega _{\rm g}}}{c} (1+\zeta \cos \theta _{\rm g})\right)\right]. \end{eqnarray} \tag{ A5 }$$

For practical pulsar timing array, we have Φ₀ = Dω_g/c ≫ 1, i.e., the GW wavelength is much smaller than the pulsar–earth distance, so one gets

$$\begin{eqnarray} &&\int _{-1}^{1} \frac{(1-\mu ^2)^2}{(1+\zeta \mu)^2}\left[1-\cos \left(\Phi _0 (1+\zeta \mu)\right)\right] d\mu\nonumber\\ &&\quad = \frac{4 \left[-4 \zeta ^3+6 \zeta +3 (\zeta ^2-1) \log \left(\frac{\zeta +1}{1-\zeta }\right)\right]}{3 \zeta ^5} + O\left(\Phi _0^{-2}\right).\qquad \end{eqnarray} \tag{ A6 }$$

Finally, the power spectra of the pulsar timing residuals then become

$$\begin{equation} S_{\rm R}(f_{\rm g})=\sum _{P=+,\times }\frac{\big|h_{c}^{P}(f_{\rm g})\big|^2}{24 \pi ^2 f_{\rm g}^3} \eta (f_{\rm g}). \end{equation} \tag{ A7 }$$

For f_g > f_cut, η is given by

$$\begin{equation} \eta (f_{\rm g})=\frac{-4 \zeta ^3+6 \zeta +3 (\zeta ^2-1) \log \left(\frac{\zeta +1}{1-\zeta }\right)}{2 \zeta ^5}\,, \end{equation} \tag{ A8 }$$

where ζ is given by Equation (A3) . For f_g ⩽ f_cut, the GWs cannot propagate, and we take η = 0.

The correlation function C(θ) between pulsars "a" and "b" defined in this paper is

$$\begin{equation} C(\theta)=\frac{\langle R_{a}(t) R_{b}(t)\rangle }{\sigma _{a} \sigma _{b}}, \end{equation} \tag{ B1 }$$

where θ is the angle between the direction to pulsar "a" and the direction to pulsar "b". The R_a(t) and R_b(t) are the GW-induced timing residuals for pulsars "a" and "b," respectively. The σ_a and σ_b are the rms values for the timing residual R_a(t) and R_b(t). One has (see Lee et al. 2008 for a similar calculation)

$$\begin{eqnarray} C(\theta)&=&\frac{1}{\sigma _{a} \sigma _{b}} \sum _{P=+,\times } \Bigg\langle \int _{0}^{\infty } \frac{\big|h_{c}^{P}(f_{\rm g})\big|^2}{64 \pi f_{\rm g}} df_{\rm g}\nonumber\\ &&\times\int d\Omega \frac{{\epsilon }^{P}_{ij}\hat{\rm \bf n}_{a}^i \hat{\rm \bf n}_{a}^j }{{\omega _{\rm g}}+c {\rm \bf {k_{\rm g}}{\,\bm \cdot\,} \hat{\rm \bf n}_{a}}} \frac{{\epsilon }^{P}_{ij}\hat{\rm \bf n}_{b}^i \hat{\rm \bf n}_{b}^j }{{\omega _{\rm g}}+c {\rm \bf {k_{\rm g}}{\,\bm \cdot\,} \hat{\rm \bf n}_{b}}} {\cal P}_{ab} \Bigg\rangle,\qquad \end{eqnarray} \tag{ B2 }$$

where ${\cal P}_{ab}$ is given by

$$\begin{eqnarray} {\cal P}_{ab}&=& 1- \cos \Phi _{a}-\cos \Phi _{b}+\cos (\Phi _{a}-\Phi _{b}), \nonumber \\ \Phi _{a}&=&\frac{D_{a}}{c} ({\omega _{\rm g}}+c{k_{\rm g}}\cdot \hat{\rm \bf n}_{a}), \nonumber \\ \Phi _{b}&=&\frac{D_{b}}{c} ({\omega _{\rm g}}+c{k_{\rm g}}\cdot \hat{\rm \bf n}_{b}). \end{eqnarray} \tag{ B3 }$$

Equation (B2) shows that the correlation function is composed of two major integrations, one over the GW frequency df_g and the one over solid angle dΩ. For the massless graviton case, these two parts are not mixed, since k_g is linear in ω_g. But k_g is nonlinear in ω_g for the case of a massive graviton. The nonlinear dependence in Equation (2) mixes the spatial integral with the frequency one. Thus, the shape of C(θ) depends on both the graviton mass and the frequency spectra of GWs, which also depends on the observation schedule. Thus, it is unlikely that one can integrate Equation (B2) analytically.