Perturbations of Cosmological and Black Hole Solutions in Massive gravity and Bi-gravity

We investigate perturbations of a class of spherically symmetric solutions in massive gravity and bi-gravity. The background equations of motion for the particular class of solutions we are interested in reduce to a set of the Einstein equations with a cosmological constant. Thus, the solutions in this class include all the spherically symmetric solutions in general relativity, such as the Friedmann-Lema\^{i}tre-Robertson-Walker solution and the Schwarzschild (-de Sitter) solution, though the one-parameter family of two parameters of the theory admits such a class of solutions. We find that the equations of motion for the perturbations of this class of solutions also reduce to the perturbed Einstein equations at first and second order. Therefore, the stability of the solutions coincides with that of the corresponding solutions in general relativity. In particular, these solutions do not suffer from non-linear instabilities which often appear in the other cosmological solutions in massive gravity and bi-gravity.


Introduction
Massive gravity is one of the potent candidates for modified theory of general relativity. As early as in 1939, a linear theory of massive gravity was proposed by Fierz and Pauli (FP) [1]. In order to avoid the inconsistency on massless limit of this theory [2,3], a nonlinear extension of the FP theory was considered [4]. Boulware and Deser, however, found that the nonlinear theory simply extended from the FP theory contains an unphysical ghost degree of freedom (BD ghost) [5]. Because of this ghost problem, a healthy theory of non-linear massive gravity had not been established for a long time.
Recently, de Rham, Gabadadze, and Tolley (dRGT) proposed a mass potential which can remove the BD ghost mode in a decoupling limit [6,7], and Hassan and Rosen have finally proven that the dRGT massive gravity theory is free from the BD ghost without taking the decoupling limit [8][9][10]. The dRGT theory of massive gravity has three parameters, graviton mass m and coupling constants of nonlinear self interactions, α 3 , α 4 . Complementary approaches of the BD ghost problem are studied in refs. [11][12][13][14].
Massive gravity includes a non-dynamical tensor field called fiducial metric, f µν , in order to construct a mass potential. For example, the FP and the original dRGT theories are constructed by adopting the Minkowski fiducial metric. Hassan and Rosen proposed an extended theory of dRGT massive gravity with a fiducial metric being dynamical as well by introducing the Einstein-Hilbert term of the fiducial metric in the action. They proved that this theory is also free from the BD ghost [15]. Since this theory contains two symmetric dynamical tensor fields of metrics, it is called bi-gravity theory.
The tests of massive gravity and bi-gravity using cosmological and black hole solutions have been explored intensively. In dRGT massive gravity, several types of exact, homogeneous, and isotropic solutions have been known so far. One example is the open Friedmann-Lemaître-Robertson-Walker (FLRW) solution with the flat fiducial metric in the FLRW slice [16]. The second-order perturbations of this solution show nonlinear instability and hence this solution is not viable unfortunately [17][18][19]. It should be noted that similar solutions with an anisotropic fiducial metric are known to be stable [20][21][22]. Another example of cosmological solutions is that with a fiducial metric which is flat but expressed in terms of nontrivial coordinates. This class of solutions can be divided into two types. The first type includes the solutions found in refs. [23][24][25], which exist for the whole parameter region of α 3 and α 4 , while the other type includes the solutions found in refs. [26,27], which exists only for a one parameter family in the parameter space (α 3 , α 4 ). Though the perturbations of the former type of solutions have already been studied in refs. [28][29][30][31], the perturbations of the latter type of solutions have not yet been investigated. Therefore, in this paper, we focus on the latter type of solutions. It should be noted that cosmological solutions with non-flat fiducial metrics are also studied in refs. [32][33][34][35].
A lot of static and spherically symmetric solutions have also been found up to now. The classification of such spherically symmetric solutions is studied in refs. [39,64,65]. The exact Schwarzschild (-de Sitter) solutions are classified to the following three classes. The first class is a solution with diagonal metric tensors [66][67][68][69] and linear perturbations of this class of solutions are studied in ref. [69][70][71][72][73] in the framework of both massive gravity and bi-gravity. The second class of solutions is a solution with a off-diagonal metric tensor and arbitrary α 3 and α 4 [74], and the perturbation of this solution is studied in ref. [69,75]. The last class is a solution with a off-diagonal metric tensor and a special choice of the parameters α 3 and α 4 [76][77][78], where linear perturbations have been studied only in massive gravity with a flat fiducial metric [79] and not in bi-gravity.
In the present article, we will give a unified and general treatment for solutions with an off-diagonal metric tensor in massive gravity and bi-gravity belonging to a one parameter family of α 3 and α 4 , which include both cosmological [26,27] and spherically symmetric black hole solutions [64,[76][77][78][79]. We will find that the equations of motion for this class of solutions exactly reduce to those of general relativity (GR) with a cosmological constant not only at the background and linear (first-order) perturbation level but also at the level of quadratic (second-order) perturbations. This result shows that massive gravity and bi-gravity can allow any spherically symmetric solution of GR including its stability, the evolution of linear perturbations, and the backreaction from linear perturbations, while it simultaneously implies that one cannot distinguish massive gravity or bi-gravity from GR by using spherically symmetric solutions and their perturbations at least up to quadratic order.
Our paper is organized as follows. In the next section, we briefly review the theory of bi-gravity (and massive gravity as a trivial case of a fixed fiducial metric) and derive the equations of motion in a general setting. In section 3, we derive a generic non-diagonal spherically symmetric background solution. Then, we investigate linear perturbations around those background solutions in section 4. There we will see that the terms coming from the mass potential must vanish in order to satisfy the Bianchi identity. In section 5, we investigate higher order perturbations and find that the same results as the linear perturbations apply for the quadratic perturbations. The final section is devoted to summary and discussion. Some details will be given in the appendices.

Review of bi-gravity
In this section, we give a brief review of bi-gravity. Bi-gravity is a theory consisting of two dynamical tensor fields, g µν and f µν , called physical and fiducial metrics, respectively. Massive gravity can be understood as a special case where the fiducial metric is fixed and non-dynamical. Its action is given by the Einstein-Hilbert term for each metric with the interaction term S mass and matter actions: is the Ricci scalar and κ represents the ratio of the effective Planck masses for g µν and f µν . In the case of κ = 0, the tensor field f µν does not have its kinetic term and hence is non-dynamical. This case corresponds to the massive gravity theory originally proposed by de Rham, Gabadadze, and Tolley [6,7]. S matter [g] and S matter [f ] are the matter actions coupled to g and f , respectively, Here we implicitly assume the matter actions possess the two general covariance with respect to g µν and f µν separately, though the full theory does not have such a symmetry. It should be noted that another type of matter coupling, which does not possess the two general covariance, is also studied by [56][57][58][59][60][61][62][63]. The energy-momentum tensors coming from S matter [g] and S matter [f ] are defined as Due to the two general covariance of matter actions we assumed, both energy-momentum tensors are conserved, that is, ∇ µ are the covariant derivatives with respect to g µν and f µν , respectively. Hereafter, we will omit the suffixes g and f when no confusion is expected. Now the interaction term S mass is tuned to be free from the BD ghost mode and given by where i-th order contributions e i (γ) are given by and m, β i being free parameters of the interaction term. Since they always appear in the combination m 2 β i , essentially there are five free parameters. The space of parameters corresponds to the one of the three parameters of the dRGT theory, m, α 3 , α 4 , and two cosmological constants, Λ (g) , Λ (f ) , and their relations are given by 1 14) Taking the variation of the action with respect to g µν and f µν , we will obtain the equations of motion for the two tensor fields. The equations of motion for g µν are given by The indices here are raised or lowered by g µν . The equations of motion for f µν are given by where G[f ] µ ν is the Einstein tensor constructed from f µν and The indices here are raised or lowered by f µν .

