New Graviton Mass Bound from Binary Pulsars

In Einstein's general relativity, gravity is mediated by a massless metric field. The extension of general relativity to consistently include a mass for the graviton has profound implications for gravitation and cosmology. Salient features of various massive gravity theories can be captured by Galileon models, the simplest of which is the cubic Galileon. The presence of the Galileon field leads to additional gravitational radiation in binary pulsars where the Vainshtein mechanism is less suppressed than its fifth-force counterpart, which deserves a detailed confrontation with observations. We prudently choose fourteen well-timed binary pulsars, and from their intrinsic orbital decay rates we put a new bound on the graviton mass, $m_g \lesssim 2 \times 10^{-28}\,{\rm eV}/c^2$ at the 95% confidence level, assuming a flat prior on $\ln m_g$. It is equivalent to a bound on the graviton Compton wavelength $\lambda_g \gtrsim 7 \times 10^{21}\,{\rm m}$. Furthermore, we extensively simulate times of arrival for pulsars in orbit around stellar-mass black holes and the supermassive black hole at the Galactic center, and investigate their prospects in probing the cubic Galileon theory in the near future.


I. INTRODUCTION
The late-time accelerated cosmic expansion poses a profound challenge for modern physics, which is known as the dark energy problem [1][2][3]. From observations, we know that dark energy manifests at lengthscales significantly larger than the Galactic size, or in field-theoretic terminology, in the infrared regime. Dark energy is often hypothesized as a cosmological constant in the standard ΛCDM model [4], but its real nature remains elusive. Astrophysical objects (in particular, the type Ia supernovae) in the relatively nearby Universe and the cosmic microwave background in the early Universe provide two classes of independent probes to measure the Hubble expansion parameter at today, H 0 . Recent observations from them, however, have inferred inconsistent values of H 0 at a significance level of 4.4 σ [5][6][7]. The discrepancy aggravates the dark energy puzzle, and in the meantime has triggered tremendous interest in searching for new physics beyond the standard paradigm.
One of the main approaches to explain the dark energy phenomena involves infrared modifications to the canonical gravity theory, general relativity (GR) [8][9][10]. Modifications usually introduce extra field contents, with a scalar degree of freedom being the simplest and the most widely investigated in literature. However, such a new scalar is likely to bring in a fifth force [8][9][10], which is stringently constrained by observations in the Solar system [11] and binary pulsars [12][13][14][15][16][17]. Therefore, to successfully account for the accelerated expansion of the Universe, we need a modified gravity theory where the theory gives rise to order-one corrections at cosmological scales but deviations from GR are extremely suppressed in the Solar system, which makes nonlinearity a crucial ingredient in the theory. For a class of such infrared modifications of gravity, including the Galileon models [18,19], this is achieved by * lshao@pku.edu.cn the Vainshtein mechanism [20,21], by which the new scalar becomes nonlinearly coupled in the local dense environment, thus suppressing the fifth force in the Solar system. The length scale within which the scalar becomes strongly-coupled is called the Vainshtein radius, r , and it is only outside of the Vainshtein radius that the linear perturbations can be trusted. It is important that when dealing with models with the Vainshtein mechanism the full nonlinear theory, as opposed to the linear theory, needs to be solved to make correct physical interpretations.
The Vainshtein mechanism is intimately related to massive gravity, and the accelerated cosmic expansion may be due to a condensate of gravitons with a Hubble-scale mass [22,23]. It was pointed out in 1970s that the unique (Lorentz invariant) linear theory of massive gravity deviates from GR by orderone corrections, known as the van Dam-Veltman-Zakharov discontinuity [24][25][26]. Vainshtein soon after suggested that this can not be used to rule out massive gravity and rather one needs to solve the nonlinear theory in environments such as the Solar system to get the right prediction [20]. In other words, while the conventional helicity-2 modes of GR become strongly coupled at the Schwarzschild radius, the extra modes of massive gravity becomes strongly coupled within a much larger Vainshtein radius for the same central mass. To extract the most strongly coupled extra modes in massive gravity, one takes the decoupling limit as per the de Rham-Gabadadze-Tolley (dRGT) tuning [27,28], where m g is the graviton mass and M Pl ≡ 1/ √ 8πG is the reduced Planck mass, and obtains a scalar effective field theory with the Galileon symmetry, where π s is the Galileon field and a and b µ are constants. The decoupling limit (1) scales away the small effects due to other modes and let us focus on the physics most important near the scale of Λ. In the presence of matter sources, we also take the limit where the energy momentum tensor T µν goes to infinity but the ratio between T µν /M Pl fixed. The Galileon scalar then couples to the trace of T µν and encapsulates the salient features, including the Vainshtein mechanism, of the extra modes in massive gravity [18,19].
In this paper, we study the simplest model that exhibits the Vainshtein mechanism, namely the cubic Galileon [18]. This model is the decoupling limit of the Dvali-Gabadadze-Porrati braneworld model [18,29], where the graviton acquires a socalled soft mass from embedding the 3-brane Universe in a 4-D bulk with an Einstein-Hilbert term. The Galileon models are also the decoupling limit of the recently discovered dRGT model [27,28], a unique nonlinear (Lorentz invariant) massive gravity with a so-called hard mass. The bi-gravity [30] or multi-gravity extension of the dRGT model also leads to a bi-Galileon or multi-Galileon theory [31][32][33]. Therefore, the cubic Galileon model is often taken as a proxy to cover all the Lorentz invariant massive gravity models, though by no means it encodes all aspects of a complete theory of massive gravity where other terms, for example the quartic Galileon term, might appear.
Horndeski theory [34], the generalized scalar-tensor theory with up to second derivatives in field equations, can be rederived in the Galileon framework [35]. It is worth to mention that, while a large class of Horndeski models have been ruled out by the coincident observation of the gravitational-wave signal [36] and the electromagnetic counterpart [37] from a binary neutron star inspiral GW170817 [38][39][40], 1 the cubic Galileon subset of Horndeski theory, though being a simple strawman model, is still comfortably alive. Consequently, it is intriguing to study the cubic Galileon in the light of recent research activities in the field.
The fifth force effects of massive gravity are often screened by the Vainshtein mechanism so effectively that the existing constraints in the dense environment can be easily evaded [42], and only at a cosmological density in the infrared regime can the theory deviate significantly from GR, accounting for the dark energy. However, de Rham et al. [43] has found that, in binary pulsar systems, the suppression factor in the extra gravitational radiation due to the Galileon mode is less than the suppression factor in the static fifth-force effect (see Sec. II). Therefore, it becomes extremely interesting to check with the existing tests related to the gravitational radiation in binary pulsars for the cubic Galileon model.
In this work, we present a thorough phenomenological study of the Galileon radiation for binary pulsar systems. We carefully choose fourteen well-timed binary pulsars to put constraints on the theory parameter of the cubic Galileon. Among these pulsars, recent observations of the double pulsar PSR J0737−3039A [44,45] give the strongest bound on the graviton mass, m g 3 × 10 −28 eV/c 2 at the 95% confidence level (C.L.). Furthermore, a combination of all fourteen pulsars in the Bayesian framework gives, with a flat prior on ln m g . It translates into a limit on the graviton Compton wavelength λ g 7 × 10 21 m. The paper is organized as follows. In the next section, we review the basics for the Galileon radiation [43]. In Sec. III, systematic studies are carried out to understand the dependence of the Galileon radiation on system parameters and the figure of merit to test it. Based on these studies, we choose fourteen binary pulsars to cast tight constraints on the graviton mass. Limits are obtained from individual pulsars as well as a combination of them. Moreover in Sec. IV, with a set of simulated times of arrival for near-future radio telescopes, we investigate the prospects in using pulsars around a stellarmass black hole (BH) companion [46][47][48] and the supermassive BH at the Galactic center (namely, the Sgr A * ) [49][50][51][52] to constrain the cubic Galileon theory. The last section presents some discussion and briefly summarizes the paper.
Throughout the paper, we implicitly assume units where = c = 1, except for a couple of places where and c are restored for readers' convenience.

