Condensate of Massive Graviton and Dark Matter

We study coherently oscillating massive gravitons in the ghost-free bigravity theory. This coherent field can be interpreted as a condensate of the massive gravitons. We first define the effective energy-momentum tensor of the coherent massive gravitons in a curved spacetime. We then study the background dynamics of the universe and the cosmic structure formation including the effects of the coherent massive gravitons. We find that the condensate of the massive graviton behaves as a dark matter component of the universe. From the geometrical point of view the condensate is regarded as a spacetime anisotropy. Hence, in our scenario, dark matter is originated from the tiny deformation of the spacetime. We also discuss a production of the spacetime anisotropy and find that the extragalactic magnetic field of a primordial origin can yield a sufficient amount for dark matter.


I. INTRODUCTION
The existence of gravitational waves was indeed confirmed by the direct detections [1,2], and their quantum counterpart is called gravitons. The gravitons are defined by perturbations around a background spacetime. The effective energy-momentum tensor of the high-frequency gravitons in General Relativity (GR) was derived by Isaacson [3,4] which enables us to treat the gravitons as massless spin-2 particles whose energy and momentum change the background geometry. Due to the nonlinear features of the Einstein equations the effective energymomentum tensor cannot be straightforwardly defined. The gravitons are well-defined when their frequencies (and their momenta) are high enough compared with the curvature scale of the background and then the energymomentum tensor is defined via a non-local operation which projects the nonlinear quantities of the gravitons onto those in low-frequency modes. However, the low energy states of gravitons, i.e., low frequency/momentum modes of gravitons, should be ill-defined in GR. This is not the case when a graviton is massive.
Although GR is now widely accepted as a low-energy effective theory of gravity, the question whether the graviton is indeed massless or not has been long discussed (see [5][6][7] for reviews and [8][9][10] for experimental constraints on the graviton mass). The linear theory of the massive spin-2 field was constructed by Fierz and Pauli in 1939 [11]. Since the gravity must be represented by a nonlinear theory of the metric tensor, the Fierz-Pauli theory requires an extension to the nonlinear theory of the metric in order to obtain the theory of the massive graviton. Generic nonlinear extension of the Fierz-Pauli theory turns to be unstable, called the Boulware-Deser ghost [12]. However, the ghost-free nonlinear massive gravity was proposed by de Rham et al. in 2010 [13,14] which was further extended into the bigravity theory [15] and the multi-gravity theory [16]. In the bigravity theory or the multi-gravity theory the gravity is still a longrange force because there exists a massless graviton in addition to the massive graviton(s). In the present paper, we focus on the bigravity theory which contains a massless graviton and a massive graviton. The effective energy-momentum tensors of both massless and massive gravitons are defined in the similar way to the case of GR [17].
The bigravity theory has received much attentions related to the discovery of dark energy and dark matter. If the graviton mass is extremely small such as m ∼ 10 −33 eV, the present accelerating expansion of the Universe can be explained by the tiny graviton mass [18][19][20][21][22][23][24][25]. Other range of the mass may explain the origin of dark matter. For instance, dark matter is originated from a matter field in the "dark sector" when m 10 −27 eV [25,26] whereas the massive graviton itself is a candidate of dark matter when 10 −4 eV m 10 7 eV [17] (see also [27,28]).
The first suggestion to dark matter in the ghost-free bigravity theory was given by [22] which found that the anisotropy of the spacetime behaves like a dust fluid as for the contribution to the Friedmann equation. However, the following questions have not been cleared: Why does the anisotropy behave as a non-relativistic fluid? Whether or not can it explain other phenomena of dark matter, e.g., the cosmic structure formation? In the present paper, thus, we explore those questions and find that the dark matter component can be regarded as the "condensate" of the massive graviton and it can give local structures of the Universe.
We shall focus on the case when the massive graviton is dominated by the zero momentum mode; that is, the configuration of the massive graviton is almost homogeneous. This configuration can be interpreted as the condensate of the massive graviton which we call the massive graviton condensate. Contrary to the case of the massless graviton, the zero momentum mode of the massive gravi-ton shows a coherent oscillation due to the mass term. Therefore, we can define the energy-momentum tensor of the zero momentum mode of the massive graviton as long as the graviton mass is larger than the curvature scale of the background spacetime. We find that the zero momentum mode of the massive graviton gives a dark matter contribution to the Fridemann equation and the tiny fluctuations around the zero momentum mode provide the cosmic structure formation. The constraint on the graviton mass to be dark matter is the same as that obtained in [17], i.e., 10 −4 eV m 10 7 eV, in general.
From the geometrical aspect the zero momentum mode of the massive graviton represents the anisotropy of the spacetime. The universe filled with the zero momentum massive gravitons is interpreted as a homogeneous spacetime. The anisotropic component of the universe acts as dark matter which is indeed shown by [22]. The metric perturbations around the homogeneous spacetime can provide the structures of the universe.
The tiny anisotropy of the universe can be produced when there is a coherent field with an anisotropic stress. A possible candidate of the source is the extragalactic magnetic field of a primordial origin (see e.g., [29][30][31][32]). Recent blazar observations implies the existence of the extragalactic magnetic field whose lower bound of the strength B 0 is about 10 −17 G [33][34][35][36][37][38][39]. This magnetic field could be produced in the early universe [40,41]. We will show that the coherent magnetic field can yield a sufficient amount of the massive graviton condensate in order to explain the present abundance of dark matter.
The paper is organized as follows. After a brief introduction about the ghost-free bigravity theory in Sec. II, we define the effective energy-momentum tensor of the coherent massive graviton in Sec. III. The homogeneous configuration of ϕ µν is studied in Sec. IV which reproduces the result obtained by [22] from a field theoretical aspect. We then study the perturbations around the homogeneous mode in Sec. V. We show that the massive graviton condensate is indeed a viable candidate of dark matter. In Sec. VI, a production of the condensate to be dark matter is discussed. We give a summary and some discussions in Sec. VII. In Appendix A, we summarize the definitions of the energy-momentum tensors of the high-frequency massive and massless gravitons in a curved spacetime. We briefly study the Bianchi I universe in bigravity in Appendix B. In Appendix C, we detail the calculations about the inhomogeneous modes of the massive graviton condensate.