Bi-spherically symmetric background solutions
Here, we attempt to classify some of the spherically symmetric solutions in bi-gravity and identify those which obey the same equations of motion as in general relativity. These classes of solutions include the cosmological and black hole solutions known so far [26, 36-39, 64, 79]. Let us consider the following bi-spherically symmetric metrics: with dΩ 2 = dθ 2 + sin 2 θdφ 2 . The matrixḡ −1f takes the following form, and from these ansatz it is straightforward to see that the square root of the above matrix is of the formγ It should be emphasized that the following discussion does not rely on the concrete expressions of a(t, r), b(t, r), c(t, r), and d(t, r), but rather relies only on the fact thatγ µ ν is of the form of eq. (3.4).
As explained earlier, we are interested in the case where the equations of motion for both metrics reduce to the Einstein equations with cosmological constants at the background level. Therefore, in order for X (g) µ ν to be a cosmological term, the non-trivial off-diagonal components,X must vanish. We focus on the case of b(t, r) = 0 or c(t, r) = 0, and A = 1, since with b(t, r) = c(t, r) = 0, we will obtain a diagonal metrics as mentioned in Sec. 1 and the perturbations of such diagonal solutions have already been studied.
For non-diagonal solutions we are interested in, the condition that eqs. (3.5) and (3.6) vanish leads to Another requirement necessary forX (g) µ ν to be a cosmological term is where, we have defined C(t, r) as With C(t, r) = 0, at least three eigenvalues ofγ µ ν are equal to A. This class of solutions includes the cosmological solutions found in refs. [23][24][25] and the Schwarzschild solutions obtained in ref. [75]. The perturbations of those solutions have already been studied in detail in refs. [28][29][30][31].
In the present study, we therefore concentrate on the case with C(t, r) = 0. In this case, the solution to eq. (3.8) is Here we have assumed that α 3 = −1. Equations (3.7) and (3.10) are consistent provided that the parameters of the theory, α 3 and α 4 , satisfy This is equivalent to the condition that the two branches of the solution (3.7) degenerate. Thus, we see that only the particular one-parameter family of α 3 and α 4 satisfying (3.11) admits the class of solutions we are focusing on. Note that when eq. (3.11) is fulfilled A can also be expressed simply as A = −β 2 /β 3 . Note also that the range of α 4 is limited as α 4 = (α 3 + 1/2) 2 + 3/4 ≥ 3/4. In this one-parameter family of α 3 and α 4 with eq. (3.11), the interaction terms for bispherically symmetric metrics (3.1) and (3.2) with eq. (3.10) are of the form of a cosmological term. For g µν the interaction term gives In the case of dRGT massive gravity, f µν is not a dynamical but a fixed metric, and hence we need not consider the equations of motion for f µν . The class of solutions with C(t, r) = 0 includes the cosmological solutions [26,27], the black hole solutions [76][77][78][79], the Lemaître-Tolman-Bondi (LTB) solution [27], and the Reissner-Nordström (RN) solution [77]. Since the equations of motion for g µν (and, in fact, those for f µν as well) reduce to the Einstein equations with a cosmological constant, any spherically symmetric solution in GR is also a solution of the one-parameter subclass (3.11) of bi-gravity and massive gravity with a suitable fiducial metric. In appendix A, we present some examples of bi-FLRW and bi-Schwarzschild-de Sitter solutions belonging to this class.

