Covariant constraints of massive gravity in metric formulation

We propose a simple method for deriving the constraints of the de Rham-Gabadadze-Tolley model in the metric and the Lagrangian formulation, as possible as keeping the Lorentz covariance. In our formulation, it is not necessary to use the Hamilton analysis, the vielbein formulation, nor the Stuckelberg trick for showing the Boulwer-Deser ghost-freeness. It realizes the Lorentz covariant expressions of the constraints in a certain parameter region.


I. INTRODUCTION
In a last decade, the understanding of massive spin-two field has been quite developed. The unique linear theory describing the massive spin-two particle in the flat spacetime has been well-known since a long time ago. It is called the Fierz-Pauli (FP) model [1]. The equations of motion (EoM) of the FP model can be expressed as follows, Although the second equations can be rewritten as ∂ µ h µν = 0 by using the third equation h = 0, we adopt this expression in (1) for a later convenience. In the Fourier space, the first equation determines a dispersion relation, and the other equations can be regarded as the constraints which reduce the components of the polarizations. Hence, in the D dimensional spacetime, the symmetric tensor h µν obeying the FP model has the (D + 1)(D − 2)/2 components. In four dimensions, the FP model has five degrees of freedom (DoF), and it coincides with the DoF of the spin-two representation of the little group SO (3). Indeed, in four dimensions, the wave equations (1) can be derived from the representation theory of the Poincare group ISO(1, 3). Hence, (1) can be regarded as a definition of the wave equation for the symmetric second rank tensor describing the massive spin-two particle. In 1970's, from a negative standpoint, it began to ask whether or not the gravity can be described by the massive spin two particle. It had been pointed out that the original FP model cannot explain the observational results in the solar system [2]. This is because the massless limit of some observables in the FP model does not coincide with those of the massless model. This discontinuity is called the vDVZ discontinuity. Although it seems that the possibility of the non-vanishing graviton mass had been excluded from the above fact, Vainshtein proposed that a class of nonlinearlization of the FP model can avoid the vDVZ discontinuity [3]. In order to clarify the latter explanation, we would like to define the "massive gravity" generally as the sum of the Einstein-Hilbert action with dynamical metric g µν and non-derivative potential terms which are general functions of the dynamical metric g µν and the flat metric η µν , i.e., Here, we assume a general nonderivative potential term √ −gV (g, η) which is expandable with respect to h µν = g µν − η µν . In the expansion, we assume also that the quadratic order terms with respect to h µν are proportional to the Fierz-Pauli mass term h 2 − h µν h µν and the first order terms vanish. Then, the action (2) can be regarded as a nonlinear extension of the FP model in flat spacetime. We should note that there are no general covariance in the action (2). But we assume that the V (g, η) which has the Lorentz covariance.
In [3], Vainshtein considered the model (2) where the potential term √ −gV (g, η) is taken to be the FP mass term, and he has shown that there are no discontinuities in the massless limit. However, the Hamiltonian corresponding to the general action (2), without any special tunings, is unbounded and contains a ghost-like mode because of the violation of a constraint corresponding to h = 0 in (1). This ghost mode is called the the Boulware-Deser (BD) ghost [4]. Although the BD ghost problem had been remaining during a long time, in 2010, the BD ghost problem had been solved by discovering the de Rham-Gabadadze-Tolley (dRGT) model [5,6]. In [5], de Rham and Gabadadze tried to tune the parameters in the massive gravity model (2) by requiring the ghost-freeness in the high-energy limit called the decoupling limit. And they concluded that the Lorentz covariant potential term V (g, η) consistent with this requirement is parametrized by only three free parameters (mass parameter, and other two dimensionless parameters) in D = 4. Although the analysis in [5] was perturbative in powers of h µν , in [6], the full-nonlinear form of the action was obtained. After that, the BD ghost-freeness of the dRGT model without taking any limits was proved by using the Hamilton analysis [7,8] with the Arnowitt-Deser-Misner (ADM) variables [10]. Furthermore, as an extension of the dRGT model, it was pointed out that the model where the flat metric η µν is replaced by an arbitrary fixed metric f µν called the fiducial metric is also ghost-free [9].
Let us summarize known results on the constraint analysis of the dRGT model. As we have mentioned, the first works in [5,6] argue the BD ghost-freeness in the decoupling limit by using the so called Stückelberg trick. And the proof of the full order model in [7][8][9] is based on the Hamilton analysis. Nowadays, in the metric formulation, there are a lot of works [11][12][13][14][15][16][17][18][19] based on the Stückelberg trick and/or the Hamilton analysis with the ADM variables. On the other hand, in the vielbein description of the dRGT model [20], there are some works based on the Lagrangian formulation [21,22]. These works derive the constraints as possible as keeping the "Lorentz covariance" (for f µν = η µν ). Thus. the analysis becomes quite tractable and realize the Lorentz covariant expressions of the constraints in a certain parameter region. However, it has still been unknown how to derive the covariant expressions of the constraints (covariant constraints) in the metric formulation except for the linearized case [23][24][25].
Our purposes is to derive the constraints with the metric and the Lagrangian formulation, as possible as keeping the Lorentz covariance. It can be predicted that this approach could realize the covariant expressions of the constraints. The difficulty in the case of the metric formulation comes from how to treat the square root matrix S µ ν ≡ g −1 f µ ν which constructs the potential terms in the dRGT model. Fortunately, there have been some works for the model linearlized around a general background solution [23][24][25]. From the results of the linearlized model, we can guess which of the linear combinations of the EoM realize the covariant constraints in the full-nonlinear level. Although the method in [23][24][25] cannot straightforwardly be extended to the nonlinear case, we can avoid this difficulty by making a good use of the algebraic relation, for the square root matrix S µ ν = g −1 f µ ν . As the results, we derive the explicit forms of the covariant constraints in a certain parameter region. Furthermore, for the general parameter region, we prove the existence of the constraints corresponding to the constraints in (1). This paper is organized as follows: In Sec.II, we summarize the fundamental properties of the dRGT model. In Sec.III, we consider the model linearized around a general background solution and obtain some suggestions for deriving the covariant constraints of the nonlinear model. In Sec.IV, we derive the explicit form of the covariant constraints for a certain parameter region. In Sec.V and Sec.VI, we prove the existence of the constraints for the general case in any dimensions. In SecV, in order to make the understanding on the constraint structure clear, we demonstrate that the existence of the constraints becomes trivial if we admit an identity. In Sec.VI, we prove the identity used in Sec.V. The last section VII is devoted to the summary.