II. BIGRAVITY THEORY
The action of the bigravity theory proposed by Hassan and Rosen [15] is given by where g µν and f µν are two dynamical metrics, and R(g) and R(f ) are their Ricci scalars. The parameters κ 2 g = 8πG and κ 2 f = 8πG are the corresponding gravitational constants, while κ is defined by κ 2 = κ 2 g + κ 2 f . The ghost-free interaction term between the two metrics is given by where {b k } (k = 0 -4) are coupling constants and the 4×4 while U k are the elementary symmetric polynomials of the eigenvalues of the matrix g −1 f . Just for simplicity, we assume that matter is coupled only to the g-metric We shall briefly discuss the case when other types of matter fields are introduced in Sec. VII and Appendix A. Our conclusion is not changed even for those cases. The fully nonlinear equations of motion are given by where T µν is the matter energy-momentum tensor while T µν (int) and T µν (int) are derived by the variations of the interaction term U with respect to g µν and f µν , respectively. The contracted Bianchi identity and the matter conservation law ∇ µ are the covariant derivatives with respect to g µν and f µν , respectively.
There is a particular vacuum solution in which two spacetimes are homothetic such that where ξ 0 is a root of the quartic equation For the homothetic solutions, we obtain thus, the constants Λ g and Λ f are effective cosmological constants for the g-spacetime and the f -spacetime, respectively. In what follows, we assume because we are interested not in dark energy but in dark matter. The equations for the homothetic spacetime are exactly reduced into those in GR which indicates that the homothetic solution contains only the massless graviton modes. The degrees of freedom of the massive graviton mode do not exist in the homothetic solution.

III. ENERGY-MOMENTUM TENSOR OF COHERENT GRAVITONS
In this section, we derive the effective energymomentum tensor of the coherently oscillating gravitons focusing on the cosmological situation. General discussion about the energy-momentum tensor of gravitons is given in Appendix A.
As is well known in GR, when we discuss some structure produced by high frequency graviational waves, we have to separate the high frequency modes from smoothed background. The length or/and time scale associated with the gravitational waves should be sufficiently shorter than the typical scale of the smooth background [3,4]. Under this setting, the energy-momentum tensor of gravitational waves is defined by the nonlinear terms of the perturbed Einstein equation averaged over a length or/and time scale. We then obtain the propagating equation for the gravitational waves and the Einstein equation for the background including the backreaction from gravitational waves. We shall apply this procedure to the cosmological setting with the coherently oscillating massive gravitons. In the coherent case, we have to take care which we perform a spatial average or a time average. We consider the homogeneous universe with tiny metric perturbations represent the inhomogeneous perturbations. Since we are interested in the coherent gravitons, the time coordinate has to be appropriately chosen in order that the t = constant hypersurfaces are given by almost homogeneous spaces. Then, on each hypersurface, the homogeneous parts can be obtained by where · · · V is the spatial average where the averaged length scale is assumed to be much lager than the scale of the inhomogeneities. The dynamics of the homogeneous spacetime in bigravity was studied in [22]. Up to the linear perturbation theory, one may directly analyze the dynamics of the perturbations under the ansatz (3.1). In the present paper, however, we consider another separation of the metrics rather than (3.1). We first summarize the strategy of our calculations and the explicit analysis are given in Sec. IV and Sec. V. The bigravity theory contains two types of dynamical degrees of freedom, the massless graviton and the massive graviton. First, we separate the metrics g µν and f µν into the massless mode and the massive mode. Up to the linear perturbations around the homothetic background, we can introduce the mass eigenstates of the gravitons. However, the definitions of the massless mode and the massive mode of the metrics would be ambiguous in the nonlinear orders in which the gravitons are no longer diagonalized. (see discussions in [28,42]). Nevertheless, the massless mode and the massive mode are still meaningful if the perturbative expansion is viable. Therefore, we only consider the situation that the spacetimes are well approximated by the homothetic solution.
We focus on the late stage of the universe such that where H is the Hubble expansion rate in which the massive graviton has too heavy mass to be excited. Hence, the amplitude of the massive graviton is suppressed and then the metrics g µν and f µν are approximated by the homothetic solution (We recall that the homothetic solution give a spacetime without the excitation of the massive graviton). We then perturbatively treat the massive mode g (massive) µν which is defined by the difference between two metrics where α := ξ 2 0 κ 2 g /κ 2 f and we assume |g (massive) µν | ≪ 1. On the other hand, the massless mode is given by As a result, the metrics g µν and f µν can be decomposed into the massless mode and the massive mode as follows: We further decompose the massless and the massive modes into the low-frequency modes and the highfrequency modes, respectively: where the high-frequency modes h µν and ϕ µν are normalized by two mass scales where · · · T is the time average 1 over some time interval T which is assumed to be Then, the metric tensors g µν and f µν are divided into four components: (0) g µν , M µν , h µν and ϕ µν . The meaning of each variables is summarized in Table I.
We briefly mention the relation between two separations (3.1) and (3.7). The variables (0) g µν , M µν , h µν and ϕ µν are divided into the homogeneous parts and the inhomogeneous parts where the homogeneous parts are defined via the spatial average · · · V as with (3.2). We then obtain It is worth noting that because h µν is massless while ϕ µν is massive. The zero momentum mode of the massless graviton cannot be high-frequency whereas that of the massive graviton can be high-frequency due to the coherent oscillation. Since we have assumed that the configuration of the fields are almost homogeneous, the massive graviton is dominated by the zero momentum modeφ µν , that is, We call this configuration of ϕ µν the massive graviton condensate because a large fraction of ϕ µν occupies the single zero momentum stateφ µν . In the separation (3.1), the "backgrounds", i.e., the homogeneous modes g , are obtained by the spatial average whereas the "background" in (3.7), i.e., the low-frequency massless mode (0) g µν , is given by the time average and then it can be inhomogeneous. An advantage of the separation (3.7) is that the high-frequency "perturbations" h µν and ϕ µν can be treated as tensor fields, propagating on the low-frequency "background" (0) g µν , with well-defined energy-momentum tensors.
The amplitude of the massive graviton is small so we have the inequalities (3.16) The amplitude of h µν /M pl is also small since h µν is a part of the inhomogeneity. As a result, we have three small quantities M µν , h µν and ϕ µν which can be treated as the tensors with respect to the "background" metric (0) g µν . We adopt the notation such that the suffices on M µν , h µν and ϕ µν are raised and lowered by (0) g µν . However, the inequality (3.16) does not suggest that the backreaction of ϕ µν to (0) g µν is also small. The orders of magnitude of the Einstein tensor of (0) g µν and the energy-momentum tensor of ϕ µν , which we denote (0) G µν and T µν G , are estimated as where T µν G is explicitly defined by (3.22) below. Thus, if |ϕ µν |/M G ∼ H/m ≪ 1, the massive graviton ϕ µν can be a dominant component of the universe. In what follows, we assume the massive graviton is the dominant component.
Just for simplicity, we consider the case This is a specific case, but this assumption is reasonable for our interest since the massless gravitons, i.e., the To discuss the dynamics of the Universe, the effect of the massless gravitons can be ignored.
To discuss dynamics of the massive graviton condensate ϕ µν , it is sufficient to include the leading and subleading contributions associated with the adiabatic expansion in terms of m −1 . 2 Up to subleading order we can ignore M µν since the amplitude of M µν is suppressed by m −2 which is the sub-subleading order (see Appendix B). The low-frequency massive mode M µν gives only negligible contributions.
Ignoring h µν and M µν , the g-spacetime metric is given by The equations for (0) g µν is given by the time-averaged Einstein equation where the "effective" energy-momentum tensors of the matter (0) T µν and that of the massive gravitons T µν G are defined by the relation and The equation of motion for ϕ µν is where the effective graviton mass m eff is defined by and δ (2) E µν include the terms of quadratic in ϕ µν which is explicitly given by