Linear perturbations
Now, we analyze linear perturbations around bi-spherically symmetric solutions given in the previous section. The two tensor fields of metrics are perturbed as The first-order perturbation, δγ µ ν , of γ µ ν is defined as which can be written in terms of the metric perturbations by solving the following equations, where δg µ ν =ḡ µρ δg ρν and δf µ ν =ḡ µρ δf ρν . For our purpose we do not need the explicit form of the solution to the above equation, though it is obtained for a general fiducial metric in ref. [80]. Actually, without the explicit form of δγ µ ν , we can directly calculate X µ ν from eq. (4.3) with eq. (3.4) as Since the Einstein tensor satisfies the Bianchi identity ∇ (g) µ G µ ν [g] = 0 and the energymomentum tensor is conserved, the tensor X (g) µ ν also satisfies ∇ µ X (g) µ ν = 0. As demonstrated in appendix B, this requirement leads to a stronger condition which yields δγ a b = 0 for a, b = 2, 3. (Note that we are interested in the case with C(t, r) = 0.) Since eq. (4.6) also implies δX (f ) µ ν = 0. Thus, the equations of motion for the linear perturbations δg µν and δf µν reduce to the linearized Einstein equations.
In order to see the implications of the equation (4.6) in more detail, we express δγ µ ν in terms of the metric perturbations. Since only the angular components of δγ µ ν enter the equation (4.6), we only have to deal with the angular components of the equation (4.4), which can easily be solved becauseγ a b = Aδ a b for a, b = 2, 3. In fact, eq. (4.4) reduces to To sum up, the equations of motion for the first-order perturbations are equivalent to the following three equations: This is one of the main results of this investigation. The equations of motion for the perturbations of the two metrics coincide with the perturbed Einstein equations, though δg µν and δf µν are subject to eq. (4.11). Then let us count the number of graviton degrees of freedom for this perturbed system. Each symmetric tensor field of metric has ten components, and there are, respectively, four constraints (the Hamiltonian and momentum constraints) in eqs. (4.9) and (4.10), since those equations are the same as the perturbed Einstein equations. Furthermore, eq. (4.11) gives three constraints among the angular components of the perturbed metrics. We have four spacetime coordinates and hence there are four gauge degrees of freedom representing the choice of coordinates. In addition to those familiar gauge degrees of freedom, it turns out that there still remains another gauge transformation retaining the equations of motion (4.9), (4.10), and (4.11), as explicitly shown in appendix C. Note that this gauge degree of freedom corresponds to the ambiguity of the linear perturbations mentioned in ref. [79] for the Schwarzschild-de Sitter solution in the dRGT theory. Thus, the number of the remaining degrees of freedom is 10 × 2 − 4 × 2 − 3 − (4 + 1) = 4, which coincides with that of two massless gravitons. We can confirm that this is consistent with the result of the Hamiltonian analysis given in appendix D: there are ten first class constraints and twelve second class constraints, and hence there are 8 (= 40 − 10 × 2 − 12) degrees of freedom in phase space.
The above analysis can be applied to dRGT massive gravity only with eqs. (4.9) and (4.11) because the derivations of these equations do not depend on the equation of motion for f µν . Since, in this case, δf µν is composed of Stückelberg fields, the condition (4.11) just determines perturbations of Stückelberg fields. The remaining variables δg µν are governed by the Einstein equations and additional gauge symmetry appears as gauge degree of freedom for Stückelberg fields.