II. DRGT MODEL
Let us start with summarizing the fundamental properties of the dRGT model. In order to make the argument as possible as general, we consider the dRGT model with an arbitrary fiducial metric f µν , Here, the symmetric polynomial e (n) (S) is defined as follows, Hence, the EoM corresponding to the action (4) can be determined as follows, Here, we define the Einstein tensor by G µν ≡ R µν − Rg µν /2. As we have mentioned, in the system where the highest order time derivative terms in the EoM is second order, the constraints can be defined as the linear combinations of the EoM independent of any second order time derivatives of the dynamical metric, ∂ 0 ∂ 0 g µ1ν1 . One of the constraints is given by taking the divergence of the EoM (15) and using the Bianchi identity ∇ µ G µν = 0, Obviously, this equation is independent of ∂ 0 ∂ 0 g µ1ν1 . Hence, (16) can be regarded as a constraint corresponding to ∂ µ h µν − ∂ ν h = 0 in (1). We should note that the special tuning of the potential term in the dRGT model (4) is not essential for the existence of the constraint (16). In other words, for a model whose potential term is different from that of the dRGT model, such as V (g, η) in (2), the constraint corresponding to (16) exists. The special tuning of the potential terms in the dRGT model (4) is essential for the existence of an additional constraint corresponding to h = 0 in (1). The covariant expression of this constraint has not been obtained and it is the purpose in this paper.