25)
2 More precisely, we use the dimensionless parameter H/m for the adiabatic expansion. We just refer to m −1 as the order of the expansion. with The symbol · · · high denotes a high-frequency projection operator which is given by for a quantity X. The functionals δ (1) R µν and δ (2) R µν are the first order and the second order of the perturbed Ricci curvatures which are explicitly shown in Appendix A.
The amplitude of the coherent oscillation decreases due to the Hubble friction which finally cause the decreasing of the energy density of the massive graviton condensate. To solve (3.23), we have to retain terms of linear in first derivatives of the metric (0) g µν . 3 On the other hand, we may ignore terms of higher orders of derivatives of (0) g µν which are sub-subleading order contributions; thus, the covariant derivatives commute Note that the quadratic terms δ (2) E µν cannot be ignored. For the homogeneous ansatz, the Friedmann equation schematically reads when the massive graviton is the dominant component of the universe (see Sec. IV for the explicit expressions). The quadratic term in (3.23) is then which yields a comparable effect to the first derivative of the metric. Therefore, we should solve the nonlinear differential equation (3.23) to discuss the dynamics of the coherent massive graviton, in general. However, we assume the Z 2 symmetry for the self-interactions of the massive graviton: the interaction terms are invariant under the Z 2 transformation ϕ µν → −ϕ µν , it prohibits appearance of δ (2) E µν and then the basic equations become much simpler. The Z 2 symmetry is realized by supposing the symmetry of the gravitational action under the replacement which is realized when In this case, ξ 0 = 1 is always a solution to the equation (2.9). For the branch ξ 0 = 1, clearly from the definition of the massive mode, the symmetry (3.31) realizes the Z 2 symmetry of the massive graviton. Indeed, the parameters (3.32) yields α = 1, β 2 = 1/2 and then As a result, the equation for the massive mode is linear since the cubic terms can be ignored for our calculations 4 . By using the normalization of the mass parameter m, we can always set in which we obtain The effective energy-momentum tensor T µν is obtained from the smoothing of the true energy-momentum tensor T µν .
Even if we assume the true energymomentum tensor is conserved, i.e., ∇ µ T µν = 0, the smoothed energy-momentum tensor is not conserved, in general, since the energy of the matter can be converted to the one of the graviton and vice versa via the equation (3.23). The contracted Bianchi identity of (3.20) reads the smoothed total energy-momentum tensor is conserved: However, in the late stage of the universe, the massive gravitons must be decoupled from the matter due to the weakness of the gravitational interaction, i.e., The parameters α and β 2 do not appear at linear order except for the right-hand side of Eq. (3.23) (we note M G = α −1/2 M pl ). As a result, the existence of the Z 2 symmetry (3.32) does not change the theory up to the linear order except for the coupling strength to the matter.
Then, the energy-momentum tensors are individually conserved: The conservation of T µν G is directly confirmed by using the equation of motion. For the freely propagating gravitons (3.38), the equation (3.23) is reduced into Using these equations, one can find which is a sufficient condition on the conservation of the graviton energy-momentum tensor (3.40). We notice, however, that two conservations (3.40) and (3.43) are not equivalent since (3.40) reads that the smoothed quantity of T µν G is conserved. Eq. (3.40) has information only about macroscopic behavior of ϕ µν while Eq. (3.43) (or (3.41) and (3.42)) involves information about microscopic behavior.
In the following sections, we will show that there exists a solution such that ϕ µν behaves as dark matter which explains not only the background dynamics of the universe but also the structure formation. Since we have assumed the inhomogeneities are smaller than the homogeneous modes, we shall linearize the expressions in terms of the inhomogeneities. For instance, the graviton energy-momentum tensor is expressed as with, in order of magnitude,