II. THEORY
Since the Hulse-Taylor pulsar provided the first indirect evidence for the existence of gravitational waves [53], many more binary pulsars have been playing an important role in probing the property of the gravitational radiation in alternative gravity theories [12,45,[54][55][56][57]. We briefly review the gravitational radiation for a binary pulsar in GR and in cubic Galileon theory in Sec. II A and Sec. II B respectively.

A. General Relativity
Consider a binary system with component masses m 1 and m 2 in an orbit with a semimajor axis a and an eccentricity e. Due to the finite propagating velocity of gravity, at the leading order the binary loses energy by radiating off gravitational waves with an emitting power [58], where the total mass M ≡ m 1 + m 2 ; the symmetric mass ratio η ≡ m 1 m 2 /M 2 ; G and c are the gravitational constant and the speed of light respectively. From the Kepler's third law for a binary system, we have where n b ≡ 2π/P b with P b the orbital period. The nonrelativistic orbital energy at the Newtonian order for the binary reads, Taking time derivatives in Eqs. (5) and (6), we havė Finally, using the energy conservation law in GR, P GR = −Ė b , we have [58],

B. Cubic Galileon
As discussed in the Introduction, (Lorentz invariant) massive gravity models in the decoupling limit essentially reduce to Galileon models plus linearized helicity-2 modes. For the cubic Galileon, we focus on the action [18,43], where h µν ≡ g µν − η µν is the perturbation of the metric, the first two terms in the integrand are the linearized Einstein-Hilbert term coupled to matter with the Lichnerowicz operator (Eh) µν ≡ − 1 2 h µν +· · · , T is the trace of the energy-momentum tensor T µν , and Λ is the strong coupling scale of the Galileon sector, related to the mass of graviton m g via Λ 3 = m 2 g M Pl . Therefore, the field equations for π s and h µν decouple, For a static system whose total mass is M, one can define the Vainshtein radius as, Within r , the fifth force from the scalar degree of freedom is strongly suppressed. In a static system, the fifth force is suppressed by a factor ∼ (L/r ) 3/2 , where L is the typical lengthscale of the system [21]. For example, for an imaginary "static" binary system, we can choose it to be the semimajor axis of the orbit, L ∼ a. However, binaries are not static. de Rham et al. [43] and Chu and Trodden [60] explicitly worked out the gravitational radiation behaviors in a time-dependent binary system at lowest orders. They found that, (i) extra Galileon radiation powers exist at the monopole, dipole, and quadrupole levels, and (ii) for the dominant radiation the suppression factor is weakened from (L/r ) 3/2 to (n b r ) −3/2 . Similar to GR, the Newtonian-order contribution of both the monopole and dipole Galileon radiation vanishes, due to the conservation of energy and linear momentum of the system respectively. Post-Newtonian order contributions do exist for them. For the quadrupole radiation, the Newtonian-order contribution is nonzero. We collect the radiation powers of these multiple moments from Ref. [43] (see also the numerical calculation in Ref. [61]) in the following for later use in this work.
• Monopole radiation. In the cubic Galileon model (10), the monopole radiation power at the leading order for a binary system is, where the constant β and the "monopole mass" M mono (also known as the reduced mass) are defined respectively as, with Γ(·) the gamma function, and the eccentricity mode function I mono n (e) is given in Eq. (24) and Fig. 1.
• Dipole radiation. The dipole radiation power at the leading order for a binary system is, where the constant c 1 and the "dipole mass" are defined respectively as , where X ≡ (m 1 − m 2 ) /M, and the eccentricity mode function I  for the quadrupole (bottom). Notice that, for clarity panels (b) and (d) have used a logarithmic scale for the vertical axes, while all the other axes are using a linear scale. For a direct check, the solid black lines in the right column of the panels reproduce the result in Figure 1 of de Rham et al. [43] for the Hulse-Taylor pulsar PSR B1913+16 whose orbital eccentricity is e 0.617 [59].
where the constant λ and the "quadrupole mass" M quad are defined respectively as, where The above mentioned eccentricity functions for monopole, dipole, and quadrupole radiations are, The behaviors of these functions are illustrated in Fig. 1 for different values of the eccentricity.
Using the same reasoning of energy balance in Sec. II A, we can get the extra contributions toṖ b from the extra Galileon radiation powers in Eqs. (14), (17), and (20), In obtaining these results, at the leading order we have used the Newtonian binding energy for the orbit in Eq. (6), and we have made use of the orbital-averaged radiating powers as discussed above. Such a simplification is sufficient for the analysis in this paper. Notice that the extra contributions tȯ P b are all proportional to the graviton mass m g , while in the linearized Fierz-Pauli theory [48,62,63], the extraṖ b ∝ m 2 g .

III. CONSTRAINTS FROM BINARY PULSARS
Now we would like to better understand the physical effect of the Galileon radiation in Eqs. (27)(28)(29) to different kinds of binary pulsars. We present some general consideration in Sec. III A, concerning the dependence on the orbital eccentricity e, the orbital period P b , and the component masses (m 1 , m 2 ). We cast constraints on the graviton mass in the cubic Galileon theory from precision timing of fourteen carefullychosen binary pulsars in Sec. III B.

A. General consideration
From Fig. 1 we see that, (i) given an eccentricity, the absolute values of the eccentricity mode functions, I n (e), increase at first, then have a peak at a particular n = n peak before they decrease; (ii) with the eccentricity increasing, n peak happens at a larger n, therefore when we deal with highly eccentric binary pulsars, more modes are to be included in order to guarantee the convergence of the sum; and (iii) the cumulative contributions of these mode functions to the Galileon radiation powers in Eqs. (14), (17), and (20) saturate after a particular n, and for e 0.95, summing up to N ∼ 15 is generally sufficient. In our following calculation, we choose N = 30 for a better accuracy. f (e) The f (e) factor inṖ b , defined in Eqs. (30)(31)(32)(33), with a normalization such that f GR = f mono = f dipo = f quad for the double pulsar [44] whose e 0.088 (dotted vertical line).
However, one should remember that, for extremely eccentric binaries with 1 − e 10 −2 , a larger cutoff N is needed.
Based on the contribution to the Galileon radiation power, we define, and These are the dependence shown in the GR radiation (9) and Galileon radiation (27)(28)(29). The proportional factors in Eqs. (31)(32)(33) can be arbitrary normalization values for the convenience of specific demonstration. The f (e) functions in Eqs. (30)(31)(32)(33) are plotted in Fig. 2 with a choice of normalization. We can see that the more eccentric binaries are emitting more gravitational waves and contributing more significantly toṖ b . This is true for the three types of Galileon radiation, as well as the quadrupole radiation in GR. If the binary is extremely eccentric with 1 − e 1, a huge amplificative factor occurs, especially for the dipole radiation (28) in the cubic Galileon. As we will see below, the quadrupole Galileon radiation is in general the main contributor toṖ b for binaries with comparable component masses [43]. The increase of f quad (e) at large eccentricities is slower than the other ones. Notice that the curves in the figure only represent one numerical factor in the radiation power, defined in Eqs. (30)(31)(32)(33), while the other dependences (for example, the dependence on the orbital period and the component masses) are omitted here, which will be investigated in the following. Worth to note that, generally the increased radiation due to eccentricity will tend to circularize the orbits (see e.g. in GR [58]), making highly eccentric binaries less likely to exist to date. It is interesting to observe that, compared withṖ GR b /P b ∝ n  29)]. It is a noteworthy feature for the gravitational radiation with the screening mechanism. With such a dependence, while more relativistic binaries are more prominent in emitting gravitational waves in GR, this may not always be true in the cubic Galileon. It resembles the dependence on n b in Lorentz-violating massive gravity [62,63] and gravity theories with a time-varying gravitational "constant" G(t) [12,48].
In Fig. 3, we plot three types of Galileon radiation, as function of the orbital period, for a binary pulsar system with component masses similar to the Hulse-Taylor pulsar, and a graviton mass m g = 10 −27 eV/c 2 . First, we observe that, the dipole radiation is orders of magnitude smaller than the other two types of Galileon radiation, thus it can be totally ignored. Though, as we recall from Fig. 2, the dependence of the f (e) factor for the dipole radiation on the eccentricity e is one of the steepest, the overall effect from dipole radiation is negligible even for e = 0.98. Then, we observe in Fig. 3 that, while for binaries with e 0.6, the quadrupole radiation dom- inates [43], for highly eccentric binaries with relativistic orbits (namely, smaller P b ), the monopole radiation might dominate instead. This happens when P b 10 3 s for e = 0.9, and P b 10 5 s for e = 0.98. It is caused by the facts thatṖ mono b increases for larger n b (smaller P b ), whileṖ quad b increases for smaller n b (larger P b ). Therefore, in particular for highly eccentric binaries to be discovered in the future, we should keep in caution of only considering the quadrupole Galileon radiation. This also applies to highly eccentric Galactic binaries (if they exist) in the mHz gravitational-wave band for the Laser Interferometer Space Antenna (LISA) [64,65]. A similar discussion will also be mentioned for pulsar-BH systems in Sec. IV. In our following calculation, we keep all three types of Galileon radiation summed, no matter of their relative strength.
In Fig. 4, we plot the relative contribution of the Galileon radiation to the quadrupole radiation in GR. First of all, we see that the relative contribution increases when the orbital period is larger. This is in accordance with the screening mechanism which works less well when the size increases [21], while the GR effects are more prominent when the orbit is more compact. Specifically it is caused by that, (i) in this regime, the quadrupole Galileon radiation dominates, which increases as ∝ P 1/6 b when P b becomes larger; (ii) the quadrupole radiation in GR decreases as ∝ P −5/3 b . With m g = 10 −27 eV/c 2 the Galileon radiation may even be larger than the GR one when P b day for nearly circular orbits. Moreover, we find that, while the absolute Galileon quadrupole radiation increases with larger eccentricity (see Fig. 2), the relative contribution to GR decreases. It can be understood from the steepness of the curves in Fig. 2 for the quadrupole radiations in GR (black dots) and in the cubic Galileon (magenta circles). The latter one is less steep, thus decreasing the relative ratio ofṖ b with larger eccentricity.
Lastly, to understand the effect of component masses, in Fig. 5 we plot the relative strength toṖ b for different mass pairs, normalized to (m 1 , m 2 ) = (1.438, 1.390) M , which is the fiducial mass pair based on the Hulse-Taylor pulsar [59]. This figure is useful to analyze binaries with given orbital period and eccentricity, but different component masses. It is well known that, for the quadrupole radiation in GR, the so-called chirp mass, M ≡ (m 1 m 2 ) 3/5 /M 1/5 , plays the key role, whose contours in the m 1 -m 2 parameter space are shown in the upper left panel in Fig. 5. The monopole and quadrupole Galileon radiations depend on the component masses in a quantitatively different, but qualitatively similar way, as shown in the right column of panels. But for the dipole radiation, shown in the bottom left panel, the dependence is totally different. Asymmetric binaries are preferred to manifest the dipole radiation. Notice that from the colorbars, the enhancement from asymmetric component masses for the dipole radiation can be as large as ∼ 10 3 , while for the other three types of radiation, the change with component masses is within a relatively limited range. This point will be clearly observed for pulsar-BH systems in Sec. IV.