III. SUGGESTIONS FROM LINEARIZED MODEL
The purpose of this section is to obtain a few suggestions for deriving the covariant constraints in the nonlinear model by investigating the covariant constraints of the model linearized around a general background solution. Just for simplicity, in this section, we consider only the case of β 0 = 0, β 1 = 0, β n = 0 (n = 2, 3, · · · , D − 1). Fortunately, for the linearized model, the linear combination expressing the additional constraint has been obtained in [23][24][25]. In the case of β 0 = 0, β 1 = 0, β n = 0 (n = 2, 3, · · · , D − 1), by denoting the linearized EoM as δE µ ν = 0, the following linear combination of δE µ ν does not contain any second order derivative of the perturbation δg µν , However, the method in [23][24][25] cannot straightforwardly be extended to the nonlinear case. Then, in this section, we show that Eq.(17) does not include ∂ µ ∂ ν δg µ1ν1 or equivalently, δφ 1 = 0 is a constraint, by a different way which can be extended to the nonlinear case.
In the case of β 0 = 0, β 1 = 0, β n = 0 (n = 2, 3, · · · , D − 1), the EoM of the nonlinear model given in (15) is expressed as follows, Here, we define S ≡ S ρ ρ . In this section, we regard this equation as a background equation. By denoting the metric perturbation as δg µν = h µν , the linearization of the EoM (15) around any background solution is given by, Here, we raise the index of the S µ ,αβ ρ by using g νρ as S µν,αβ ≡ g νρ S µ ,αβ ρ . We should note that the superscripts µν in the S µν,αβ is not symmetric with respect to the permutation µ ↔ ν although the matrix S µ ρ g ρν ≡ S νµ satisfies S µν = S (µν) . The antisymmetric part S [µν],µ1ν1 is determined as follows, from the linearization of the relation S µν = S νµ . We have used the relation (20) to obtain Eq. (19). On the other hand, the symmetric part S (µν),µ1ν1 is determined from the following relation, which is just the relation (12) expressed in terms of S µ ,αβ ν defined in (19). Here, we define the matrix multiplication [AB] µ ν ≡ A µ ρ B ρ ν for any matrices A µ ν , B µ ν . In [23,24], the explicit form of the symmetric part S (µν),µ1ν1 has been obtained as a function of S µ ν , g µν by solving the relation (21), and the constraint (17) is obtained by using the explicit form of the S µν,µ1ν1 . However, the explicit form of the S µν,µ1ν1 is not so simple and it is not obvious how to extend the linear analysis to the nonlinear case. Although more tractable method is demonstrated in [25] by redefining the massive spin two field h µν , this method cannot also be extended to the nonlinear case straightforwardly. Then, in this paper, we consider to show the existence of the constraint (17), by making a good use of the relation (21), without solving the relation (21) explicitly.
Let us investigate the constraint structure of the linearized equation (19). The divergence of Eq. (19), ∇ µ δE µ ρ = 0, denotes the vector constraints. We should note that the divergence of the linearized Einstein tensor ∇ µ δG µ ρ is not equal to zero for general backgrounds because of the non-commutativity of the covariant derivative ∇ µ and the variation operator δ. Indeed, for general backgrounds, the linearization of the Bianchi identity ∇ µ G µ ρ = 0 becomes, We should note that the second line of the above equation (22) becomes equal to zero only in the case of the Einstein manifold where the background metric g µν satisfies the Einstein equation with a cosmological constraint Λ because the Einstein tensor is proportional to metric, i.e. G µν = Λg µν . In other cases G µν = Λg µν , ∇ µ δG µ ρ is not equal to zero, ∇ µ δG µ ρ = 0, generally. Hence, by using the linearized Bianchi identity (22), we obtain, In the second line, we use the background equation (18) in order to eliminate the Einstein tensor G µν . Because this quantity (23) does not contain any second order derivative terms ∂ µ ∂ ν h µ1ν1 , under the EoM δE µ ν = 0, the relations ∇ µ δE µ ν = 0 can be regarded as constraints. Those are just the linearization of the vector constraints (16). Next, let us show that the linear combination (17) does not depend on any second order derivatives of h µν by using the relation (21). By multiplying the inverse matrix S −1 to the vector constraints (23), the quantity S −1νρ ∇ µ δE µ ρ in the linear combination (23) can be expressed as follows, In order to confirm that the linear combination (17) does not depend on ∂ µ ∂ ν h µ1ν1 , we have to express the coefficient tensor in front of the second order derivative term ∇ µ ∇ ν h µ1ν1 in ∇ ν (S −1ν ρ g σρ ∇ µ δE µ σ ) of Eq.(17) as in terms of S µν , g µν . Then, we have to deform the tensor S µν,µ1ν1 in (24) by using the relation (21). Let us show that it is not necessary to derive the explicit solution of (21) for S µν,µ1ν1 nor redefine the field h µν as done in [23][24][25]. There is a simple way showing that the linear combination (17) does not depend on ∂ µ ∂ ν h µ1ν1 . From the relation (21), we can easily show the following relations, We find that Eq.(25) determines the trace part S ρ ,µ1ν1 ρ in Eq.(24) as in terms of S µ ν . Furthermore, Eq.(26) determines the symmetrization of S −1ν ρ S (µρ),µ1ν1 with respect to µν. Although the antisymmetrization of S −1ν ρ S (µρ),µ1ν1 with respect to µν are not easily determined from the relation (26), this fact does not affect to our analysis as we will see below.
Then, by using the relations (25) and (26), Eq. (24) can be rewritten as follows, Here, the bracket [µ| · · · |ν] denotes the antisymmetrization with respect µν, e.g., We find the important fact that all the coefficient matrices in Eq.(27) except for g µν(µ1ν1) are antisymmetric with respect to µν. These terms do not affect to the constraint structure because the covariant derivative in these terms turn into the Riemann curvatures by taking the divergence of Eq. (27), i.e., Then, the only quantity with the second order derivative of h µν is g µν(µ1ν1) ∇ ν ∇ µ h µ1ν1 . However, we can cancel this term by the trace of the EoM (19), Therefore, the following linear combination does not contain any second derivatives of the massive spin-two field, Therefore, under the EoM δE µ ν = 0, the eqution δφ 1 = 0 can be regarded as a constraint.