IV. MASSIVE GRAVITON CONDENSATE AS DARK MATTER
In this section, we consider the homogeneous modē g µν andφ µν . We assume the flat Friedmann-Lemaître-Robertson-Walker (FLRW) background and a simple ansatz forφ µν where a andφ are functions of t. Note that the ansatz (4.2) trivially satisfies the constraints (3.42).
The massive graviton originally appears from the metric perturbations. The present set up (4.2) corresponds to considering the axisymmetric Bianchi type I universe in bigravity which we will detail in Appendix B (see [22] for more details). One may worry about thatḡ µν should be also given by the the axisymmetric Bianchi type I universe rather than FLRW universe (4.1). However, as we will see just below, the averaged graviton energy-momentum tensor is indeed isotropic and (4.2) is consistent with (4.1) (see also [43]). Furthermore, even if one replaces (4.1) with the Bianchi type universe, its anisotropy decreases as a −6 and then the anisotropic part inḡ µν is quickly ignored.
The equation (3.41) including up to the first derivatives of the metric reads where H =ȧ/a. By using (4.3), we find and other components are zero. The approximative so- whereφ 0 and θ 0 are integration constants. Since the initial phase θ 0 is not important for the discussion, we set θ 0 = 0. After averaging over the time interval T ≫ m −1 , the graviton energy-momentum tensor is calculated as where the energy of the massive graviton condensate is Therefore, the massive graviton condensate behaves as a dust fluid. This guarantees the ansatz (4.1).
To explain the abundance of dark matter, the amplitude ofφ is required to bē where H 0 is the present Hubble parameter. Since the physical spacetime is given by (3.19), the universe is filled with the coherent "gravitational waves" ϕ µν /M pl whose dimensionless amplitude and the frequency are Hz . (4.10) The oscillations have too small amplitude and too high frequency and thus there should be no constraint on the existence of ϕ µν at present.

V. COSMIC STRUCTURE FORMATION
To study the cosmic structure formation, we then introduce small inhomogeneity to the metric and the massive graviton ϕ µν . Note that the low-frequency "background" (0) g µν is not necessary to be homogeneous (see (3.11)). We treat (0) g µν including the inhomogeneity δg µν as the low-frequency massless mode which is verified as long as the momentum of the inhomogeneity is smaller than the graviton mass: where k is the comoving momentum of the inhomogeneity defined by (5.16) later.
For the calculations, we use the adiabatic expansion in terms of the graviton mass inverse m −1 (see Section 2.5 in [44] for calculations in the case of the scalar condensate and [45] for the vector case). We set the orders of both g µν (t) and δg µν (t, x) as O(m 0 ). Since the time derivatives acting on the low-frequency modes do not change the order of magnitude m −1 , i.e., ∂/∂t = O(m 0 ), the Friedmann equation leads toρ G = O(m 0 ). Hence, the homogeneous mode of the massive gravitonφ µν is of order O(m −1 ). On the other hand, the amplitude of the inhomogeneous mode of the massive graviton δϕ µν is of order O(m 0 ) as we will show later.
To evaluate the inhomogeneous parts of T µν G the coherent background including the sub-leading order is required which is given bȳ whereφ 1 ,φ 2 are slowly varying functions with To determine the explicit form ofφ 2 we need to solve the equation (4.3) up to the order O(m −1 ); however, the equation (4.3) is valid up to the order O(m 0 ). Nevertheless, the explicit form ofφ 2 is not necessary to evaluate δT µν G .
Since we take the adiabatic expansion in terms of not k/m but H/m, the spatial derivatives do or do not change the order m −1 depending on the scales of the inhomogeneities. For the large scales such that the spatial derivatives acting on the variables do not change the order of m −1 , i.e., in which the consistency of the Einstein equation leads to On the other hand, for the small scales the order of k is Then, the graviton energy-momentum tensor is evaluated as As we will show, this classification of the scales corresponds to the scales beyond or below the Jeans scale. The graviton mass should be m 10 −4 eV since we have not detected any deviations from the Newtonian gravitational law in the laboratory scales. 5 In that case the scale a/k ∼ (mH) −1/2 is quite small compared with the structures of the Universe. Therefore, the case (5.7) is irrelevant to the cosmic structure formation. Nevertheless, we shall discuss both scales (5.4) and (5.7), for completeness. Furthermore, the discussion about the small scale (5.7) will make the properties of the massive graviton condensate clear.
We have ignored sub-subleading terms to derive the equations (3.41) and (3.42). Supposing δϕ µν is of order O(m 0 ), these equation showing the order of errors explicitly are written as The corrections to the equations come from the approximations Eqs. (3.28) and (3.38) as well as the nonlinear quantities of ϕ µν . It is worth noting that, as long as the expressions are linear in inhomogeneities, the accuracy of the calculations are same in both cases (5.4) and (5.7). As for the graviton energy-momentum tensor, although we have ignored higher order corrections to δT µν G such as m 2φ3 δϕ, they are of order O(m −1 ). Up to the order O(m 0 ), the higher order terms are negligible.
In the standard cosmological perturbation theory, general perturbations can be decomposed into scalar type, vector type, and tensor type perturbations and they are decoupled at the linear order due to the spatial homogeneity and isotropy of the background. In the present case, the background configurationφ µν breaks the spatial rotational symmetry and then the different modes could couple in the scalar-vector-tensor decomposition (see [46] for perturbations around the anisotropic universe). However, there still exists the rotational symmetry in the y−z plane. The perturbations can be decomposed into the even parity perturbations and the odd parity perturbations associated with the rotation in the y − z plane: where explicit forms of perturbations are shown in Appendix C. Furthermore, due to the spatial translation symmetry of the background all variables can be transformed into the momentum space, e.g., and the different momentum modes do not couple.
Henceforth we use variables in the momentum space and the notation For the calculations, we decompose the perturbations into the odd parity perturbations and the even parity perturbations. However, to clarify the physical meaning of the results, we shall divide the odd parity perturbations and the even parity perturbations into the scalar type, the vector type, and the tensor type components. We define three dimensional harmonic scalar Y S , vectors and where ∂ 2 = ∂ i ∂ i and i, j = (x, y, z). The indices i, j are raised and lowered by δ ij and δ ij . The quantities Y S , Y i V , Y ij T are even parity quantities associated with the two dimensional rotation in the y − z plane while Y i V , Y ij T are odd parity quantities. The suffixes S, V , and T are attached to classify the quantities into the scalar, the vector, and the tensor types associated with the three dimensional rotation. We further introduce the quantities as then we obtain nine harmonics which are summarized in Table II TABLE II. The classifications of the even parity perturbations and the odd parity perturbations. scalar vector tensor Using the gauge condition (see Appendix C for details), the perturbations for the low-frequency mode δg µν can be given by We will calculate the perturbed graviton energymomentum tensor in the following subsections. In the momentum space, the perturbed graviton energymomentum tensor is given by the form from which we can read the energy, the velocity, the isotropic pressure, and the anisotropic stress regarding the massive graviton condensate as a fluid. The velocity for the even parity perturbation is decomposed into whereas the anisotropic stresses are decomposed into Since calculations for general perturbations are complicated, we just show particular solutions. The generic solutions are given in Appendix C.
First, the odd parity perturbations are not important for the structure formation since these modes contain only the vector type and the tensor type perturbations. We discuss only the even parity perturbations in this section. Furthermore, we find the vector type components are always decoupled up to the subleading order whereas the scalar type and the tensor type components are coupled in the small scales. Since the vector type perturbations represent the rotational modes and decay in time, the vector modes are not important. Hence, it is sufficient for the structure formation to consider the even parity perturbations without the contributions from the vector type components.
The irrotational solution for the even parity perturbations can be found under the ansatz where a, b = (y, z) whose indices are raised and lowered by δ ab and Y x and Y a are defined by (C3). Since the nondiagonal components of δϕ ij only have decaying modes, we just set δϕ ij = 0 for i = j. Note that we do not use the scalar-vector-tensor type harmonics to represent the components of the massive graviton in order that the expression is written in a clear form. Under the adiabatic expansion, the variables for the massive graviton can be given by We shall show the particular solutions for the large scales (5.4) and the small scales (5.7) in order.