Second-order perturbations
In the previous section, we have shown that the first-order perturbations obey the perturbed Einstein equations and hence the behavior of the perturbations coincides with that of GR, though δg µν and δf µν are subject to eq. (4.11). One may then ask the question as to how one can discriminate this class of solutions in bi-gravity from the corresponding solutions in GR.
One possibility is to take into account the back reaction on the physical metric g µν from δf µν at second order. For this purpose, we incorporate second-order perturbations as follows: The perturbed metrics now give rise to the second-order perturbations of γ µ ν as whereγ µ ν is the background quantity defined in eq. (3.4) and δγ µ ν satisfies eq. (4.6), and hence δγ a b = 0 for a, b = 2, 3. The interaction term in the equations of motion for g µν can be calculated explicitly even at second order, and is given by This tensor satisfies the conditions assumed in appendix B, which, together with the Bianchi identity, yield Even at second order, leading to X (f ) (2)µ ν = 0 as well. These conditions provide the relation between g (2) ab and f (2) ab as follows: for a, b = 2, 3 and A, B, C = 0, 1. Thus, the metric perturbations obey the perturbed Einstein equations also at second order, and the number of graviton degrees of freedom coincides with that of two massless gravitons even at second order. This fact indicates that one cannot discriminate this class of solutions from the corresponding solutions in GR even at second order, unfortunately. On the other hand, this fact, fortunately, implies that our solutions are free from non-linear instabilities even in cubic action, which plague many cosmological solutions in massive gravity, such as the diagonal open FLRW solution [16][17][18], flat FLRW solution [24,28], and de Sitter solution [23,28].