IV. COVARIANT CONSTRAINTS
From this section, we investigate the additional constraint of the full-nonlinear model described by the EoM (15). In this section, by focusing on the case of β 0 = 0, β 1 = 0, β 2 = 0 , β n = 0 (n = 3, · · · D − 1), we obtain the Lorentz covariant expression of the additional constraint (for f µν = η µν ). As we will see in the following section V, this is the most general case where the covariant expression of the constraint exists.
Let us start with the case of β 0 = 0, β 1 = 0, β n = 0 (n = 2, 3, · · · D − 1) where the EoM (15) is given by, From the covariant constraint of the linearized model (30), we can easily predict that the following linear combination of E µν may be independent of any second order derivative terms ∂ µ ∂ ν g µ1ν1 , Indeed, by linearizing the above quantity (32) and using the background equation E µν = 0, we obtain the quantity (30). In this part, we show that the quantity (32) does not contain any second order derivative terms ∂ µ ∂ ν g µ1ν1 in the full-nonlinear level, i.e., we will see that the equation φ 1 = 0 (under the EoM E µν = 0) can be regarded as a covariant constraint. Although we suffer again from the problem due to the square root matrix, we can avoid this problem by using a method analogous to the linear case.
In the case of β 0 = 0, β 1 = 0, β n = 0(n = 2, 3, · · · D − 1), the quantity S −1ν ρ ∇ µ E µρ in (32) is given by, 1 The definitions of the covariant derivative of the "matrices" S µ1ν1 and f µ1ν1 are not changed from the definition of a usual second rank tensor, i.e., We should note that the covariant derivative of the square root matrix, ∇ µ S µ1ν1 , depends on the first derivative ∂ λ g µν through not only the Levi-Civita connection but also the partial derivative of the square root matrix, ∂ µ S µ1ν1 . If we substitute the expression (33) into the quantity (32) straightforwardly, we find that there is the second order covariant derivative ∇ µ ∇ ν S µ1ν1 which contains second order partial derivative ∂ µ ∂ ν S µ1ν1 . However, the dependence of ∂ µ ∂ ν S µ1ν1 on ∂ µ ∂ ν g µ1ν1 is not obvious and we cannot identify the dependence of (32) on ∂ µ ∂ ν g µ1ν1 in this way.
On the other hand, the covariant derivative of the fiducial metric, ∇ µ f µ1ν1 , depends on the first derivative ∂ λ g µν only through the Levi-Civita connection because the partial derivative ∂ µ f µ1ν1 is independent of ∂ µ g µ1ν1 . Hence, if we can rewrite ∇ µ ∇ ν S µ1ν1 in terms of ∇ µ ∇ ν f µ1ν1 , we can identify the dependence of (32) on ∂ µ ∂ ν g µ1ν1 . Then, we try to express the quantity ∇ ν S −1νσ ∇ µ E µ σ in (32) so that the second order covariant derivatives are only acting on the fiducial metric f µν .
Let us eliminate the covariant derivative from the square root matrix in Eq.(33). From the relation, which can be derived by taking the covariant derivative of the relation (7), we can show the following relations, We find that the quantity ∇ µ S in (33) can perfectly be rewritten in terms of ∇ µ f αβ by using (36). On the other hand, from (37), we can determine the symmetric part S −1(ν ρ ∇ µ S µ)ρ of S −1ν ρ ∇ µ S µρ in (33). Hence, we obtain the following expression, Although there remains a covariant derivative acting on S in the last term in (38), we should note that this covariant derivative turn into curvatures by taking an additional divergence because of the anti-symmetry of the superscripts µν. Hence, we obtain, Here, Φ (1) is the term which does not depend on any second order derivatives ∂ µ ∂ ν g αβ . The other terms in (39) depend on second order derivatives ∂ µ ∂ ν g αβ . The terms with curvatures, which are the second and the third terms in (39), are contributions from the anti-symmetric part S We find that the only term depending on the second order covariant derivatives, the first term in (39), is expressed in terms of ∇ ν ∇ µ f µ1ν1 . Because ∂ µ ∂ ν f µ1ν1 in ∇ ν ∇ µ f µ1ν1 does not depend on ∂ µ ∂ ν g αβ , we can pick up the terms depending on the second order derivative ∂ µ ∂ ν g αβ from Eq.(39). The parts depending on the second order derivative ∂ µ ∂ ν g αβ are given by, Here, we define the lowered connection Γ ρ,µν ≡ g ρσ Γ σ µν . We find that the r.h.s of (41) and the r.h.s of (42) are canceled with each other in Eq.(39). Hence, the following quantity Ψ (1) does not depend on ∂ µ ∂ ν g αβ , Finally, by substituting (43) into (39), we obtain, Because the quantities Φ (1) , Ψ (1) are independent of ∂ µ ∂ ν g αβ , the only term depending on ∂ µ ∂ ν g αβ is R/2. It is obvious that R/2 can be canceled by the trace E µ µ , and we find that the following quantity does not depend on ∂ µ ∂ ν g µ1ν1 , which can be derived by using (44) and (31). In the case of f µν = η µν , this linear combination can be regarded as a Lorentz scalar function. Therefore, in the case of β 0 = 0, β 1 = 0, β n = 0 (n = 2, 3, · · · D − 1), the equation φ 1 = 0 (under the EoM E µν = 0) can be regarded as a covariant constraint.
B. β0 = 0, β1 = 0, β2 = 0 , βn = 0 (n = 3, · · · D − 1) case Next, we derive the covariant constraint for the case of β 0 = 0, β 1 = 0, β 2 = 0 , β n = 0 (n = 3, · · · D − 1) where the EoM (15) is given by, This is the most general case where the covariant expression of the additional constraint is possible, because, in the following section V, we will see that there are no covariant form of the additional constraint in other cases. Let us calculate the following quantity, as we have done in the previous part IV A. The β 1 -term ∇ ν S −1ν ρ ∇ µ Y µρ (1) (S) in the above equation (47) have been calculated in (44) of the previous part IV A. Then, we now rewrite the β 2 -term ∇ ν S −1ν ρ ∇ µ Y µρ (2) (S) so that the second order covariant derivatives are only acting on f µν . By using the expansion formula given in (A.11), (2) (S) can be expressed as, We find that the second and the third terms of the right hand side of the above equation (49) can be rewritten by using Eqs.(36) and (37) given in the previous part IV A. Furthermore, the forth term of the right hand side of Eq.(49) is trivially rewritten as S −1ν ρ ∇ µ f µρ . On the other hand, the first term of the right hand side of Eq.(49) can be rewritten by using the following relation, which can be derived by using the identity obtained by multiplying S −1 Y (n) (S) ρ µ on the both sides of the relation (35) and use the commutativity [Y (n) (S), S] µν = 0, as we have shown Eq. (13). Hence, by using the identities (36), (37), and (50), we obtain, The divergence of the above relation (51) is given by, Here, we define the Φ (2) which does not depend on ∂ µ ∂ ν g αβ . Furthermore, we express the equation (52) by using the quantity Ψ (1) defined in (43). We have succeeded to express the second order covariant derivative terms in the quantity ∇ ν S −1ν ρ ∇ µ Y µρ (2) (S) by using ∇ µ ∇ ν f µ1ν1 . As we have done in Eqs. (41) and (42), by decomposing R and ∇ µ ∇ ν f µ1ν1 in Eq.(52) into the terms depending on ∂ µ ∂ ν g αβ and the other terms, the dependence of On the other hand, the dependence of S µν E µν on ∂ µ ∂ ν g αβ reads, Therefore, from the relations (44), (46), (47), (54) and (55), we find that the following linear combination of E µν does not depend on ∂ µ ∂ ν g αβ , Here, the quantity Ψ (2) does not depend on ∂ µ ∂ ν g αβ because of Eqs.(54) and (55). Therefore, in the case of β 0 = 0, β 1 = 0, β 2 = 0 , β n = 0 (n = 3, · · · D − 1), there is a covariant constraint φ cov = 0. As we will see in the following section V, this is the most general case where the covariant expression is possible.