A. Large scale inhomogeneity
In the large scales (5.4), the consistency of the equations leads to that the amplitudes of the variables have to be The order O(m −2 ) quantities give just sub-subleading contributions thus we can ignore them. Eqs. (5.11) and (5.12) yield the constraint equations which determine B 1 , C 1 and ψ 1 as Eq. (5.10) gives We do not find other equations within our accuracy. After using the above equation and taking the time average, we obtain δρ G = 6m 2 φ 1 δϕ 1 +φ 2 δϕ 2 +φ 2 1 2Ψ + k 4 k 4 H T , from the perturbed graviton energy-momentum tensor. Therefore, the graviton energy-momentum tensor is given by a form of pressureless perfect fluid without the vector type components (i.e., an irrotational fluid). Although the irrotational property is obtained because of our specific ansatz (5.31), the pressureless property is hold even if we consider the general solution.
Note that the evolution of δϕ 1 (and alsoφ 2 ) has not been determined within our accuracy. We can, however, determine the dynamics of δρ G by choosing the combinations δρ G as an independent variable instead of δϕ 1 . The dynamics of δρ G is determined by the conservation law of the averaged graviton energy-momentum tensor 6 which reads Although the definition of δρ G contains the tensor mode H T , Eq. (5.39) indicates that the dynamics of δρ G is independent from the tensor mode H T . Therefore, if we focus on only the macroscopic behaviour of the massive graviton condensate (i.e., we focus on only the dynamics of T µν G T ), the scalar modes and the tensor mode are decoupled. Needless to say, at the microscopic level, the scalar modes and the tensor modes should be coupled. For example, the dynamics of δϕ 1 would depend on the dynamics of the tensor mode H T as well as the scalar modes Φ, Ψ. Furthermore, the couplings between the scalar-vector-tensor modes would appear when we consider more higher order corrections. The present calculations are justified up to the sub-subleading order.
As a result, the massive graviton condensate behaves as a dust fluid in the large scales (5.4) even if a small inhomogeneity is introduced. The massive graviton condensate can cluster and then explain the cosmic structure formation. 6 We notice again that although the conservation law of the nonaveraged graviton energy-momentum tensor

B. Small scale inhomogeneity
Next, we discuss the solution in the small scales (5.7) in which the amplitudes are given by Φ, Ψ, H T , δϕ 1,2 = O(m 0 ) , The constraint equations (5.11) and (5.12) yield and whereas (5.10) gives Note that there are ten independent equations for ten independent variables {δϕ 1,2 , B 1,2 , C 1,2 , φ 1,2 , ψ 1,2 }. Hence, the dynamics of the massive graviton are completely determined within our accuracy, differently from the previous case. We notice that the equations (5.44) and (5.45) yield the Schrödinger equation in the cosmological background. The "wavefunction" u(t, x) is defined by the relation where δϕ(t, x) is the variable in the real space which is used only here. Then, we obtain where mΦu ≃ mΦφ 1 since we have considered the linearized theory. The wavefunction u is dominated by the coherent mode u 0 =φ 1 which suggests that the almost homogeneous configuration of ϕ µν represents the condensate of the massive graviton. A large fraction of massive gravitons occupies the state u 0 except for the tiny perturbations δϕ µν . Note that u 0 is also a solution to the Schrödinger equation since the gravitational potential Φ is zero for the homogeneous configuration.
We obtain and v V = π V = 0 .
where δ := δρ G /ρ G = 2δϕ 1 /φ 1 is the relative perturbation of the energy density. Clearly, in the scales beyond the Jeans scale, this equation admits a growing mode solution δ ∝ a due to the Jeans instability in the dust dominant universe. On the other hand, we find for k 2 ≫ k 2 J with a ∝ t 2/3 . The massive graviton condensate shows the acoustic oscillation.
The tensor mode is not decoupled from the scalar mode. The gravitational wave H T is sourced by the anisotropic stress π T which is related with the energy density of the massive graviton condensate. The acoustic oscillation of the massive graviton emits the gravitational waves.