B. Constraints
After having looked at the general features for the Galileon radiation in the last subsection, now we turn to realistic binary pulsar systems. We carefully choose fourteen binary pulsars that are known to be precisely timed. They are chosen from the pool of millisecond pulsars in the second Data Release of the International Pulsar Timing Array program [70], as well as several other well-known systems with measurement oḟ P b [44,45,54,59,68,69,[71][72][73][74][75][76][77][78]. The relevant parameters for the test of the Galileon radiation are listed in Table I. Because binary pulsars have a large variety in terms of system parameters and observational characteristics, for detailed description of these systems, readers are referred to the original timing papers which are given for each pulsar in the table. As one can see from the last subsection, the Galileon radiation depends on a set of physical parameters of the binary pulsar system. These pulsars in Table I represent, to our best knowledge and resource, a set of the most suitable binary pulsar systems to date to perform the test in this paper.
For the chosen binary pulsars, post-Keplerian parameters, other thanṖ b , were used to calculate the component masses m 1 and m 2 , assuming the validity of GR expressions [12,67]. The component masses are listed in the fourth and fifth columns in Table I. In the strictest sense, the component masses should be calculated consistently using cubic Galileon theory instead of GR. Nevertheless, as the Vainshtein suppression of the fifth force is more significant than the Galileon radiation [43,79], it is safe to use the GR formulae to extract these parameters.
As we can see in the penultimate column in the table, most of the chosen pulsars have uncertainties for the value of measuredṖ b at the level of σṖobs b ∼ 10 −15 s s −1 . PSR B1534+12 has an even better measurement with σṖobs b 3×10 −16 s s −1 [75]. The excellent timing precision is attributed to the long-term observation of these pulsars and the continuous improvements of the instruments at large radio telescopes [67,80,81]. However, these values cannot be directly used due to the astrophysical contribution and imperfect knowledge about their distances as well as the Milky Way's gravitational potential [66,67]. The most significant contribution comes from the kinematic Shklovskii effect and the Galactic acceleration contribution. These contributions need to be subtracted using the measurement of the proper motion and the modeling of the Galactic potential [66]. The subtraction introduces extra uncertainties in the intrinsicṖ b parameter. The uncertainty after the subtraction, denoted as σṖint b , is listed in the last column of Table I for each pulsar.
All the pulsars in the table have passed the tests of GR [12]. In particular, the measured orbital decay rates, after subtracting the kinematic Shklovskii effect and the Galactic acceleration contribution, agree with the GR prediction in Eq. (9) within uncertainty for all binary pulsars that we consider. Therefore, here we do not look for evidence of the Galileon radiation, instead we put upper bounds on the graviton mass via constraining the extra Galileon radiation. By plainly assuming that the Galileon radiation is smaller than the uncertainty inṖ b , we obtain upper bounds for the graviton mass for each pulsar. These bounds are plotted in Fig. 6 for cases using σṖint b and σṖobs b , as function of the orbital period P b and the orbital eccentricity e.
The scenario of using σṖint b should be taken as providing the current bounds on m g . In this scenario, the best bound is given by the agreement ofṖ b to the GR prediction within 10 −3 from the recent observation of the double pulsar PSR J0737−3039A [45]. The bound on m g reads, m g 3 × 10 −28 eV/c 2 (95% C.L.) , which is two times better than the bound from the Hulse-Taylor pulsar PSR B1913+16, which gives m g 6 × 10 −28 eV/c 2 (95% C.L.) [59]. Although the larger eccentricity of PSR B1913+16 (e 0.617) is beneficial to the test, compared with a mildly small eccentricity for PSR J0737−3039A (e 0.088), the astrophysical contribution introduces a significant uncertainty to the intrinsicṖ b of PSR B1913+16 [66]. The double pulsar, being relatively close to the Earth, is not yet limited by the astrophysical contribution inṖ b [44].
The scenario of using σṖobs b , on the other hand, should be taken as over-optimistic, representing the ability of clean binary pulsar systems in constraining the graviton mass. This is only possible when the kinematic Shklovskii effect is measured and the Galactic acceleration is modeled both to be better than the observational uncertainty. In such an idealized case, the best pulsar would be PSR B1534+12 [75], which gives m g 6 × 10 −29 eV/c 2 (95% C.L.) by itself alone. These over-optimistic results are also plotted in Fig. 6 with orange circles, but they are only considered to be indicative of technology capability and not used in the following analysis.
Because the theory parameter, m g , is universal to different pulsars, we can combine the constraints from individual pul-  Table I).
sars into a combined constraint. We assume that the measurements for different binary pulsars are independent, thus the covariant matrix is diagonal for them. We make use of a simple (logarithmic) likelihood ln L ≡ − 1 is a sum of Galileon radiations in Eqs. (27)(28)(29), σṖint b is given in Table I, and the summation is over all binary pulsars considered in this paper. We can also assign prior knowledge to the graviton mass within the framework of Bayesian statistics [82]. In Fig. 7 we plot the combined cumulative probability distribution of the graviton mass for two different sets of prior knowledge. The combined bound is dominantly influenced by PSRs J0737−3039A and B1913+16, while the other twelve pulsars only play a minor role. From the figure, we obtain m g 3 × 10 −28 eV/c 2 (95% C.L.) , when a uniform prior on m g ∈ (10 −29 , 10 −27 ) eV/c 2 is assumed, and when a uniform prior on ln m g for m g ∈ (10 −29 , 10 −27 ) eV/c 2 is assumed; the high end of the prior range, namely 10 −27 eV/c 2 , comes from the analysis in Ref. [43]. While the former bound (35) is very robust, the latter bound (36) using a uniform prior on ln m g depends on our choice of prior range. It is a generic feature of Bayesian analysis. In a cosmologically favored reasoning, one might expect the graviton mass to be around the current Hubble scale, namely m g ∼ H 0 ∼ 10 −33 eV/c 2 . If we had used a uniform prior on ln m g for m g ∈ (10 −34 , 10 −27 ) eV/c 2 , we would have obtained a much tighter bound of m g . However, such a bound is dominated by the prior knowledge. It means that binary pulsars are not yet is the measurement uncertainty of the timing parameterṖ b , while σṖint b is the uncertainty ofṖ b after accounting for the kinematic Shklovskii effect and the Galactic acceleration contribution [66,67]. Parenthesized numbers represent the 1-σ uncertainty in the last digit(s) quoted. For more details, interested readers are referred to the cited references alongside the pulsar name for each pulsar.  , respectively. The Hulse-Taylor pulsar PSR B1913+16 [59] and the double pulsar PSR J0737−3039A [45] are annotated. sensitive to a cosmologically small graviton mass. Therefore, we stick to our relatively conservative result in Eq. (36). In contrast, the above change in the prior range does not affect the bound in Eq. (35) when a uniform prior on m g is adopted.
The results on m g here improve the one in Ref. [43], m g

