Non-local massive gravity

We present a general covariant action for massive gravity merging together a class of"non-polynomial"and super-renormalizable or finite theories of gravity with the non-local theory of gravity recently proposed by Jaccard, Maggiore and Mitsou (Phys. Rev. D 88 (2013) 044033). Our diffeomorphism invariant action gives rise to the equations of motion appearing in non-local massive massive gravity plus quadratic curvature terms. Not only the massive graviton propagator reduces smoothly to the massless one without a vDVZ discontinuity, but also our finite theory of gravity is unitary at tree level around the Minkowski background. We also show that, as long as the graviton mass $m$ is much smaller the today's Hubble parameter $H_0$, a late-time cosmic acceleration can be realized without a dark energy component due to the growth of a scalar degree of freedom. In the presence of the cosmological constant $\Lambda$, the dominance of the non-local mass term leads to a kind of"degravitation"for $\Lambda$ at the late cosmological epoch.


I. INTRODUCTION
The construction of a consistent theory of massive gravity has a long history, starting from the first attempts of Fierz and Pauli [1] in 1939. The Fierz-Pauli theory, which is a simple extension of General Relativity (GR) with a linear graviton mass term, is plagued by a problem of the so-called van Dam-Veltman-Zakharov (vDVZ) discontinuity [2]. This means that the linearized GR is not recovered in the limit that the graviton mass is sent to zero.
The problem of the vDVZ discontinuity can be alleviated in the non-linear version of the Fierz-Pauli theory [3]. The non-linear interactions lead to a well behaved continuous expansion of solutions within the so-called Vainshtein radius. However, the nonlinearities that cure the vDVZ discontinuity problem give rise to the so-called Boulware-Deser (BD) ghost [4] with a vacuum instability.
A massive gravity theory free from the BD ghost was constructed by de Rham, Gabadadze and Tolley (dRGT) [5] as an extension of the Galileon gravity [6]. On the homogenous and isotropic background, however, the selfaccelerating solutions in the dRGT theory exhibit instabilities of scalar and vector perturbations [7]. The analysis based on non-linear cosmological perturbations shows that there is at least one ghost mode (among the five degrees of freedom) in the gravity sector [8]. Moreover it was shown in Ref. [9] that the constraint eliminating the BD ghost gives rise to an acausality problem. These problems can be alleviated by extending the original dRGT theory to include other degrees of freedom [10][11][12] (like quasidilatons) or by breaking the homogeneity [13] or isotropy [14,15] of the cosmological background. * Electronic address: lmodesto@fudan.edu.cn † Electronic address: shinji@rs.kagu.tus.ac.jp Recently, Jaccard et al. [16] constructed a nonlocal theory of massive gravity by using a quadratic action of perturbations expanded around the Minkowski background. This action was originally introduced in Refs. [17,18] in the context of the degravitation idea of the cosmological constant. The resulting covariant nonlinear theory of massive gravity not only frees from the vDVZ discontinuity but respects causality. Moreover, unlike the dRGT theory, it is not required to introduce an external reference metric.
Jaccard et al. [16] showed that, on the Minkowski background, there exists a scalar ghost in addition to the five degrees of freedom of a massive graviton, by decomposing a saturated propagator into spin-2, spin-1, and spin-0 components. For the graviton mass m of the order of the today's Hubble parameter H 0 , the vacuum decay rate induced by the ghost was found to be very tiny even over cosmological time scales. The possibility of the degravitation of a vacuum energy was also suggested by introducing another mass scale µ much smaller than m.
In this paper we propose a general covariant action principle which provides the equations of motion for the non-local massive gravity [16] with quadratic curvature terms. The action turns out to be a bridge between a class of super-renormalizable or finite theories of quantum gravity [19][20][21][22][23][24][25] and a diffeomorphism invariant theory for a massive graviton.
The theory previously studied in Refs. [19][20][21][22][23][24][25] has an aim to provide a completion of the Einstein gravity through the introduction of a non-polynomial or semipolynomial entire function (form factor) without any pole in the action. In contrast, the non-local massive gravity studied in this paper shows a pole in the classical action making it fully non-local. However, the Lagrangian for massive gravity can be selected out from the theories previously proposed [19][20][21][22][23][24][25] once the form factor has a particular infrared behavior. The non-local theory resulting from the covariant Lagrangian is found to be unitary at tree level on the Minkowski background. Moreover, the theory respects causality and smoothly reduces to the massless one without the vDVZ discontinuity.
We will also study the cosmology of non-local massive gravity on the flat Friedmann-Lemaître-Robertson-Walker (FLRW) background in the presence of radiation and non-relativistic matter 1 . Neglecting the contribution of quadratic curvature terms irrelevant to the cosmological dynamics much below the Planck scale, the dynamical equations of motion reduce to those derived in Ref. [16]. We show that, as long as the graviton mass m is much smaller than H 0 , the today's cosmic acceleration can be realized without a dark energy component due to the growth of a scalar degree of freedom.
Our paper is organized as follows. In Sec. II we show a non-local covariant Lagrangian which gives rise to the same equation of motion as that in non-local massive gravity with quadratic curvature terms. We also evaluate the propagator of the theory to study the tree-level unitarity. In Sec. III we study the cosmological implications of non-local massive gravity in detail to provide a minimal explanation to dark energy in terms of the graviton mass. We also discuss the degravitation of the cosmological constant induced by the non-local mass term. Conclusions and discussions are given in Sec. IV.