VI. PRODUCTION OF MASSIVE GRAVITON CONDENSATE
As shown in the previous section, the massive graviton condensate is indeed a candidate of dark matter. We thus consider a production mechanism of the massive graviton condensate and discuss whether the massive graviton condensate can be the dominant component of dark matter.
To generate the massive graviton condensate we need an anisotropic source which is coherent on the cosmological scale 7 . A candidate is the cosmological scale magnetic field. The blazar observations suggests the lower bound of the strength of the extragalactic magnetic field B 0 is about 10 −17 G. The upper bound is obtained from the CMB observations which is about 10 −9 G [47]. Since the production and the evolution of the cosmological scale magnetic field are subject to discussion (see [48][49][50] for reviews), we consider the simplest scenario.
We assume that the coherent magnetic field is generated in the early universe, e.g., in the inflationary regime of the universe. Here, we do not discuss the initial spectrum of the dark matter perturbations which should depend on the details of the coherent magnetic field. We just estimate the produced amount of the coherent massive gravitons. If the magnetic field adiabatically evolves, the energy density decreases as a −4 . On the other hand, the energy density of the massive graviton decreases as a −3 . Therefore, the production of the massive graviton condensate by the magnetic field can be ignored in the late stage of the universe (m ≫ H). The separation (3.7) is not justified in H m. To discuss the early universe we have to directly analyze the Bianchi universe which is summarized in Appendix B. By solving Eqs. (B9)-(B13), a typical behavior ofφ produced by the coherent magnetic field is shown in Fig. 1.
The amplitude of ϕ µν does not grow in H ≫ m since the Hubble friction is too large. As a result, the dominant production of the massive graviton condensate should occur at H ∼ m. The produced amplitude of the condensate is estimated as where the asterisk represents the quantities at the production and B is the strength of the magnetic field. Ω B and Ω r are the present density parameters of the coherent magnetic field and radiation, respectively. To obtain the final expression we use the universe is dominated by radiation at the production, M 2 pl H 2 * ∼ ρ * r , and both B 2 and ρ r decrease as a −4 . Once the gravitons are produced, ϕ µν is decoupled from the magnetic field because the contribution from the interaction quickly decreases (see Eq. (B29)). In order that ϕ µν is the dominant component of dark matter, the present amplitude has to bē ϕ 0 ∼ M pl H 0 /m where the present and the produced amplitudes are related by a where z eq is the redshift of the equality time. The graviton mass should be 10 −4 eV m 10 7 eV where the lower bound is obtained from laboratory-scale experiments of gravity while the upper bound is given by the lifetime of the massive graviton. Hence, a consistent scenario is constructed when the present magnetic field is 10 −12 G B 0 10 −10 G , (6.3) which is indeed a viable region of the coherent magnetic field.