Conclusions and discussion
In the present study we have investigated the perturbations of a class of spherically symmetric solutions in massive gravity and bi-gravity. First, we classified spherically symmetric solutions in massive gravity and bi-gravity and identified the specific class for which the background equations of motion are identical to a set of the Einstein equations with a cosmological constant. These solutions are allowed only with the one-parameter family of α 3 and α 4 satisfying eqs. (3.11). This class of solutions includes many known solutions, e.g., the FLRW solutions in ref. [26,27], the Schwarzschild(-de Sitter) solutions in ref. [76][77][78][79], the LTB solution in ref. [27], and the RN solution in ref. [77]. In fact, any spherically symmetric solution in GR is included in this class with a suitable choice of the fiducial metric f µν .
Next, we have investigated linear perturbations on this class of solutions. We have found that the interaction terms in the equations of motion for both metrics, δX (g) µ ν and δX (f ) µ ν , vanish thanks to the Bianchi identities, and hence the equations of motion reduce to eqs. (4.9)-(4.11), which are the perturbed Einstein equations with the relation (4.11).
We have also found that, in addition to the usual gauge symmetry associated with spacetime coordinate transformation, there is another gauge symmetry of the linear perturbations given by eqs. (C.8)-(C.11), which has already been known for the perturbations of the Schwarzschild de Sitter solution in dRGT massive gravity [79].
We have shown that the above result applies to second-order perturbations as well. Thus, one cannot distinguish this class of solutions in massive gravity and bi-gravity from the corresponding solutions of GR up to second order. This fact, however, implies that this class of solutions do not suffer from the non-linear instabilities, which often appear in the other cosmological solutions in massive gravity and bi-gravity. These aspects would suggest that massive gravity or bi-gravity with this one parameter family in (α 3 , α 4 ) may have additional fully non-linear symmetry, which may be responsible for the stability. Further investigations are necessary to clarify this point.
In this article, only spherically symmetric background solutions are discussed. So, it is an interesting and open question whether the results obtained in this article hold for more general background solutions. Our analysis on the background solutions can at least be applied to anyγ having the form of eq. (3.4) in any basis vectors, because the background equations of motion (3.12) can be obtained in an algebraic way from eq. (3.4) irrespective of a concrete expression forγ. Extending the above analysis to linear and non-linear perturbations of more general solutions, however, is a non-trivial issue simply because such analysis accompanies the derivatives. We will address these issues in a future publication.

A.1 Bi-cosmological solutions
First we consider a family of bi-FLRW solutions, in which physical metric takes the following FLRW form:ḡ µν dx µ dx ν = −dt 2 + a 2 (t) dr 2 1 − Kr 2 + r 2 dΩ 2 . (A.1) Comparing this metric with eq. (3.1) yields R(t, r) = a(t)r. We assume that f µν takes the same FLRW metric but in a coordinate (t(t, r),r(t, r), θ, φ) different from that of g µν , (A.5) In order to apply the results of the main body, the radial coordinater is determined to satisfy the following relation,r while the time coordinatet is arbitrary. In this case, the equations of motion for both metrics become Einstein equations with cosmological constants so that a(t) and b(t) obey the Friedmann equation with respect to each (cosmic) time, t ort. This kind of bi-FLRW solution becomes a slight generalization of that found in ref. [36], in which a specific choice of the coordinatet is adopted. For b = 1,K = 0, and Λ (f ) eff = 0, the fiducial metric f µν becomes the flat Minkowski one and hence this bi-cosmological solution includes that obtained in ref. [26,27] in dRGT massive gravity with the flat fiducial metric.