II. SUPER-RENORMALIZABLE NON-LOCAL GRAVITY
Let us start with the following general class of nonlocal actions in D dimension [19][20][21][22][23][24][25], where κ = √ 32πG (G is gravitational constant), |g| is the determinant of a metric tensor g µν , is the d'Alembertian operator with M = /M 2 , and M is an ultraviolet mass scale. The first two lines of the action 1 Note that cosmological consequences of non-local theory given by the Lagrangian Rf ( −1 R) have been studied in Refs. [26][27][28][29][30][31][32]. In this case the function f ( −1 R) can be chosen only phenomenologically from the demand to realize the late-time cosmic acceleration and so on.
consist of a finite number of operators multiplied by coupling constants subject to renormalization at quantum level. The functions h 2 (z) and h 0 (z), where z ≡ − M , are not renormalized and defined as follows for general parametersã n andb n , while The form factor V (z) −1 in Eq. (3) is made of two parts: (i) a non-local operator ( + m 2 )/ which goes to the identity in the ultraviolet regime, and (ii) an entire function e H(z) without zeros in all complex planes. Here, m is a mass scale associated with the graviton mass that we will discuss later when we calculate the two-point correlation function. H(z) is an entire function of the operator z = − M , and p γ+N+1 (z) is a real polynomial of degree γ + N + 1 which vanishes in z = 0, while N = (D − 4)/2 and γ > D/2 is integer 2 . The exponential factor e H(z) is crucial to make the theory super-renormalizable or finite at quantum level [19][20][21][22][23][24][25].
Let us expand on the behaviour of H(z) for small values of z: For the most simple choice p γ+N+1 (z) = z γ+N+1 , H(z) simplifies to In particular lim z→0 H(z) = 0. We will expand more about the limit of large z in Sec. II B, where we will explicitly show the power counting renormalizability of the theory.

A. Propagator
In this section we calculate the two point function of the gravitational fluctuation around the flat space-time. For this purpose we split the g µν into the flat Minkowski metric η µν and the fluctuation h µν , as Writing the action (1) in the form S = d D x L, the Lagrangian L can be expanded to second order in the graviton fluctuation [33] where Under this shift the graviton field is transformed as h µν → h µν − ξ µ,ν − ξ ν,µ . The presence of this local gauge invariance requires for a gauge-fixing term to be added to the linearized Lagrangian (9). Hence, if we choose the usual harmonic gauge (∂ µ h µν = 0) [21,34] the linearized gauge-fixed Lagrangian reads where the operator O is made of two terms, one coming from the linearized Lagrangian (9) and the other from the gauge-fixing term (11). Inverting the operator O [33], we find the following two-point function in the momentum space (with the wave number k), where we omitted the tensorial indices for O −1 . The operators {P (2) , P (1) , P (0) ,P (0) }, which project out the spin-2, spin-1, and two spin-0 parts of a massive tensor field, are defined by [33] where ω µν = k µ k ν /k 2 and θ µν = η µν − k µ k ν /k 2 . These correspond to a complete set of projection operators for symmetric rank-two tensors. The functions α(k 2 ) and β(k 2 ) are achieved by replacing → −k 2 in the definitions (10). By looking at the last two gauge-invariant terms in Eq. (13), we deem convenient to introduce the following definitions, Through these definitions, the gauge-invariant part of the propagator greatly simplifies to