IV. PROJECTED CONSTRAINTS WITH PULSAR-BH SYSTEMS
A well-timed pulsar around a BH companion is a longsought holy grail in the pulsar astronomy. The discovery of these systems will enable a couple of unprecedented tests of gravity theories, in particular on the aspects related to the property of BH solutions [17,46,48,80,81,83]. Up to now, despite extensive dedicated searches, no conclusively convincing candidate has been found. The uncertainty in the estimation of the number of these potential systems pertains mainly to their formation channels. We will not discuss further the involved astrophysics here. Nevertheless, several studies have shown a good potential for the discovery of pulsar-BH systems in the near future [47,49,51,52,84]. If such a pulsar-BH system is discovered, it will provide a completely new playground to perform interesting tests of gravity, including the Galileon radiation [48].
In Fig. 8, we plot the expected contribution from the Galileon radiation to the orbital decay rate for a pulsar -10 M BH system. In this figure, the graviton mass is assumed to be m g = 10 −28 eV/c 2 , which is roughly the best bound from the combination of the current binary pulsars, obtained in Eqs. (35) and (36). Because of the asymmetric masses, for extreme (unrealistic) systems with P b 1 minute and e 0.95, the Galileon dipole radiation could be relevant to the total radiation. However, even if there exist such kind of systems, they are unlikely to be detected with radio telescopes due to the large orbital acceleration and limited computation resources [47]; however, fast imaging and imaging searches based on significant circular polarization or scintillation might help [85][86][87]. The expected number for the LISA detector is very low as well, due to their short lifetime before the merger [64]. For binaries with P b 1 day, the Galileon quadrupole radiation is still the dominant contributor.
Depending on the eccentricity, now the Galileon monopole radiation could play an essential role for binaries with P b 1 day. Discovery of such binaries might be realistic for current and upcoming radio telescopes [47,80]. Liu et al. [47] conducted extensive time-of-arrival (TOA) simulation for a pulsar-BH system with masses (m 1 , m 2 ) = (1.4, 10) M . Similar to Ref. [47], we have conducted extensive mockdata simulations, for different pulsar-BH configurations. For all these simulations we have assumed one observing session per week with ten TOAs of given uncertainty σ TOA , over a period of five years. We further assume that the TOAs follow a Gaussian distribution and are uncorrelated (white noise). Our simulations cover an orbital period range from 0.2 to about 100 days. We find that, with given orbital period, the precision in measuringṖ b , denoted by the (dimensionless) quantity σṖ b , is only weakly dependent on the orbital eccentricity. The dependence on P b is nicely fitted by,  of TOAs are typical for precision timing observations in pulsar astronomy (see e.g. Ref. [70]). The expected precision for the three different values of σ TOA is shown as shaded regions in Fig. 8. We see that pulsar-BH systems will have a great potential to improve the current best bound on m g . For example, if σ TOA about 0.1 µs (10 µs) is achieved, an eccentric pulsar with P b 1 month (P b 1 day) can probe m g down to the level of 10 −28 eV/c 2 . We would like to emphasise that the timing precision assumed here can generally only be achieved for recycled pulsars, even with large radio telescopes like the Five-hundred-meter Aperture Spherical Telescope (FAST) [88,89] or the upcoming Square Kilometre Array (SKA) [80,81,83]. Such pulsar-BH systems could in principle be the result of an exchange encounter in regions with high stellar density, like globular clusters and the Galactic centre region (see e.g. Ref. [90]). The actual timing precision will, in addition, depend on other parameters like pulsar luminosity, pulse shape, etc.
In Fig. 9, we plot the projected bounds on m g using theṖ b precision in Eq. (37) with σ TOA = 0.1 µs. As we can see from the figure, if we can time a near-circular pulsar-BH binary with P b smaller than a few days, the bound (36) in this paper can be improved. If the binary is highly eccentric, then an orbital period P b smaller than a few months is sufficient to improve the bound, as was indicated in Fig. 8. On the other hand, if one wants to probe the cosmologically interested range for m g ∼ 10 −33 eV/c 2 [43], a sub-minute circular binary, or a subhour highly eccentric binary, is needed. We consider such cases highly unlikely to be discovered with near-future technologies, leaving aside the fact that the existence of such a system in our Galaxy is almost certainly excluded, due to its short merger time.
In another direction, extensive searches for pulsars around the Sgr A * , the supermassive BH at the centre of our Galaxy, are ongoing [49,51,52,91]. In Fig. 10 we plot the Galileon radiation for a pulsar -Sgr A * BH system for different eccentricities with m g = 10 −28 eV/c 2 . For eccentric systems with P b 1 yr, the Galileon monopole radiation prevails over the Galileon quadrupole radiation for e 0.6. Therefore, there might be an opportunity to test the Galileon monopole radiation within this class of systems. To address this question, we perform mock-data simulations to investigate the precision one can expect forṖ b for a pulsar in a suitable orbit around Sgr A*. Similar to Ref. [49], we have created mock data with one TOA of 100 µs precision every week, over a time span of five years. Furthermore, we have assumed that the pulsar orbit is unperturbed and therefore our parameter estimation is based on a phase-connected timing solution that provides a perfect fit over the whole time span of observations. Even under such optimistic assumptions, we find that, for an orbit with P b = 0.5 yr and e = 0.8, it is unlikely to get a σṖ b better than 10 −12 s s −1 . Therefore, only in the event of a pulsar in a highly eccentric orbit with P b 1 day being discovered, we might be able to improve the bound (36). However, we consider the existence of a pulsar in such an orbit extremely unlikely.
The tests in this section depend sensitively on the actual eccentricity of the pulsar-BH system, as shown in the Figs. 8-10. In fact, the analysis with Fig. 8 and Fig. 10 is conservative, because the Galileon radiation power was calculated averaging over the orbital timescale [43]. In reality, during the periastron passage the gravitational radiation will be max-imized, and produce prominent features. These features are believed to provide even better distinguishable signals. An analysis resolving the orbital timescale is beyond the scope of this paper. Because all three kinds of the Galileon radiations are proportional to m g , the results discussed in Fig. 8 and 10 can be rescaled easily with different graviton mass for different binary systems.