A.2 Bi-Schwarzschild de Sitter solutions
Our results are applied to the following bi-Schwarzschild de Sitter metrics as well: where r (g) andr (f ) represent Schwarzschild radii, and Λ

B Bianchi identity
In this appendix, we will show that a symmetric tensor satisfying a condition given below must vanish as long as it obeys Bianchi identity and the background metricḡ µν takes the matrix form of eq. (3.1).
Let us consider the following symmetric tensor X µν : where ǫ denotes the order of perturbations and Λ is a constant. The goal of this section is to show that X (n)µ ν vanishes if the Bianchi identity, ∇ µ X µ ν = 0, is imposed. The tensorḡ µρ X (n)ρ ν is symmetric because and both of X µν and Λg µν are symmetric. Then, from the property of the background metric g µν , it is characterized by three arbitrary functions as follows: or equivalently, On the other hand, from the eq. (B.1), the Bianchi identity reads where∇ µ is the covariant derivative with respect toḡ µν . Then, zero-th and first components of this equation are which yields the following solution when R is not a constant, X The remaining components of this equation are given by where we have used the relation (B.10). Removing X (n) 23 from these equations leads to the following equation for X (n) 33 : Since this is just the Laplace equation on a sphere, its solution is constant over the sphere: (B.14) By plugging this solution into eqs. (B.11) and (B.12), we obtain Thus, the solution of the Bianchi identity is given by (

C Additional gauge symmetry of linear perturbations
The linear perturbations have an additional gauge symmetry, which is combination of gauge transformation of g µν and f µν separately but keeping the equation (4.11). In this appendix, we will give a concrete form of such coordinate transformation. For this purpose, let us consider infinitesimal gauge transformation generated by x µ → x µ − ξ µ for g µν and x µ → x µ − (ξ µ + δξ µ ) for f µν 2 . We denote the difference A 2 δg ab − δf ab in eq. (4.11) under this transformation by ∆ ab , that is, To determine the gauge transformation, one establish a bi-tangent bundle T 2 M i.e. a fibre bundle locally isomorphic to M × T The additional gauge symmetry is characterized by ∆ ab = 0. The (2, 2) component of this condition is given by The remaining (2,3) and (3,3) components are given by where we have used eq. (C.2). One can easily find, similarly to eq. (B.13), that these equations reduce to the Laplace equation on a sphere: whose solution becomes To sum up, this additional gauge symmetry is characterized by Ξ(t, r, θ, φ), P (t, r), Q(t, r) as δξ 0 = Ξ(t, r, θ, φ), (C.8) r)Ξ(t, r, θ, φ) + R(t, r)Q(t, r) cos θ ∂ r R(t, r) , (C.9) δξ 2 = Q(t, r) sin θ, (C.10) One may regard R(t, r) itself as a radial coordinate and, in the new coordinates (t, R, θ, φ), the above transformation (C.8)-(C.11) with P (t, r) = Q(t, r) = 0 simply reduces to the transformation of the time coordinate. We can directly observe this symmetry in the action. Actually, the quadratic action of the mass term for the linear perturbations becomes where C(t, r) is the function defined in eq. (3.9), Λ

D Hamiltonian analysis of linear perturbations
We will count the number of graviton degrees of freedom of linear perturbations by means of the Hamiltonian analysis. So, we omit the matter action in this appendix. For this purpose, it is useful to decompose the perturbations in terms of spherical harmonics Y m l as done in ref. [81]. Due to the spherical symmetry of the background metrics, the modes with different eigenvalues of rotation (l, m) or parity (odd or even) develop independently, and the dynamics of each mode does not depend on m. Hence, we may suppose that m is equal to zero, without loss of generality.