B. Power counting super-renormalizability
The main properties of the entire function e H(z) useful to show the super-renormalizability of the theory are the following, where we assumed p γ+N+1 (z) = z γ+N+1 . The first limit tells us what is the leading behaviour in the ultraviolet regime, while the second limit confirms that the next to the leading order goes to zero faster then any polynomial. Let us then examine the ultraviolet behavior of the theory at quantum level. According to the property (18), the propagator and the leading n-graviton interaction vertex have the same scaling in the high-energy regime [see Eqs.
(2), (4), (15), (17), and (18)]: vertex : In Eq. (20) the indices for the graviton fluctuation h µν are omitted and h i (− M ) is one of the functions in Eq. (2). From Eqs. (19) and (20), the upper bound to the superficial degree of divergence is In Eq. (21) (1), i.e. κ,λ, a n , b n together with the finite number of couplings that multiply the operators O(R 3 ) in the last line of Eq. (1). We now assume that the theory is renormalized at some scale µ 0 . Therefore, if we set a n = a n ( in Eq. (2), the functions (15) and (16) reduce tō Thus, in the momentum space, only a pole at k 2 = m 2 occurs in the bare propagator and Eq. (17) reads The tensorial structure of Eq. (24) is the same as that of the massless graviton and the only difference appears in an overall factor 1/(k 2 −m 2 ). If we take the limit m → 0, the massive graviton propagator reduces smoothly to the massless one and hence there is no vDVZ discontinuity. Assuming the renormalization group invariant condition (22), missing the O(R 3 ) operators in the action (1), and settingλ to zero, the non-local Lagrangian in a D dimensional space-time greatly simplifies to On using the function α( ) = 2(V ( ) −1 − 1)/(κ 2 ), the Lagrangian (25) can be expressed as If we are interested only in the infrared modifications of gravity, we can fix H(− M ) = 0. This condition restricts our class of theories to the non-local massive gravity.

C. Unitarity
We now present a systematic study of the tree-level unitarity [33]. A general theory is well defined if "tachyons" and "ghosts" are absent, in which case the corresponding propagator has only first poles at k 2 − m 2 = 0 with real masses (no tachyons) and with positive residues (no ghosts). Therefore, to test the tree-level unitarity, we couple the propagator to external conserved stress-energy tensors, T µν , and we examine the amplitude at the poles [35]. When we introduce the most general source, the linearized action (12) is replaced by where g is a coupling constant. The transition amplitude in the momentum space is Since the stress-energy tensor is conserved, only the projectors P (2) and P (0) will give non-zero contributions to the amplitude. In order to make the analysis more explicit, we expand the sources using the following set of independent vectors in the momentum space [33,[35][36][37]: where ǫ i are unit vectors orthogonal to each other and to k. The symmetric stress-energy tensor reads where we introduced the notation The conditions k µ T µν = 0 and k µ k ν T µν = 0 place the following constraints on the coefficients a, b, d, e i , f i [33]: where k 2 := k 2 0 − k 2 . The conditions (31) and (32) imply while the condition (33) leads to Introducing the spin-projectors and the conservation of the stress-energy tensor k µ T µν = 0 in Eq. (28), the amplitude results where T := η µν T µν . The residue at the pole k 2 = m 2 reads Res A If we assume the stress-tensor to satisfy the usual energy condition, then the following inequality follows Using the conditions (35), (36), and (40) in Eq. (39), we find that for D ≥ 3. This shows that the theory is unitary at tree level around the Minkowski background. As we see in Eq. (38) the contribution to the residue from the spin-0 operator P (0) is negative, but the spin-2 operator P (2) provides a dominant contribution with a positive sign of Res A k 2 =m 2 . Hence the presence of the spin-2 mode is crucial to make the theory unitary.

D. Equations of motion
Let us derive the equations of motion up to curvature squared operators O(R 2 ) and total derivative terms [17,[38][39][40][41]. The action of our theory is S = d D xL, where the Lagrangian is given by Eq. (25). The variation of this action reads where we omitted the argument − M of the form factor V −1 . We also used the relations ∇ µ g ρσ = 0, ∇ µ G µν = 0, and The action is manifestly covariant in general. Hence its variational derivative (the left hand side of the modified Einstein equations) exactly satisfies the Bianchi identity Taking into account the energy-momentum tensor T µν , the equation of motion at the quadratic order of curvatures reads Except for the very high-energy regime the quadratic curvature terms should not be important in Eq. (45). Neglecting the O(R 2 µν ) terms and setting e H(− M ) = 1 in Eq. (45), it follows that which is the same equation as that studied in Ref. [16] in the context of non-local massive gravity with the graviton mass m.
If we apply Eq. (46) to cosmology, the d'Alembertian is of the order of ∼ d 2 /dt 2 ∼ ω 2 , where ω is the characteristic frequency of a corresponding physical quantity. Provided ω ≫ m the term m 2 −1 G µν in Eq. (46) is suppressed relative to G µν , so that the Einstein equation G µν ≃ 8πGT µν is recovered. In order to realize the standard radiation and matter eras, it is expected that m should not be larger than H 0 . At the late cosmological epoch, the effect of the non-local term m 2 −1 G µν can be important to modify the dynamics of the system.
If we take the derivative of Eq. (46) by exerting the operator , it follows that This equation is invariant under the symmetry which realizes the Afshordi-Smolin idea [42] for the degravitation of the cosmological constant. Equation (47) does not admit exact de Sitter solutions. There exist de-Sitter solutions characterized by G dS µν = 8πG ρ eff Λ g µν for the modified model in which the operator in Eq. (47) is replaced by + µ 2 , where µ is a small mass scale [16]. If the energy-momentum tensor on the right hand side of Eq. (47) is given by T (Λ) µν = ρ Λ g µν , we obtain the effective cosmological constant ρ eff Λ = ρ Λ µ 2 /(m 2 + µ 2 ). For µ much smaller than m, it follows that ρ eff Λ ≪ ρ Λ . In the limit µ → 0, the effective cosmological constant disappears completely.
The crucial point for the above degravitation of ρ Λ is that both G dS µν and and T (Λ) µν vanish at de Sitter solutions. For the background in which the matter density ρ varies (such as the radiation and matter eras), the two d'Alembertians in Eq. (47) give rise to the contributions of the order of ω 2 . In other words, the above degravitation of ρ Λ should occur at the late cosmological epoch in which ω drops below µ [16].
A detailed analysis given in Sec. III shows that, even for ρ Λ = 0 and µ = 0, a late-time cosmic acceleration occurs on the flat FLRW background. This comes from the peculiar evolution of the term m 2 −1 G µν in Eq. (46), by which the equation of state smaller than −1 can be realized. Even in the presence of the cosmological constant, the non-local term eventually dominates over ρ Λ at the late cosmological epoch. In the following we focus on the theory based on the field equation (46), i.e., µ = 0.

III. COSMOLOGICAL DYNAMICS
We study the cosmological dynamics on the fourdimensional flat FLRW background characterized by the line element ds 2 = −dt 2 + a 2 (t)(dx 2 + dy 2 + dz 2 ), where a(t) is the scale factor with the cosmic time t. Since we ignore the O(R 2 µν ) terms and set H(− M ) = 0 in Eq. (45), our analysis can be valid in the low-energy regime much below the Planck scale.
We introduce a tensor S µν satisfying the relation by which the second term on the left hand side of Eq. (46) can be written as m 2 −1 G µν = m 2 S µν . In order to respect the continuity equation ∇ µ T µν = 0 of matter, we take the transverse part S T µν of the symmetric tensor S µν , that is, ∇ µ S T µν = 0. Then, Eq. (46) can be written as We use the fact that S µν can decomposed as [16,43] where the vector S µ has the time-component S 0 alone in the FLRW background, i.e., S i = 0 (i = 1, 2, 3). From Eq. (51) we have where u ≡ S 0 0 and v ≡ S i i , and a dot represents a derivative with respect to t. In the presence of the matter energy-momentum tensor T µν = (ρ, a 2 P δ ij ), the (00) and (ii) components of Eq. (50) are where H =ȧ/a. Taking the divergence of Eq. (51), it follows that 2∇ µ S µν = ∇ µ (∇ µ S ν + ∇ ν S µ ). From the ν = 0 component of this equation we obtain The (00) and (ii) components of Eq. (49) givë which can be decoupled each other by defining In summary we get the following set of equations from Eqs. (53)-(57): From Eqs. (59)-(61) one can show that the continuity equationρ + 3H(ρ + P ) = 0 holds. For the matter component we take into account radiation (density ρ r , pressure P r = ρ r /3), non-relativistic matter (density ρ m , pressure P m = 0), and the cosmological constant (density ρ Λ , pressure P Λ = −ρ Λ ), i.e., ρ = ρ r + ρ m + ρ Λ and P = ρ r /3 − ρ Λ . Each matter component obeys the continuity equationρ i + 3H(ρ i + P i ) = 0 (i = r, m, Λ).
In order to study the cosmological dynamics of the above system, it is convenient to introduce the following dimensionless variables where r H ≡ H ′ /H, and a prime represents a derivative with respect to N = ln(a/a i ) (a i is the initial scale factor). From Eq. (59) it follows that We define the density parameter of the dark energy component, as Ω DE ≡ Ω Λ + Ω NL . From Eqs. (59) and (60) the density and the pressure of dark energy are given respectively by Then, the dark energy equation of state w DE = P DE /ρ DE can be expressed as From Eq. (60) the quantity r H = H ′ /H obeys by which the effective equation of state of the Universe is known as w eff = −1 − 2r H /3. On using Eqs. (60)-(63) and the continuity equation of each matter component, we obtain the following differential equations In Eq. (71) the derivative of r H is given by Let us first study the case in which the cosmological constant is absent (Ω Λ = 0). We assume that m is smaller than the today's Hubble parameter The field S is amplified by the force term on the right hand side of Eq. (71). Meanwhile the homogeneous solution of Eq. (71) decays for r H = −2 and −3/2. Then, for |V | ≫ 1, the field S grows as which behaves as S ≃ Since r H is constant, the growth rates of the Hubble parameter squared H 2 are the same as those of V and S. Hence we have the super-inflationary solution H ∝ a √ 57/6−1/2 approaching a big-rip singularity. However the above study neglects the contribution of the O(R 2 ) terms, so inclusion of those terms can modify the cosmological dynamics in the high-curvature regime.
In order to confirm the above analytic estimation, we numerically integrate Eqs. (69)-(73) with the initial conditions U = U ′ = 0, V = V ′ = 0, and S = S ′ = 0 in the deep radiation era. In Fig. 1 we plot the evolution of w DE and w eff as well as the density parameters Ω NL , Ω m , Ω r versus the redshift z r = 1/a − 1. Clearly there is the sequence of radiation (Ω r ≃ 1, w eff ≃ 1/3), matter (Ω m ≃ 1, w eff = 0), and dark energy (Ω NL ≃ 1, w eff ≃ −1.5) dominated epochs. We identify the present epoch (z r = 0) to be Ω NL = 0.7. As we estimated analytically, the dark energy equation of state evolves as w DE ≃ −1.791 (radiation era), w DE ≃ −1.725 (matter era), and w DE ≃ −1.506 (accelerated era).
Notice that, even with the initial conditions V = V ′ = 0, the growing-mode solution to Eq. (70) cannot be eliminated due to the presence of the term −2r H . Taking into account the decaying-mode solution to V in the radiation era, the coefficient c 2 of Eq. (76) corresponding to V = V ′ = 0 at N = 0 (i.e., a = a i ) is c 2 = √ 33/132 + 1/4 ≃ 0.29. Up to the radiation-matter equality (a = a eq ), the field evolves as V ≃ 0.29(a/a i ) ( for a eq < a < 1. In the numerical simulation of Fig. 1 the initial condition is chosen to be a i = 2.6×10 −7 with a eq = 3.1×10 −4 . Since the cosmic acceleration starts when the last term in Eq. (68) grows to the order of 1, we have m 2 V 0 /(8H 2 0 ) ≈ 1, where V 0 is the today's value of V . Using the analytic estimation (80), the mass m is constrained to be m ≈ 10 −7 H 0 . In fact, this is close to the numerically derived value m = 1.5 × 10 −7 H 0 .
Thus, the mass m is required to be much smaller than H 0 to avoid the early beginning of cosmic acceleration. If the onset of the radiation era occurs at the redshift z r larger than 10 15 , i.e., a i 10 −15 , the analytic estimation (80) shows that the mass m needs to satisfy the condition m 10 −17 H 0 to realize the successful cosmic expansion history.
If we consider the evolution of the Universe earlier than the radiation era (e.g., inflation), the upper bound of m should be even tighter. On the de Sitter background (Ḣ = 0) we have r H = 0, in which case the growth of V can be avoided for the initial conditions V = V ′ = 0. However, inflation in the early Universe has a small deviation from the exact de Sitter solution [44] and hence the field V can grow at some extent due to the nonvanishing values of r H . For the theoretical consistency we need to include the O(R 2 ) terms in such a high-energy regime, which is beyond the scope of our paper. In the presence of the cosmological constant with the energy density ρ Λ , the cosmological dynamics is subject to change relative to that studied in Sec. III A. During the radiation and matter eras we have 3H 2 ≫ |(m 2 /4)(U + 3V − 4Ṡ 0 )| in Eq. (59) and hence 3H 2 ≃ 8πGρ. In order to avoid the appearance of ρ Λ in these epochs, we require the condition ρ Λ 3H 2 0 /(8πG). The non-local mass term finally dominates over the cosmological constant because the equation of state of the former is smaller than that of the latter. If the condition 8πGρ Λ ≫ |(m 2 /4)(U + 3V − 4Ṡ 0 )| is satisfied today, the non-local term comes out in the future. The case (a) in Fig. 2 corresponds to such an example. In this case, the dark energy equation of state is close to −1 up to z r ∼ −0.9. It then approaches the asymptotic value w DE = −1.506.
For smaller values of ρ Λ , the dominance of the nonlocal term occurs earlier. In the case (b) of Fig. 2 the energy densities of the non-local term and the cosmological constant are the same orders today (Ω NL = 0. 36 and Ω Λ = 0.34 at z r = 0). In this case the dark energy equation of state starts to decrease only recently with the today's value w DE = −1.39.
In the case (c) the transition to the asymptotic regime w DE = −1.506 occurs even earlier (around z r ∼ 100). Observationally it is possible to distinguish between the three different cases of Fig. 2. In the limit that ρ Λ → 0, the evolution of w DE approaches the one shown in Fig. 1. For smaller ρ Λ the graviton mass m tends to be larger because of the earlier dominance of the non-local term. For the cases (a), (b), and (c), the numerical values of the mass are m = 8.37 × 10 −9 H 0 , m = 1.22 × 10 −7 H 0 , and m = 1.46 × 10 −7 H 0 , respectively. In Fig. 3 we plot the evolution of Ω NL and Ω Λ for the initial conditions corresponding to the case (b) in Fig. 2.
After Ω NL gets larger than Ω Λ today, Ω NL approaches 1, while Ω Λ starts to decrease toward 0. This behavior comes from the fact that, after the dominance of Ω NL , the terms on the left hand side of Eq. (59) balance with each other, i.e., 3H 2 + (m 2 /4)(U + 3V − 4Ṡ 0 ) ≃ 0. Then the cosmological constant appearing on the right hand side of Eq. (59) effectively decouples from the dynamics of the system. This is a kind of degravitation, by which the contribution of the matter component present in the energy density ρ becomes negligible relative to that of the non-local term.

IV. CONCLUSIONS AND DISCUSSIONS
In this paper we showed that the field equation of motion in the non-local massive gravity theory proposed by Jaccard et al. [16] follows from the covariant non-local Lagrangian (25) with quadratic curvature terms. This is the generalization of the super-renormalizable massless theory with the ultraviolet modification factor e H(− M ) .
Expanding the Lagrangian (25) up to second order of the perturbations h µν on the Minkowski background, the propagator of the theory can be expressed in terms of four operators which project out the spin-2, spin-1, and two spin-0 parts of a massive tensor field. The propagator (24) smoothly connects to that of the massless theory in the limit m → 0 and hence there is no vDVZ discontinuity. We also found that the theory described by (25) is unitary at tree level, by coupling the propagator to external conserved stress-energy tensors and evaluating the residue of the amplitude at the pole (k 2 = m 2 ).
In the presence of a conserved energy-momentum tensor T µν , the non-local equation of motion following from the Lagrangian (25) is given by Eq. (45). In the lowenergy regime much below the Planck scale the quadratic curvature terms can be negligible relative to other terms, so that the equation of motion reduces to (46) for H(− M ) = 0. We studied the cosmological dynamics based on the non-local equation (46) in detail on the flat FLRW background.
The tensor field S µν , which satisfies the relation (49), can be decomposed into the form (51). In order to respect the continuity equation ∇ µ T µν = 0 for matter, the transverse part of S µν needs to be extracted in the second term on the left hand side of Eq. (46). Among the components of the vector S µ in Eq. (51), the three vector S i (i = 1, 2, 3) vanishes because of the symmetry of the FLRW space-time. In addition to the vector component S 0 , we also have two scalar degrees of freedom U = S 0 0 + S i i and V = S 0 0 − S i i /3. Among these dynamical degrees of freedom, the scalar field V exhibits instabilities for the cosmological background withḢ = 0. Even in the absence of a dark energy component, a late-time accelerated expansion of the Universe can be realized by the growth of V . In order to avoid an early entry to the phase of cosmic acceleration, the graviton mass m is required to be very much smaller than the today's Hubble parameter H 0 . We showed that the equation of state of this "dark" component evolves as w DE = −1.791 (radiation era), w DE = −1.725 (matter era), and w DE = −1.506 (accelerated era), see Fig. 1.
While the above property of the non-local massive gravity is attractive, the evolution of w DE smaller than −1.5 during the matter and accelerated epochs is in tension with the joint data analysis of SNIa, CMB, and BAO [45]. In the presence of the cosmological constant Λ (or other dark energy components such as quintessence), the dark energy equation of state can evolve with the value close to w DE = −1 in the deep matter era (see Fig. 2). In such cases the model can be consistent with the observational data. In the asymptotic future the non-local term dominates over the cosmological constant, which can be regarded as a kind of degravitation of Λ.
Recently, Maggiore [46] studied the modified version of the non-local massive gravity in which the second term on the left hand side of Eq. (46) is replaced by m 2 (g µν −1 R) T , where T denotes the extraction of the transverse part. In this theory the −2Ḣ term on the right hand side of Eq. (63) disappears, in which case the growth of V can be avoided for the initial conditions V =V = 0 (i.e., decoupled from the dynamics). Since the growth of the fields U and S 0 is milder than that of the field V studied in Sec. III A, w DE evolves from the value slightly smaller than −1 during the matter era to the value larger than −1 [46]. It will be of interest to study whether such a theory can be consistently formulated in the framework of the covariant action related to the super-renormalizable massless theory.
While we showed that the theory described by the covariant Lagrangian (25) is tree-level unitary on the Minkowski background, it remains to see what happens on the cosmological background. This requires detailed study for the expansion of the Lagrangian (25) up to sec-ond order in cosmological perturbations about the FLRW background. We leave such analysis for future work.