V. CONSTRAINT STRUCTURE
From this section, we extend the argument in the previous section IV to the general case (β n = 0, n = 0, 1, · · · D−1). Although the covariant expression is impossible in the general case, we can show the existence of the additional constraint.
For the first step, in this section, the existence of the additional constraint can be shown by using the following identity, Although the general proof of this identity will be given in the following section VI, we can confirm that the identity (57) is valid for the case of n = 1, 2 from the results (44) and (54) in the previous section IV. By using the identity (57), we can easily show the existence of the additional constraint.

B. General case
Although the additional constraint cannot be expressed as a covariant form in the general case β k = 0 (3 ≤ k), we can easily show the existence of the additional constraint.
Here, we use the identity (A.8), and θ µν is the projection operator living in D − 1 dimensional space orthogonal to the time direction g 0 µ . We lower indexes of θ µν by using g µν and raise indexes by using g µν . We find that the quantity ∇ ν S −1ν ρ ∇ µ Y µρ (n) (S) does not depend on g 0µ g 0ν ∂ 0 ∂ 0 g µν nor g 0µ θ ρν ∂ 0 ∂ 0 g µν . In other words, the quantity ∇ ν S −1ν ρ ∇ µ Y µρ (n) (S) depends only on θ µ α θ ν β ∂ 0 ∂ 0 g µν . On the other hand, from (59), the second order "time" derivative terms in EoM (15) are expressed as, The EoM (64) can be solved with respect to the θ α µ θ β ν ∂ 0 ∂ 0 g αβ . In order to do that, we introduce an "inverse matrix" θ −1 µν,µ1ν1 of the operator θ (µν)(µ1ν1) so that it satisfies the following condition, These conditions are uniquely solved as, Using the above operator (66), we obtain, From the identity (63) and the above relation (67), we find that the following linear combination of E µν does not depend on any second time derivative ∂ 0 ∂ 0 g µν , or, Hence, there is the additional constraint ψ = 0 or ψ ′ = 0. In the case of f µν = η µν , this combination (68) (or (69)) is not Lorentz invariant but invariant for spacial rotations. This result is consistent with known several results. [25] has concluded that the covariant constraint in the linearized model with metric formulation for D = 4 exists only in the case of β 3 = 0. [21] have argued the covariant constraints in the vielvein formulation for D = 4 and concluded that the covariant constraint exists only in the case of β 3 = 0. These results are consistent with our result that the covariant constraint exists only in the case of β 0 = 0, β 1 = 0, β 2 = 0 , β n = 0 (n = 3, · · · D − 1).