V. DISCUSSION
In this paper we systematically studied the Galileon radiation in cubic Galileon theory [43] in the context of pulsar timing. Because the Galileon radiation is screened differently than its fifth-force counterpart, such a study is essential to better understand the basic role of the Vainshtein mechanism in screening the gravitational radiation.
From a collection of fourteen well-timed binary pulsars, we have obtained a new bound on the theory parameter for cubic Galileon, namely the graviton mass, m g 2 × 10 −28 eV/c 2 (95% C.L.) , when a uniform prior on ln m g for m g ∈ (10 −29 , 10 −27 ) eV/c 2 is used. This improves a previous bound from the Hulse-Taylor pulsar [43] by a factor of five. Though the bound (38) is weaker than a few other bounds such as bounds from the Earth-Moon-Sun system and dark matter clusters [23], it is nevertheless a robust bound from a completely different regime, namely from the dynamic gravitational radiation instead of the static environments. It is also immune from uncertain assumptions about the dark matter distributions and the virialization of the gravitating systems. Therefore, we consider this bound complementary to the existing bounds. Finally, de Rham et al. [79] have discussed radiations from a generic Galileon theory with all allowed interactions in the 4-dimensional spacetime. The inclusion of quartic or quintic Galileon complicates the calculations considerably. The authors found that, the naive perturbation theory predicts divergent results in the radiation power when these higher-order terms are considered, meaning that the perturbations themselves are nonlinear. Partial results were obtained for binary systems with specific assumptions about the screening lengthscales. In particular, only circular binary orbits were analyzed, thus not applicable to most of the systems that we consider in this paper. For circular orbits that were considered in Ref. [79], meaningful bounds can be derived when there is a hierarchy of strong coupling scales. We wish to perform a more complete analysis with higher-order Galileon interactions in a future study.