D.1 Odd mode perturbations
Non-vanishing components of the odd mode perturbations with m = 0 are given by and where P l is the Legendre polynomial. In this subsection, hereafter, we omit the suffix l and the summation with respect to l for brevity. From the perturbed Einstein-Hilbert action with the mass term (C.12), the conjugate momenta of h where λ := (l − 1)(l + 2),M 2 pl := l(1+l) 1+2l M 2 pl π, and √ −ḡ represents the determinant of only 0, 1 components: We schematically decompose the Hamiltonian density as follows: where H odd GR,(g/f ) represents the contribution from each Einstein-Hilbert term and the effective cosmological term, which is the second term (or third term) in the right hand side of eq. (C.12). H odd mass represents the contribution from the first term in the right hand side of eq. (C.12). This decomposition is justified because S mass does not include time derivative of g µν and f µν . From the expression of the action (C.12), H odd mass is explicitly calculated as It should be noted that, for l = 1 mode, H odd mass vanishes, which implies that dynamics of l = 1 mode coincides with that of GR. Therefore, there should be additional gauge symmetry, under which each metric transforms independently. This transformation, actually, corresponds to the arbitrary function P (t, r) in eq. (C.11).
From now on, we focus on l ≥ 2 modes. The Hamiltonian density from the Einstein-Hilbert term is calculated as are some functions of t, r and C (g) is given by The primary constraints of this system are and then, the total Hamiltonian is Time evolution of the primary constraints is given bẏ ], (D. 23) which generate the following two secondary constraints, Time evolution of C (g) is given bẏ and that of C These equations impose another constraint, From time evolution of C (2) , we obtain yet another constraint, determines the combination of multipliers, A 2 v (g) − v (f ) . Then, no further constraints are generated.
Since one multiplier remains undetermined, one can easily find that there is gauge symmetry in this system. More explicitly, one can confirm that there are two first class constraints (and four second class constraints) in this system through the presence of two zero eigenvalues of 6 × 6 matrix {C I , C J }, where C I represent all of the six constraints. These two gauge symmetries correspond to the ones which the theory originally possesses. To summarize, the number of graviton degrees of freedom in this system is and completely coincides with the case of two massless gravitons. For the l = 1 mode, there are four variables (eight variables in phase space), h 0 (g/f ) , h 1 (g/f ) . As mentioned above, the interaction term S mass vanishes for l = 1 mode, and hence the action reduces to decoupled two Einstein-Hilbert action. Then, there are four first class constraints and four gauge symmetries which correspond to the general covariance of g µν and f µν separately. These four gauge symmetries can be arranged into the ones of the full theory and the additional ones described by P (t, r) in eq.(C.11). To summarize, the number of degrees of freedom of the odd l = 1 mode is , H (g/f ) 5 , H (g/f ) 6 , P (g/f ) 3 , P (g/f ) 4 , P (g/f ) 5 , P . On the other hand, H even mass is given by The time evolutions of these constraints are given bẏ +linear terms of H (g/f ) 4 , H (g/f ) 5 , H (g/f ) 6 , P (g/f ) 3 , P (g/f ) 5 , P , time development of these constraints only determines two of the multipliers v I (g/f ) and hence no more constraint appears. One can see that four of the multipliers v I (g/f ) remain undetermined, which implies that this system has corresponding gauge symmetry. Concrete calculation shows that this system has eight first class constraints (and eight second class constraints) through eight non-zero eigenvalues of 16 × 16 matrix {C I , C J }, where C I represent all of the sixteen constraints. These eight constraints are composed of six gauge symmetry of full theory and two additional symmetry described by Ξ in eqs. (C.8) and (C.9). To summarize, the number of graviton degrees of freedom for even modes can be estimated as 1 and six gauge degrees of freedom, that is, there are six first class constraints and four second class constraints. Four gauge degrees of freedom come from the ones of the full theory and two come from the additional ones described by Ξ in eqs. (C.8),(C.9). Then, the number of dynamical degrees of freedom is 1 and ten gauge degrees of freedom, that is, there are ten first class constraints and four second class constraints. Then, the number of dynamical degrees of freedom is It should be noted that six gauge symmetries correspond to the one of full theory, two gauge symmetries correspond to Ξ in eqs. (C.8),(C.9), and the other two gauge symmetries correspond to Q(t, r) in eq. (C.10).