VI. PROOF OF IDENTITY
In the previous section V, we have seen that the identity (57) plays a crucial role for the existence of the additional constraint. In this section, we complete our argument by proving the identity (57) for the general case.
Let us explain how to show the identity (57). By using the expansion formula, whose derivation is given in (A.11) of the Appendix A, we can express the quantity S −1ν ρ ∇ µ Y µρ (n) (S) in the left hand side of the identity (57) as follows, As we have done in the Sec.IV, in order to confirm the identity (57), we try to rewrite ∇ µ ∇ ν S αβ in ∇ ν S −1νσ ∇ µ Y (1) µ σ in terms of ∇ µ ∇ ν f αβ by using the relation, By multiplying the matrix S n (n ∈ Z) on the above equation (72) from the left hand side, and multiplying the matrix S m g (m ∈ Z) from the right hand side, we obtain more general formula, Here, we define the matrix multiplication [AB] µν ≡ A µ ρ B ρν , and we omit the metric g in the definition of the matrix multiplication. In particular, in the case of n = m, the left hand side of the above equation (73) can be expressed as a symmetrization of a matrix product, As we have derived Eq.(50), we obtain the following formula from Eq.(73), Furthermore, we can rewrite the quantity S −1ν ρ ∇ µ S 2n+1 µρ (n = 0, 1, 2, 3 · · · ) from the relation (74), whose derivations are given in Appendix B. Here, we separate the set of n into the even case n = 2m (m = 0, 1, 2, · · · ) and the odd case n = 2m+1 (m = 0, 1, 2, 3 · · · ). In the above equations (75), (76) and (77), all the covariant derivative terms except for the anti-symmetric terms with respect to µν in (76) and (77) are expressed in terms of ∇ λ f µν . This fact means that we can rewrite all the terms with ∂ µ ∂ ν g αβ in the quantity ∇ ν S −1ν ρ ∇ µ Y µρ (n) (S) as in terms of ∇ µ ∇ ν f αβ and Riemann curvatures because the covariant derivative of S 2m ∇ µ SS 2m−1 [µν] and S 2m+1 ∇ µ SS 2m [µν] in Eqs.(76) and (77) will turn into curvatures by taking an additional divergence. Therefore, in principle, we can pick up the terms depending on ∂ µ ∂ ν g αβ by using the identities (75), (76) and (77) from ∇ ν S −1ν ρ ∇ µ Y µρ (n) (S) in the identity (57).
By substituting (75) into the divergence of (71), we obtain, Here, we define the following quantity, By considering the several cases K = 4m, 4m + 1, 4m + 2, 4m + 3 (m = 0, 1, 2, 3 · · · ) for T (K) , we obtain, whose derivation is summarized in Appendix B. Here, we use Eqs.(76) and (77) for calculating Eqs. (80) and (82), respectively. Furthermore, we define the notation "∼" which means the equivalence up to terms independent of ∂ µ ∂ ν g αβ , and define the bi-linear function, for any a, b ∈ Z. We should note that this tensor is not symmetric for the permutation µ ↔ ν, but satisfies the following relation,T Then, in the case of a = b, this tensor becomes symmetric with respect to µν,T µνµ1ν1 (S a , S a ) =T (µν)µ1ν1 (S a , S a ). Next, we identify the terms depending on ∂ µ ∂ ν g αβ in Eqs (80) . Thus, we can identify these dependences on ∂ µ ∂ ν g αβ by using the following relations, which can be shown by straightforward calculation. Here, we define the new bilinear function, which is just the symmetrization ofT µνµ1ν1 with respect to µν, i.e., T µνµ1ν1 =T (µν)µ1ν1 . Hence, T µνµ1ν1 (S n , S m ) satisfies the following identity, Furthermore, we should note the relationT µνµ1ν1 (S m , S m ) = T µνµ1ν1 (S m , S m ). By using the identities (86) and (88), the dependence of Eqs.(80), (81), (82) and (83) on ∂ µ ∂ ν g αβ can be identified as follows, We should note that all the above terms (89), (90), (91) and (92) can be summarized as follows, for K = 0, 1, 2, 3 · · · . By substituting the above relation (93) into (78), we obtain, (94) Finally, using the expansion formula (A.15) in Appendix A, we obtain the identity, (96)