VII. SUMMARY AND DISCUSSIONS
We provide a scenario in which a tiny deformation of the spacetime is dark matter in the ghost-free bigravity theory. This deformation is interpreted as the "condensate" of the massive gravitons. Differently from the case of the massless graviton, the zero momentum state of the massive graviton is well-defined when m 2 ≫ H 2 . We find that the zero momentum massive graviton with small fluctuations is a viable candidate of dark matter.
We have also studied a production mechanism of the coherent massive gravitons with the mass range 10 −4 eV m 10 7 eV and shown the coherent magnetic field with 10 −12 G B 0 10 −10 G yields a sufficient amount of massive gravitons for dark matter. When the present value of the coherent magnetic field is determined by a future observation, we can fix a suitable value of the graviton mass to be dark matter.
Although we discussed the magnetic field as a source of the massive graviton condensate, another source to produce the condensate could exist. In general, if there exits a coherent anisotropic stress π coh whose coherent scale is L m −1 and the density is π coh /ρ r ∼ 10 −10 in the age H ∼ m, the massive graviton condensate is produced and becomes dark matter. Since the gravitons universally couple to matter fields, the source is not necessary to be a standard model particle. Any matter field can be a source of the gravitons.
If the anisotropic stress is a random field instead of the coherent field, the stochastic massive gravitons are produced which have been discussed in [17]. Even for the stochastic case, we obtain a viable scenario of the massive graviton dark matter. Hence, if the anisotropic stress existed in the early universe, it inevitably yields the stochastic or the coherent massive gravitons and then the massive gravitons can be dark matter.
One may interpret that the Z 2 symmetry is a fundamental symmetry of the massive graviton. Although we have assumed the Z 2 symmetry only for the selfinteractions of gravitons, the symmetry can be introduced into the matter-graviton interactions as well. The Z 2 symmetry holds if the action is invariant under the replacement g ↔ f in bigravity. At low energy scales, matter fields can couple to both metrics g µν and f µν via an effective composite metric [51][52][53][54][55][56][57][58][59]. Hence, when all matter fields couple to the composite metric g eff µν defined by g eff µν = the matter-massive graviton interactions respect the Z 2 symmetry. Indeed, expanding the metrics under (3.7) we obtain thus, the matter action is manifestly invariant under ϕ µν → −ϕ µν . In this case, the Yukawa interaction of the massive graviton does not appear and then the present graviton mass constraints cannot be applied. The details about the Z 2 symmetric bigravity theory are under investigation.
We have introduced the Z 2 symmetry to the selfinteractions of the massive graviton in order to simplify the calculations. However, the massive graviton condensate can be dark matter even without the Z 2 symmetry because the nonlinear terms are always sub-leading contributions and then may not affect the dynamics at leading order. The leading order expression of T µν G would be unchanged. In order that the massive graviton condensate is dark matter, the Z 2 symmetry of the selfinteractions should not be required.
Our scenario can be directly confirmed when we observe the coherent anisotropic oscillation of the Universe. The frequency of the oscillation is unfortunately too high to detect the oscillation as a "gravitational wave" by the present and future gravitational wave detectors. However, if the graviton mass can be sufficently light due to, for example, the Z 2 symmetry, the coherent oscillation will be detectable.
Finally, we comment on an interesting remaining question: Is the almost homogeneous configuration of ϕ µν the Bose-Einstein condensate of the massive graviton? As is well-known, a coherent massive scalar field, for example axion, is a viable dark matter candidate [60][61][62][63][64][65][66][67][68][69][70][71][72][73][74]. This coherent scalar field is interpreted as the Bose-Einstein condensate. Our result would be a generalization of the Bose-Einstein condensate of the massive scalar field to that of the massive tensor field. However, the present argument is completely classical and we have not discussed any quantum aspect of the massive graviton. Therefore, it would be interesting to study a connection to the quantum theory of the gravitation, but this is beyond the scope of the present paper. In this section, we summarize the definitions of the graviton energy-momentum tensors for generic cases. For completeness, we introduce both the g-matter fields ψ g and the f -matter fields ψ f and do not assume (3.32). The matter action is given by and we denote the energy-momentum tensors of the gmatter field and the f -matter field as T µν and T µν , respectively.
To define the energy-momentum tensor of gravitons, we shall decompose the metric into the "background" and the "perturbations". However, in general, the deomposition into the perturbations and the background may not be well-defined because the perturbations and the background interact with each others and then the equations may not be separable when the backreaction from the perturbations to the background is included. To decompose the Einstein equations, we assume the perturbations contain only high-frequency modes whereas the background consists of only the low-frequency modes. In this case, we obtain two independent eqautions for the low-frequency background and the high-frequency perturbations via a low-frequency projection · · · low and a high-frequency projection · · · high , respectively. Therefore, we assume the metrics are expressed by the lowfrequency backgrounds with the high-frequency perturbations: (A2) with |δg The high-frequency mode and the low-frequency mode are defined by which is same for f µν . It is worth noting that we only assume the perturbations are high-frequency modes but do not assume the perturbations are high-momentum modes. In GR, the situations with the high-frequency waves and the situations with the high-momentum waves are equivalent since the graviton is massless. However, these situations are not equivalent in bigravity due to the existence of the massive graviton. As already mentioned, the definitions of the massless mode and the massive mode are ambiguous in general. They can be defined when the curvature scale of the metrics are smaller than the graviton mass squared In this case, although the low-frequency modes of the massive gravitons can be excited by some source including matter as well as the backreactions from highfrequency modes, the amplitudes of the low-frequency massive modes must be tiny (Indeed, in Appendix B, we will see the amplitude of the low-frequency massive mode is suppressed by m −2 for the homogeneous matter distributions.). Hence, we may assume the background is approximated by the homothetic solution. Then, the spacetime are expressed as spacetimes are expressed as where (0) g µν and M µν are the low-frequency massless mode and the low-frequency massive mode while h µν and ϕ µν are the high-frequency massless mode and the low-frequency massive mode, respectively (see Table I).
The assumption that the homothetic background is good approximation means |M µν | ≪ | (0) g µν |. In this case, M µν , h µν and ϕ µν can be treated as the tensors with respect to The Ricci tensor for the g-spacetime is expanded as where ∆g µν = M µν + δg (high) µν . The high-frequency modes of g µν and f µν are written in terms of the mass eigenstates as The functionals δ (1) R µν [χ] and δ (2) R µν [χ] are defined by for a tensor χ µν . The linear quantities in M µν and δg (high) µν are decomposed whereas the quadratic quantity δ (2) R µν [∆g] have the cross terms between M µν and δg (high) µν . The quadratic quantity is given by After taking the high/low-frequency projections, we obtain The other quantities are expanded in the similar way and then the equations are decomposed into ones for the low-frequency modes and for the high-frequency modes, respectively. Up to the linear order the high-frequency mode equations are decomposed into the massless one and the massive one where we define (A20) Note that Substituting this into the trace of (A17), we find a constraint equation on the trace ϕ α α From the low-frequency mode equations, we obtain the Einstein equation for the homothetic background where (0) G µν is the Einstein tensor for is interpreted as the total energy-momentum tensor including the gravitons as well as the matters. The effective matter energy-momentum tensors for matters and the graviton energy-momentum tensors are defined by and (3.22).
Finally, we derive the equation for the massive mode M µν which is given by where and ∆ (2) Note that the source term ∆S µν is given by the difference between two matter energy-momentum tensors whereas the source term for the massless mode is given by the sum of energy-momentum tensors. The massive mode can give an anti-gravity since the positiveness of the source is not guaranteed even if all energy of the sources are positive definite.
Choosing the gauge N g = 1, we find following equations: the Friedmann equations the constraint equation and the equations for anisotropies where we define and Eqs. (B12) and (B13) yield Hence, when the anisotropic stress is ignoredπ = 0, the sum of the shears σ g and σ f decreases as a −3 g which is the same as the standard decaying law of the shear in GR. On the other hand, the difference between them does not decreases as a −3 g due to the "potential" U β . Instead, σ g − σ f (and also β g − β f ) decrease as a −3/2 g as shown in Fig. 2 and then acts as the "dark matter" component of the universe (see also [22]). at the late time where we assume Λ g | ξ0=1 = 0. Therefore, even if the matter component is not introduced, the same behavior as the dust dominant universe can be obtained.
We then assume β ≪ 1 and consider the regime m 2 ≫ H 2 g , H 2 f . Note that the smallness of β do not suggest that contributions of those to the Friedmann equations are also small because the quantity m 2 β 2 can be comparable to H 2 g and H 2 f . We introduce a small quantity ǫ := H g /m and set β = O(ǫ). Eq. (B11) gives The evolutions of N f , ξ and the rescaled Hubble expansion rate (which is scale free due to the scale factor). We set the same parameters as Fig. 2. where we have used σ g , σ f H g , H f . Since the spacetime can be approximated by the homothetic spacetime in the regime m 2 ≫ H 2 g , H 2 g , the ratio ξ can be expanded around ξ 0 . Supposing the normal branch such that the Friedmann equations yield The deviation of the lapse function is then given by where we notice that, althoughξ and H g have quantities of order O(ǫ), they are canceled and then N f − ξ 0 is of order O(ǫ 2 ). Using m 2 eff ≫ Λ g and including only leading order contributions, the Friedmann equation is expressed by The variablesφ andh are the normalized massive mode and the normalized massless mode of the anisotropies, or, following the notion of the main text, they can be interpreted as the "massive graviton condensate" and the "massless graviton condensate", respectively, which obeÿ It is worth noting that although we have assumed the inequalities β ≪ 1 and m 2 ≫ H 2 g , H 2 g we have not used the high-frequency and the low-frequency projections in the present calculations.
When we ignore the anisotropic stress, the (averaged) energy densities ofh andφ decrease as The effect of the homogeneous mode of the massless graviton can be ignored in time. The averaged differences ξ − ξ 0 T and N f − ξ 0 T correspond to the low-frequency massive mode M µν which are given by (B31) As we expected, M µν is suppressed by m −2 which is just a sub-subleading contribution.
Finally, we discuss the production of the anisotropies. We assume the matter field is composed of radiation and the coherent magnetic field: whereB is a constant and the strength of the magnetic field is given byBe 2βg /a 2 g . A typical behavior ofφ is shown in Fig. 1 8 . The dominant production occurs just after H g = m eff in which Eq. (B29) can be barely used. The produced amplitude is then estimated as In the main text, since we have assumed M G = M pl and used the normalization m eff = m, we obtain (6.1).

