{'fraction_non_alphanumeric': 0.02955372583434427, 'fraction_numerical': 0.0126474459902434, 'mean_word_length': 4.48005657708628, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Signatures of a Graviton Mass in the Cosmic Microwave Background

There exist consistent low energy effective field theories describing gravity in the Higgs phase that allow the coexistence of massive gravitons and the conventional 1/r potential of gravity. In an effort to constrain the value of the graviton mass in these theories, we study the tensor contribution to the CMB temperature anisotropy and polarization spectra in the presence of a non-vanishing graviton mass. We find that the observation of a B-mode signal consistent with the spectrum predicted by inflationary models would provide the strongest limit yet on the mass of an elementary particle -- a graviton -- at a level of m\lesssim 10^(-30) eV\approx(10 Mpc)^(-1). We also find that a graviton mass in the range between (10 Mpc)^(-1) and (10 kpc)^(-1) leads to interesting modifications of the polarization spectrum. The characteristic signature of a graviton mass in this range would be a plateau in the B-mode spectrum up to angular multipoles of l\sim 100. For even larger values of the graviton mass the tensor contribution to the CMB spectra becomes strongly suppressed.


Introduction
The possibility of a nonzero graviton mass is intriguing and has attracted the attention of theorists for a long time, see e.g. [1,2,3,4,5,6] (see [7] for a recent review). Massive gravitons have a number of peculiar properties, such as the van Dam-Veltman-Zakharov discontinuity [2,3], ghost instabilities [5] and strong coupling effects at unacceptably low energy scales [6] which complicate the construction of a sensible theory. However, in light of the cosmological constant problem and coincidence problems between baryonic matter, dark matter, and dark energy, there has been increased interest in massive gravity theories, and it was found that models of modified gravity with Lorentz-violating graviton mass terms may avoid all of these problems [8,9,10]. The effective field theories remain valid up to reasonably large energy scales provided a large enough subgroup of the full diffeomorphism group is left unbroken by the graviton mass.
From a phenomenological point of view, the class of models characterized by the residual local symmetry x i → x i + ξ i (t) is of particular interest. Among the possible choices of residual subgroups required by consistency, this is the only choice giving rise to massive gravitational waves [10]. 1 Some phenomenological and cosmological consequences of these models were studied in Refs. [11,12,13,14]. 2 In particular, late time cosmological attractors have been found where an additional dilatation symmetry t → λt, x i → λ −γ x i gets restored, where γ is a real constant. The main properties of the system in the vicinity of these attractors can be summarized as follows.
(i) The evolution of the background cosmology is described by the usual Friedmann equation.
For suitable choices of the Lagrangian in the symmetry breaking sector that leads to the graviton mass there will be an additional "dark energy" component that may contribute to the observed accelerated expansion of the universe.
(ii) The equations describing the evolution of the scalar and vector fluctuations coincide with those of General Relativity.
(iii) The equation describing the evolution of the tensor fluctuations is modified by the presence of a mass term. In other words, the dispersion relation for gravitational waves takes the form ω 2 = p 2 + m 2 g .
As a consequence of (i) and (ii), the strongest bounds on the graviton mass, m g , in this class of models come from direct or indirect observations of gravitational waves. Since these observations are very limited so far, this opens up the possibility for m g to be quite large. The best constraint 1 Note that in the Lorentz-violating theories a presence of graviton mass terms does not yet imply that the tensor modes are massive. For instance, the graviton mass term m 2 00 h 2 00 just fixes the gauge h 00 = 0. 2 A related class of bigravity models was considered in [15,16]. on the graviton mass in these models currently comes from indirect evidence for the emission of gravitational waves from binary pulsars timing [17] yielding an upper limit of m g 10 4 pc −1 ≈ 3 × 10 −15 cm −1 ≈ 6 × 10 −20 eV. (1) Relatively large graviton masses m g 0.1 pc −1 could be detected by observing the characteristic signature of a large graviton mass, a strong monochromatic signal in gravitational wave detectors due to relic gravitons at a frequency equal to the graviton mass [11]. This signal might be observed either by LISA or using millisecond pulsar timing data. (See, e.g. [18] for a recent analysis.) Graviton masses close to the bound (1) can also be found by using higher frequency gravitational wave detectors to measure a time delay between optical and gravitational wave signals from a distant source.
Finally, these theories may give rise to non-universality of high multipoles of the galactic black hole metric. If present, these may be detected by LISA [14].
The characteristic energy scale Λ of the symmetry breaking sector that leads to a massive graviton is of order Λ ∼ M P l m g . This scale is of the same order as the energy density scale of the Universe now if m g ∼ H 0 where H 0 ≈ 0.7 × (3 Gpc) −1 is the current Hubble constant. 3 In other words, if the observed acceleration of the Universe were due to some modification of gravity that gives rise to a graviton mass, one would expect the mass to be of order H 0 , far below the reach of the experiments mentioned above. This motivates us to look for signatures that are sensitive to much smaller values of the graviton mass.
One expects that a small graviton mass will leave an imprint in the temperature anisotropy and polarization spectra of the cosmic microwave background (CMB). We study the contribution of tensor perturbations to the CMB spectra in a modified gravity theory with the properties outlined above and show that this is indeed the case. We treat the graviton mass as a phenomenological parameter. Our results will therefore be valid for any theory in which a modification of gravity amounts to massive gravitational waves with the ΛCDM cosmological background being unchanged.
The paper is organized as follows. After a summary of the basic equations, we start with an analytic discussion of the CMB spectrum in a massive gravity theory in subsection 2.1. After an analytic discussion of the general properties of the B-mode spectrum, we analytically calculate the contribution to the CMB B-type polarization for low multipole coefficients ignoring the effects of reionization and show that the characteristic feature in the low ℓ range is a plateau. For a range of masses above the Hubble rate at the time of recombination this contribution even dominates over the contribution from reionization and is a good approximation to the full numerical spectrum for low multipole coefficients. We then calculate the contributions to the CMB temperature anisotropy and polarization from a spatially homogeneous tensor mode. The existence of such a contribution is a very unusual property of massive gravity theories. In the massless case, tensor modes are frozen as long as their wavelength is larger than the Hubble scale H −1 where H(t) ≡ a −1 (da/dt) and a(t) is the scale factor of an isotropic Friedmann-Robertson-Walker (FRW) cosmological model. Since temperature anisotropies get generated only by a time varying tensor mode, this implies that very long wavelengths cannot contribute. 4 If the graviton is massive, however, when the expansion rate of the universe drops below the graviton mass, these modes acquire an oscillatory time dependence with a frequency set by the graviton mass, and a decreasing amplitude due to Hubble friction. As a consequence, long wavelength modes in a massive gravity theory will generate a temperature anisotropy quadrupole that will get converted into a polarization quadrupole during recombination if the mass is in the right range and more efficiently once the universe becomes reionized. The zero mode does not contribute to the B-mode spectrum, but does contribute to the temperature T , polarization E-mode and T E cross-correlation quadrupoles.
We then proceed to a numerical treatment in subsection 2.2. The most interesting result is that a graviton mass would strongly modify the shape of the B-mode spectrum for ℓ < 100. If such modifications are observed in experiments such as CMBPol, it would provide strong support for massive gravity theories. If on the other hand the observed spectrum is consistent with General Relativity, this would imply an upper bound on the graviton mass of m g < (10Mpc) −1 ≈ 10 −30 eV . We conclude with a brief summary of the signatures of a graviton mass in the CMB in Section 3.

Tensor Contribution to the CMB for a Massive Graviton
Since its discovery about 45 years ago, the cosmic microwave background radiation has greatly improved our understanding of the very early universe. With the help of current and future experiments this trend is likely to continue. Especially the detection of a B-mode signal would provide valuable information that could be used to verify and explore the details of inflation. The B-mode signal would also provide us with a powerful tool to detect or constrain cosmic strings, and, as we shall see, would also help us put tight bounds on massive gravity theories. In the near future the best constraints on or possible detections of a B-mode signal are expected to come from experiments such as Spider, PolarBeaR, EBEX, BICEP/SPUD, C ℓ OVER, from the Planck satellite, and in the slightly more distant future hopefully from a mission like CMBPol [19].
The properties of the early universe are encoded in correlations of the temperature anisotropies and polarization patterns at different points in the sky. The quantities most commonly used to represent the two-point correlations are the TT as well as the TE, EE, and BB multipole coefficients. The contribution of the tensor fluctuations to them is given by These formulas are equivalent to those of Zaldarriaga and Seljak [20] up to integration by parts. 5 While this makes no difference in the massless case, in the massive case the formulas as written here are better suited for numerical calculations.
In these equations, q = pa(t) is the comoving momentum; τ = dt/a(t) is the conformal time; T 0 = 2.725 K is the microwave background temperature at the present conformal time is the probability distribution of last scattering (or visibility function), withκ(τ ) the photon collision frequency (or differential optical depth); τ 1 is any time taken early enough before recombination so that any photon present at τ 1 would have collided many times before the present; and Ψ(q, τ ) is the "source function," which is customarily calculated 5 We also use slightly different conventions. Their gravitational wave amplitude h and power spectral function P h (k) are related to our gravitational wave amplitude D q (τ ) by h √ P h = D/2. In consequence, their function Ψ √ P h is 1/4 times our source function Ψ.
from a hierarchy of equations for partial-wave amplitudes [21,22] (the dot means the derivative with respect to τ ):∆ with Alternatively, the source function can be calculated from an integral equation [23,24], and we have used both approaches.
Here D q (τ ) is the gravitational wave amplitude (apart from terms that decay outside the horizon), defined by with β(q, ±2) and e ij (q, ±2) the stochastic parameter and polarization tensor for helicity ±2, normalized so that and forq in the 3-direction Finally, d(q, τ ) is the quantity In massive gravity, the evolution of the gravitational wave amplitude is described by the solution to the equation for a minimally coupled massive scalar field [11,12] In conventional cosmological perturbation theory, the solution of Eq. (13) which remains finite for a → 0 and has a wavelength larger than the Hubble scale H −1 is not observable locally and does not contribute to the CMB anisotropy and polarization spectra. On the other hand, the other, decaying graviton mode produces locally measurable effects even in the k → 0 limit. Since it does decay, it is usually assumed to be negligible by the time of recombination.
In the massive gravity case, Eq. (13) implies that a homogeneous metric perturbation starts to oscillate with a frequency equal to m g when the expansion rate H drops below the graviton mass m g . This is completely analogous to what happens for light scalar fields (e.g. axion, moduli, . . . ). An important difference, however, is that in the case of the graviton these oscillations may directly affect the CMB spectra (or, if the mass m g is high enough, may be observed by gravitational wave detectors), similar to a super-Hubble decaying mode or a generic sub-Hubble tensor perturbation in the massless limit. Indeed, the presence of a zero mode implies that, superimposed upon a conventional Hubble expansion, spatial metric components experience anisotropic (but homogeneous) high frequency oscillations with a small amplitude. The effect of such oscillations on the CMB spectra can easily be understood analytically and we will return to this at the end of subsection 2.1.

Analytic results for low multipole coefficients
Before presenting the results of the numerical calculations in the next subsection, let us start with a brief analytic discussion of the spectrum. For the most part, we will limit ourselves to the contribution to the spectrum generated during recombination and ignore the effects of reionization. In the massless case, the effects of reionization give the dominant contribution to the spectrum for ℓ 20 but leave the higher multipole coefficients unchanged (or rather change them trivially by an overall rescaling by e −2τ reion , where τ reion = 0.087 ± 0.017 is the optical depth of the medium due to reionization and despite the clash of notation should not be confused with conformal time), so that ignoring the effects of reionization is good as long as one is interested in ℓ > 20. In the massive case, our numerical results indicate that for a range of masses above the Hubble rate at recombination the contribution from recombination provides a good approximation to the spectrum even at low ℓ providing additional motivation for this simplifying assumption.
We will focus on the B-mode spectrum because it is the most interesting one from an experimental point of view. The discussion could be straightforwardly extended to include the TT, TE, and EE spectra, but we limit ourselves to the contribution of the zero mode to those.
Depending on their comoving momentum, the modes fall into one of two classes or one of three classes depending on whether the mass is smaller than the Hubble rate at recombination or larger than that. Class II corresponds to modes that are non-relativistic as they enter the horizon and during their subsequent evolution. Class III corresponds to modes that enter the horizon when they are relativistic but become non-relativistic before recombination.
For masses below the Hubble rate at recombination, the first possibility is that modes are relativistic at the time they enter the horizon. In this case they will still be relativistic during recombination. These modes are essentially unaffected by the graviton mass, and the spectrum for the values of ℓ these modes contribute to is expected to agree with the one in the massless case. The second possibility is that the modes enter the horizon when they are already non-relativistic. The multipole coefficients these modes contribute to will be different from the ones for the massless case and we will discuss those in more detail below.
For masses above the Hubble rate at recombination there is a third option. The modes can be relativistic as they enter the horizon but become non-relativistic by the time of recombination.
The results are summarized in terms of the range of multipole coefficients the different classes affect for a given mass in Figure 1.
Let us now discuss these regimes in more detail. Both for masses below and above the expansion rate during recombination, the short-wavelength modes that are in the relativistic regime during recombination, are not affected by the graviton mass and will lead to the same spectrum as in the massless case.
In terms of the multipole number ℓ ≈ q(τ 0 −τ rec ) this transition to the massless regime corresponds to where z r ≈ 1088 is the redshift at recombination, and we have used the five-year WMAP values for the cosmological parameters [27]. In particular, equation (15) implies that, ignoring the contribution generated during reionization, the B-mode spectrum is not modified for masses smaller than ∼ 300H 0 , which is the scale corresponding to the size of the visible patch of the Universe during recombination. For larger masses, but still smaller than the expansion rate at recombination, m g < H(τ r ) ≈ 2 × 10 4 H 0 , the modes that are affected by the graviton mass are superhorizon during recombination because they satisfy q a(τr) m g < H(τ r ). As a consequence these modes do not oscillate during recombination. Nevertheless, just like in the massless case, the corresponding source termḊ q in (6) is non-zero (though small) and some amount of polarization is being generated. In contrast to the massless case, however, where the value ofḊ q at recombination is determined by the value of q alone, and goes to zero as q → 0, in the massive caseḊ q depends on (q 2 + m 2 g a(τ r ) 2 ) and is independent of q for long wavelength modes leading to an enhancement of the spectrum for ℓ < ℓ 0 .
For values of the graviton mass larger than the Hubble rate at recombination, m g H(τ r ), all modes start to oscillate before recombination. The modes with long wavelengths start to oscillate as soon as the expansion rate of the universe drops below the graviton mass, i.e. at a time τ m < τ r , such that H(τ m ) = m g .
In particular, for (q/a) 2 ≪ m g H at the matter dominated stage, where t ∝ τ 3 . Shorter modes start to oscillate when they enter the horizon just like in the massless case. The transition between these two regimes happens at q m = m g a(τ m ).
To a good approximation, all modes with momenta smaller than q m have the same evolutionthey are frozen until τ m , and oscillate afterwards with a frequency set by the mass. This value of comoving momentum, q m , corresponding to the transition between class II and III translates to a value in ℓ-space of Here z m is the redshift corresponding to τ m and is determined by condition (16) which can be written more explicitly as At m g = H(τ r ) the multipole number ℓ m coincides with ℓ 0 and takes a value of around ℓ m ∼ 65. At higher masses these two scale are different. According to (15) ℓ 0 grows linearly with mass. On the other hand, ℓ m grows more slowly, ℓ m ∝ m 1/3 g for masses that become relevant during matter domination (i.e., for m g 1.5 × 10 5 H 0 ), and ℓ m ∝ m 1/2 g at higher m g . For masses much larger than this, all modes oscillate rapidly during recombination. As a consequence the polarization signal gets averaged out and becomes strongly suppressed.
Modes corresponding to angular scales between ℓ m and ℓ 0 enter the horizon when they are still relativistic, but become non-relativistic before recombination. As a result, they are still expected to exhibit the conventional oscillation pattern in the angular spectrum, but the phase of oscillations is different because the oscillations at late times are driven by the mass rather than the spatial momentum.
After this discussion of various regimes, let us take a more detailed look at the spectrum. As discussed, for the modes referred to as class I in Figure 1 that are relativistic during recombination the spectrum to a good approximation agrees with the one in the massless case. While a number of analytic results for the temperature anisotropy and polarization have been found for this case, we will not review those here and refer the interested reader to the literature [28,29,30,31,32].
As can be seen by inspection of equation (13), for the modes referred to as class II in Figure 1 the dependence of the gravitational wave amplitude on comoving momentum is trivially given by that of the power spectrum because they are frozen as long as they are outside the horizon and the mass already dominates by the time they enter. This does not in general guarantee that the same is true for the source function as it will generically develop its own q-dependence. To a good approximation the momentum dependence of the source function during recombination is the same as that of the gravitational wave amplitude provided the comoving momentum of the mode is less than the duration of recombination in conformal time, i.e. q∆τ rec ≪ 1. This is satisfied for all modes in class II for the range of masses we are interested in and hence does not provide an additional constraint.
For modes in class II, the dependence of the gravitational wave amplitude and that of the source function on comoving momentum are then trivially given by that of the power spectrum, implying e.g. for a standard inflationary scenario that Ψ(q, τ )q 3 2 − n T 2 is q-independent. This allows us to evaluate the expression for C T BB,ℓ given by equation (2) analytically. Conventionally, one first evaluates the integrals over conformal time and then integrates over momentum. For us it will be more convenient to perform the integral over momentum first. This is possible by rewriting the square of the integral over time as an integral in a plane and using the identity: where J ν (ρ) is the Bessel function of the first kind. The resulting four integrals over q can then be done exactly using an integral known as the Weber-Schafheitlin integral (see e.g. [33]). Dropping terms of order and higher in the terms in the integrals over conformal time arising from the integral over comoving momentum, and setting n T = 0 for simplicity, one finds the following expression for the power spectrum: The ℓ-dependence is now explicit and it is easy to see that for ℓ 10 this becomes independent of ℓ. The ℓ-independent quantity I is defined as 7 and it encodes the dependence of the spectrum on the mass through the dependence of the source function on the mass. The quantity I 2 as a function of mass is shown in Figure 2 for a scalar amplitude of ∆ 2 R = 2.41 × 10 −9 , a tensor-to-scalar ratio r = 1. The oscillatory features seen in the plot can be understood from the fact that the tensor perturbations of the metric take the form given in equation (17). In particular, they are regular in the limit t → 0 and the "decaying" mode, which diverges in this limit, is absent. In turn, this property (which also takes place for larger values of q, (q/a) 2 ≥ m g H) is a consequence of local isotropy of the Universe at very early times. We assume the latter to be produced by inflation and use the inflationary prediction for the primordial power spectrum, but the existence of the oscillations in I, as well as the existence of oscillations in the multipole power spectra of polarization and temperature anisotropy seen in Figures 3 and 4 below ("primordial peaks", similar to the well known acoustics peaks produced by scalar perturbations but with approximately twice less asymptotic period in ℓ, T ℓ = π(τ 0 − τ r )/τ r ≈ 140 [34]), is a more general phenomenon not depending on how this early time isotropy was achieved. As we will derive shortly, what enters into the source function in a crucial way is the tensor perturbation of the metric evaluated at the time of recombination. This quantity viewed as a function of the graviton mass oscillates around zero, implying that the integral does, too. After squaring, this will give rise just to what is seen in Figure 2.
Ignoring terms higher order in the quantity (21) is typically only a good approximation for ℓ 30, but our numerical calculations show that the plateau persists to higher values of ℓ. For the values of masses where I is close to zero, the higher order terms in this expansion give the leading contribution and have to be included even for low ℓ. Furthermore, for these ranges of masses the effects of reionization become important. Especially for masses below the expansion rate at recombination, reionization is important as it will increase the sensitivity of polarization measurements for masses near the low end of the accessible mass range m g ∼ 300H 0 by about one order of magnitude to m g ∼ 30H 0 .
After discussing the B-mode signal in a massive gravity theory, let us briefly discuss another unusual feature in these theories, a contribution of modes of extremely long wavelength to the TT, EE, and TE quadrupoles. In the limit of vanishing momentum, the Boltzmann hierarchy, i.e. equations (6), (7) become very simple. Only∆

(T )
T,0 (q, τ ) and∆ (T ) P,0 (q, τ ) get generated while all others remain zero. Similar equations were considered in Ref. [35], where CMB polarization and anisotropy in the Kasner universe were studied. To calculate the contribution to the CMB temperature anisotropy and polarization in this limit, it is convenient to write equations (6), (7) as where we have defined the integral optical depth κ(τ ) as From equations (3),(4), and (5) we see that the contribution of the zero mode to the TT, EE, and TE quadrupoles is thus given by It is straightforward to solve equations (6) and (7) for∆ where I 1 and I 2 are the following integrals It is convenient to discuss the behavior of the integrals (28), (29) first neglecting the contribution from the reionization epoch. In the absence of reionization, the functions e −κ(τ ) , e − 3 Let a(τ r )∆τ ∼ 35 kpc be the characteristic width of these step functions (or the duration of the recombination epoch). The relevant parameter which determines the behavior of integrals (28), (29) is then δ = m g a(τ r )∆τ .
For small masses, such that δ ≪ 1 (but, of course, assuming m g > H 0 ) one has Consequently, in this case polarization is negligible, while∆ T,0 (τ 0 ) = 2D(τ r ) gives the temperature anisotropy quadrupole. Note that for a fixed initial amplitude of the metric perturbation, the anisotropy is smaller for larger m g . For large masses, δ ≫ 1, the mode oscillates rapidly even during recombination. As a result, both integrals I 1 and I 2 are very small (∝ e −δ ) and both temperature and polarization quadrupoles are negligible.
The largest amount of polarization is generated for δ ∼ 1. In this case recombination cannot be treated as instantaneous, so that there is no cancellation between the two terms in the expression for the polarization∆ (T ) P,0 (τ 0 ). On the other hand, metric oscillations during recombination do not wash out the whole effect yet, and one gets comparable contributions to polarization and temperature anisotropy.
Finally, let us include the effect of reionization. In the presence of reionization Eq. (30) does not hold even for small masses. Instead, one has 8 where τ reion is again the optical depth of the medium due to reionization. As a result, both integrals get somewhat suppressed. On the other hand, the two terms in Eq. (27) no longer cancel, and the contributions to polarization and temperature anisotropy can be of the same order.

Numerical results
To calculate the angular power spectra for values of ℓ > 50, we use CMBfast [20,36] with the evolution equation for the tensor perturbations modified according to Eq. (13). For the low multipole coefficients, ℓ ≤ 50, we use the CMBfast source function, but perform the line of sight integration in Mathematica using equations (2), (3), (4), and (5) because the CMBfast results become unreliable for ℓ ≤ 50 at least for large graviton masses. BB,ℓ multipole coefficients for the range of masses that lead to the most interesting signal in the CMB. The masses are given by m g = µ × 3000H 0 , where µ is given in the legend. Longer dashes correspond to larger mass. All plots are for a scalar amplitude ∆ 2 R = 2.41 × 10 −9 , a tensor-to-scalar ratio, r = 1, and a tensor spectral index n T = 0. For the remaining cosmological parameters parameterizing the background, we use the five-year WMAP values [27].
The issue arises because of rapid oscillations of the source function for all values of comoving momentum. After the integration by parts as implemented in CMBfast, eq. (29) of [37] involves second derivatives of the source function, which are unpleasant to deal with numerically. As a result, the line of sight integration, as implemented in CMBfast, produces unreliable results. The problem exists for the BB spectrum as well but is especially severe for the low-ℓ parts of the EE and TE spectra, where it appears at masses of order m g ∼ 3000H 0 .
Using independent Mathematica code, we also checked that the source function as produced by the modified CMBfast is reliable. We assume a scale-invariant power spectrum for the tensor perturbations, n T = 0. The results for a range of masses are shown in Figure 3.
In all the plots we drop the quadrupole, because its value depends on the IR cutoff at low momenta, as follows from the discussion in section 2.1 (see also section 3 for more details). We have used ∆ 2 R = 2.41 × 10 −9 and have set the tensor-to-scalar ratio, r, to unity. In agreement with our estimates, the effect of the mass is rather mild for masses much below the Hubble rate during recombination and is present only for very low ℓ.
For masses approaching the Hubble rate during recombination, the spectrum is significantly modified up to ℓ ∼ 100, and at m g = 3 × 10 4 H 0 the characteristic plateau at ℓ 100 is fully developed. As we increase the mass further, the height of the plateau increases up to values of µ ≡ mg 3000H 0 ≈ 25 but starts to decrease beyond that in agreement with the oscillations we saw in our semi-analytic result in Figure 2. The origin of the oscillations is that depending on the mass the metric perturbation enters recombination in different phase.
In agreement with our qualitative arguments summarized in Figure 1, for masses µ 10 we see that a transition region appears between the multipole moment ℓ m where the plateau ends and the multipole moment ℓ 0 where the massless spectrum is approached.
We are not showing the spectra at higher masses, because the polarization signal becomes strongly suppressed for µ ≫ 150 because of rapid oscillations of the metric during recombination in agreement with Figure 2.
So far we have only shown the C T BB,ℓ multipole coefficients as they are the most interesting from a phenomenological point of view. In Figure 4, we show a comparison of all four CMB spectra (temperature anisotropy, E-and B-type polarization and T E cross-correlations) for µ = 10 and the massless case.
We see that the mass affects both E-and B-type polarization in a rather similar way. The effect of the mass on the temperature anisotropy is rather mild and the shape of the spectrum for massive gravity is very similar to the massless case. Unlike polarization, the temperature anisotropy receives contributions not only from recombination and reionization, but from all times. As a result, the contribution at horizon crossing dominates and one obtains a plateau reflecting the flatness of the primordial spectrum in both the massive and the massless case (see [28] for the analytic expression describing this plateau in the massless case which, as is seen from the upper left plot in Figure 4, produces a good approximation to the massive case, too). While the polarization spectra look rather different from the ones in the massless case, one should have in mind that we only show the tensor contribution to the signal. The main component for the temperature anisotropy as well as the TE and EE spectra comes from the scalar perturbations which are identical to the ones in general relativity. One may still wonder how large a tensor signal one could tolerate in the massive case given the existing data and what has already been ruled out. To this end, we perform a Markov chain Monte Carlo study for a single value of mass m g = 3 × 10 4 H 0 , or equivalently µ = 10. We use the publicly available CosmoMC code [38,39] to sample the parameter space together with CAMB [40,41] to generate the spectra for a given set of cosmological parameters, and we use a modified version of the WMAP likelihood code that is now available on the LAMBDA website [42] to evaluate the likelihood function for a given spectrum.
In addition to varying the six parameters of the ΛCDM model and marginalizing over the Sunyaev-Zeldovich amplitude, we allow the tensor-to-scalar ratio to vary but keep the mass and the tensor spectral index fixed. We do not implement the slow-roll consistency condition but set n T = 0 as one may expect the consistency condition to be modified in these theories, but this does not significantly change the results.
We find that the tensor-to-scalar ratio for this value of mass is constrained to r < 0.11 at 95% confidence level by the five-year WMAP data alone. The results are shown in Figure 5.
Different from the massless case [27], where adding additional data sets like BBN, supernovae, or baryon acoustic oscillations significantly lowers the allowed tensor-to-scalar ratio, the results do not change much for the massive case. In the massless case, the reason for the significant  Figure 5: These plots show marginalized likelihood plots obtained from a Markov chain Monte Carlo study of a massive gravity model with a mass m g = 3 × 10 4 H 0 , or equivalently µ = 10 using the five-year WMAP data. The dark and light blue contours correspond to 68% and 95% confidence level, respectively.
strengthening of the bound is that additional data sets constraining the baryon or dark matter abundance break a chain of degeneracies. The main constraint on the tensor-to-scalar signal in the massless case currently comes from the low-ℓ TT spectrum, where the tensor signal has a plateau as can be seen e.g in Figure 4. Raising n S lowers the power in the scalar spectrum at low ℓ thus allowing for larger r. On the other hand raising n S becomes possible only because of a degeneracy between Ω B and n S so that adding the BBN priors on the baryon abundance eliminates this possibility and substantially lowers the bound on the tensor-to-scalar ratio in the massless case. In the massive case, however, the polarization data is starting to constrain the model, eliminating the degeneracy in a different way. Including a BBN prior then does not significantly lower the bound in this case and all in all the bound is roughly a factor of two stronger for m g = 3 × 10 4 H 0 than it is in the massless case. As one might expect from looking at the spectra in the massive case, there is now a degeneracy between the tensor-to-scalar ratio and the optical depth, however. This is shown in Figure 5. We interpret these results as telling us that the model is not in conflict with present data but not much more. In a more serious analysis, the mass should certainly not be taken as fixed but be thought of as an unknown parameter that has to be extracted from the data. We leave a more systematic study for when more data becomes available.

Discussion
To summarize, we see that a detection of the CMB B-mode signal either with Planck or with next generation CMB measurements such as CMBPol [19], in addition of opening a new observational window on inflation, will also provide a sensitive probe of the graviton mass. We showed that the most interesting consequence of the graviton mass for the cosmic microwave background is probably the characteristic plateau in the B-mode spectrum for multipoles with ℓ 100. This plateau is most pronounced for masses a few times the Hubble rate at recombination, but in principle CMB polarization measurements are capable of constraining the graviton mass down to m −1 g ∼ 10 Mpc. Taking into account that large graviton masses m −1 g ≪ 10 kpc lead to a strong suppression of the tensor contribution to the CMB spectra, we conclude that the observation of B-mode with the conventional inflationary spectrum would provide by far the tightest bound on the mass of an elementary particle -a graviton -at a level of m g 10 −30 eV.
Of course, even more exciting would be to find out that a graviton mass is actually nonzero. It is worth stressing that gravitational waves after being produced during inflation remain practically undisturbed throughout the later evolution of the Universe [43] and secondary sources of the B mode, such as the weak lensing contribution, are negligibly small at ℓ < 100 [44]. It also appears extremely hard to mock up the effect of the mass by modifying an inflationary model or invoking another mechanism generating gravitational waves such as cosmic strings. Consequently, the detection of a B-mode signal with the shape discussed above would provide an unambiguous signal of a graviton mass.
We have also seen that in a massive gravity theory superhorizon tensor perturbations are physical and contribute to the quadrupole of TT, EE, and TE spectra. As a consequence, the amplitude of these quadrupoles are IR sensitive and can be significantly enhanced over the rest of the spectrum. The ratio between the quadrupole and the rest of the spectrum is model dependent. It provides a probe of the total duration of the inflation provided the tensor perturbations are naturally set to zero at early times in a given model. This is the case, for instance, if the graviton mass is not constant during inflation and in the beginning is bigger than the expansion rate. One way this possibility may be realized is if the dilatation symmetry t → λt, x i → λ −γ x i is not present in the full theory, but only gets restored during inflation as one approaches a cosmological attractor discussed in [12]. In this particular scenario, the quadrupoles measure the number of e-folds during the period for which the graviton is lighter than the Hubble rate.

Acknowledgments
It is a pleasure to thank Eiichiro Komatsu, Peter Tinyakov and Matias Zaldarriaga for useful discussions. SD thanks the Theory Group at the University of Texas for hospitality during the early stages of this project. RF would like to thank the Stanford Institute for Theoretical Physics for hospitality during the early stages of this work and the Aspen Center for Physics during the late stages. The work of SD and RF has been partially supported by the National Science Foundation under Grant No. PHY-0455649. AS was partially supported by the grant RFBR 08-02-00923 and by the Scientific Programme "Elementary Particles" of the Russian Academy of Sciences.