VII. SUMMARY
We have shown the existence of an additional constraint of the dRGT model in the Lagrangian formulation and the metric formulation. We found an identity (57) which plays a crucial role for the existence of the additional constraint, and we have proved the identity for any dRGT potential terms. This identity realizes the covariant expression of the additional constraint (56) in a certain parameter region.
There is a possibility that the analysis in this paper may be extended to the Bimetric gravity [31], which is obtained by regarding f µν as a dynamical field and adding the Einstein-Hilbert term (4). In [25], the Bimetric gravity linearized around a general background solution has been considered. They have deriven the linear combination of the linearized EoM which does not contain any second order time derivative terms of δg 0µ and δf 0µ (Here, we denote the perturbation of the metrics g µν , f µν as δg µν , δf µν , respectively.). This linear combination in the Bimetric gravity corresponds to the identity (57). Therefore, by extending the analysis in this paper to the Bimetric gravity, we may develop the tractable method for deriving the constraints in the Bimetric gravity.
Expansion of Y µν (n) : Let us define the matrices Y µν (n) (S) as follows, Although the tensor S µν denotes the square root matrix given in (4), the following properties are valid for any symmetric matrices S µν = S (µν) . Y µν (n) (S) can be expanded in powers of [S n ] µν as follows, This relation can be obtained by solving the recursion relation, which can easily be shown by using the following expansion relation, : Furthermore, we derive a expansion formula for the quantity Y µνµ1ν1 (n−1) (S) defined as follows, Y µνµ1ν1 (n−1) (S) ≡ 1 (n − 1)! g µνµ1ν1µ2ν2···µnνn S µ2ν2 · · · S µnνn . (A.14) The symmetrized quantity Y µν(µ1ν1) (n) (S) can be expanded as follows,