Appendix C: General perturbations
In this section, we summarize the calculations for the general perturbations around the homogeneous solutions (4.1) and (4.2). We use the specific choice of the coupling constants κ 2 g = κ 2 f and (3.36) to simplify the calculations.

Harmonics expansion
The homogeneity of the background solution leads to that all variables can be transformed into the momentum space. Furthermore, due to the rotational symmetry in the y-z space, the perturbations are decomposed into the even parity perturbations and the odd parity perturbations and there is no interaction between the even and the odd parity perturbations at linear order. The even and odd parity perturbations are parametrized as where and with Y S = e ikx . For our study, it is useful to define the harmonics associated with the three dimensional Euclidean space. We define a three dimensional vector and a tensor as and which satisfy (5.17) and (5.18). Using the three dimensional harmonics, any three dimensional vector (or tensor) is uniquely decomposed into the scalar and the vector (and the tensor) quantities. For example, the even parity perturbations and the odd parity perturbations of the graviton energy-momentum tensor are given by (5.25) and (5.26), respectively.

Perturbations for massless mode
The low-frequency mode (0) g µν has the gauge symmetry by which we can obtain (5.23) and (5.24).
The general form of the even parity metric perturbations is given by with and where we note and Under the gauge transformation with the perturbations transform as Φ g → Φ g +ξ 0 , and then where k 2 ξ L = k 2 x ξ 1 + k 2 ξ 2 .
Therefore one can fix the gauge so that in which B i and H ij are given by where B V and H T are functions of time.
The odd parity metric perturbations are with and where we notice Under the gauge transformation with ξ µ (odd) = (0, 0, ξ odd Y a ) , the metric perturbations transform as and then One can choose the gauge in which B i g and H ij g corresponds to the vector perturbation and the tensor perturbation, respectively:

Adiabatic expansion
In addition to the linearization of the inhomogeneities, we shall take the adiabatic expansion in terms of m −1 . To verify the linearization the amplitudes of the inhomogeneities are of O(m n ) with n ≤ 0.
The equation (5.10) yields that, up to the subleading order, the massive graviton can be expressed as δϕ µν = δϕ 1µν cos[mt] + δϕ 2µν sin[mt] , where δϕ 1µν and δϕ 2µν are slowly varying function in time 9 . We use the suffixes 1 and 2 to represent the slowly varying functions in front of cos[mt] and sin[mt] (see (5.32)). The orders of δϕ 1µν and δϕ 2µν are determined to be consistent with the equations.

Odd parity
We first study the odd parity perturbations in which (5.12) is trivially satisfied. We discuss the large scales (5.4) and the small scales (5.7) in order.

a. Large scales
The consistency of the equations yield and other variables are of order O(m 0 ). Form the equations of motion (5.10) and (5.11) we obtain the constraint equations and the dynamical equationṡ while the dynamical equation arė The solutions are when a ∝ t 2/3 (the dust dominant universe).
We then obtain the velocity Since D 1,2 and E 1,2 are independent, we can choose D 1,2 and D 1,2 + 2E 1,2 as new variables which obey the same equation as (C48) and (C49). The vector quantities v (odd) and π (odd) V evolve independently from the tensor quantity π (odd) T . As a result, the vector and the tensor modes are decoupled.

Even parity a. Large scales
In the large scales (5.4), we obtain We cannot find other equations within our accuracy. We note that the amplitudes of D 1,2 and E 1,2 always decreases since the gravitational potential Φ does not affect the dynamics of D 1,2 and E 1,2 . The Jeans instability does not lead to the growth of D 1,2 and E 1,2 .
After taking the oscillation average, we find that the pressures are zero. The energy density and v S are given by (5.36) and (5.37), respectively. The vector mode of the velocity is The massive graviton condensate is a form of a pressureless perfect fluid. Clearly from the definitions of δρ G , v S and v V , the vector mode is decoupled from other modes. As already mentioned in the main text, the evolution of δϕ 1 is not determined. However, the dynamics of δρ G is determined by the conservation law of the averaged graviton energy-momentum tensor.

b. Small scales
Finally, we study the even parity perturbations in the small scales (5.7). The consistency leads to