Appendix B: Detailed calculations
In this appendix, we summarize the detailed calculation of Sec.VI. Derivation of (76) and (77): The derivations of Eqs.(76) and (77) are given by, The coefficient operator of the first term can be decomposed as follows, In the first line, we separate the set of k into even numbers k =2k ′ and odd numbers k =2k ′ + 1. In the second line, we separate the region 0 ≤ k ′ ≤ 2m − 1 into 0 ≤ k ′ ≤ m − 1 and m ≤ k ′ ≤ 2m − 1, and redefine the indexes. On the other hand, by using the formula (76), the second term of (B.3) can be expressed as follows, Here, we use the relation "∼" which means the equivalence up to terms without ∂ µ ∂ ν g αβ . By substituting (B.4) and (B.5) into (B.3), we obtain (80). As a similar way, we can derive (82) by using (77). Derivation of (81) and (83): Let us derive Eq.(81). In the case of K = 4m + 1, (79) becomes, As we have derived Eq.(B.4), the first term of the right hand side of the above equation (B.6) can be deformed as, On the other hand, the second term of the right hand side of Eq.(B.6) can be expressed as, In the second line, we separated the region 0 ≤ k ≤ 2m into 0 ≤ k ≤ m − 1 and m ≤ k ≤ 2m, and redefined the indexes. By substituting (B.7) and (B.8) into (B.6), we obtain Eq.(81). As a similar way, we obtain Eq.(83).