ournal of C osmology and A stroparticle hysics J Dipolar dark matter with massive bigravity

. Massive gravity theories have been developed as viable IR modiﬁcations of gravity motivated by dark energy and the problem of the cosmological constant. On the other hand, modiﬁed gravity and modiﬁed dark matter theories were developed with the aim of solving the problems of standard cold dark matter at galactic scales. Here we propose to adapt the framework of ghost-free massive bigravity theories to reformulate the problem of dark matter at galactic scales. We investigate a promising alternative to dark matter called dipolar dark matter (DDM) in which two diﬀerent species of dark matter are separately coupled to the two metrics of bigravity and are linked together by an internal vector ﬁeld. We show that this model successfully reproduces the phenomenology of dark matter at galactic scales (i.e. MOND) as a result of a mechanism of gravitational polarisation. The model is safe in the gravitational sector, but because of the particular couplings of the matter ﬁelds and vector ﬁeld to the metrics, a ghost in the decoupling limit is present in the dark matter sector. However, it might be possible to push the mass of the ghost beyond the strong coupling scale by an appropriate choice of the parameters of the model. Crucial questions to address in future work are the exact mass of the ghost, and the cosmological implications of the model.


Introduction
In the last century, cosmology has progressively developed from a philosophical to an empirical scientific discipline, witnessing high precision cosmological observations, which culminated with the standard model of cosmology, the Λ-CDM model [1]. The standard model is based on General Relativity (GR) and is in great agreement with the abundance of the light elements in the big bang nucleosynthesis, the anisotropies of the cosmic microwave background (CMB), the baryon acoustic oscillations, lensing and the observed large scale structures. It notably relies on the presence of dark energy in form of a cosmological constant Λ, giving rise to the accelerated expansion of the universe.

Massive gravity context
In the prevailing view, the cosmological constant corresponds to the constant vacuum energy density which receives large quantum corrections. Unfortunately, so far there is no successful mechanism to explain the observed unnatural tiny value of the cosmological constant, see e.g. [2,3]. This difficulty has sparked a whole industry studying modifications of gravity in the infra-red (IR) invoking new dynamical degrees of freedom.
One tempting route is massive gravity, motivated by the possibility that the graviton has a mass. It was a challenge over forty years to construct a covariant non-linear theory for massive gravity. The foregathered remarkable amount of effort finally rose to the challenge [4][5][6][7][8][9][10][11], to construct the potential interactions in a way that got rid of the Boulware-Deser (BD) ghost [12]. This theory has been further extended to more general models by adding additional degrees of freedom. A very important extension is massive bigravity theory [13].
Once the interactions in the gravitational sector were guaranteed to be ghost-free, a natural follow up question was how to couple this theory to the matter fields without spoiling the ghost-freedom. First attempts were already discussed in the original paper of bigravity [13]. If one couples the matter fields to both metrics simultaneously, this reintroduces the BD ghost [14,15]. Furthermore, the one-loop quantum corrections detune the special potential interactions at an unacceptable low scale and hence this way of coupling is not a consistent -1 -

JCAP12(2015)026
one [15]. The safer way is to couple the matter fields to just one metric. In this way, the quantum corrections give rise to contributions in form of cosmological constants. One could also try to couple the matter field to the massless mode, which is unfortunately also not ghost free [16]. Another possible way of coupling the matter fields to both metrics is through a composite effective metric constructed out of both metrics [15,[17][18][19], which is unique in the sense that it is the only non-minimal matter coupling that maintains ghost-freedom in the decoupling limit [20][21][22]. Furthermore, the quantum corrections are guaranteed to maintain the nice potential structure. Other important consequence of this new matter coupling is the fact that it helps to evade the no-go result [23] for the flat Friedmann-Lemaître-Robertson-Walker (FLRW) background. A detailed perturbed ADM analysis revealed the existence of a BD ghost originating from an operator involving spatial derivatives [15,17]. Therefore, the ghost will probably reappear for highly anisotropic solutions. Since the ghost remains absent up to the strong coupling scale, the matter coupling can be considered in an effective field theory sense at the very least till the strong coupling scale. The precise cut-off of the theory or mass of the ghost has still to be established.
The absence of the BD ghost is not only important at the classical level, but also at the quantum level. For the classical ghost-freedom the relative tuning of the potential interactions is the key point. Therefore one has to ensure that the quantum corrections do not detune the potential interactions. Concerning the decoupling limit, it is easy to show that the theory receives no quantum corrections via the non-renormalization theorem due to the antisymmetric structure of the interactions [24] (in fact the same antisymmetric structure of the Galileon interactions protect them from quantum corrections [25][26][27]). Beyond the decoupling limit, the quantum corrections of the matter loops maintain the potential interactions provided that the above criteria are fulfilled. Concerning the graviton loops, they do destroy the relative tuning of the potential interactions. Nevertheless, this detuning is harmless since the mass of the corresponding BD ghost is never below the cut-off scale of the theory [28]. The bimetric version of the theory shares the same property [29].
Massive gravity/bigravity theory has a rich phenomenology. The decoupling limit admits stable self-accelerating solutions [30], where the helicity-0 degree of freedom of the massive graviton plays the role of a condensate whose energy density sources self-acceleration. Unfortunately these solutions suffer from strong coupling issues due to the vanishing kinetic term of the vector modes [30,31]. With the original massive gravity and the restriction of flat reference metric, one faces the no-go result, namely that there are no flat FLRW solution [23]. One can construct self-accelerating open FLRW solutions [32], which however have three instantaneous modes [33] and suffer from a nonlinear ghost instability [34]. The attempt to promote the reference metric to de Sitter [35,36] also failed due to the presence of the Higuchi ghost [35]. In order to avoid these difficulties, one either gives up on the FLRW symmetries [23], or invokes new additional degrees of freedom [13,37,38].
Thanks to the freedom gained in the inclusion of the second kinetic term in the bimetric extension [13], there now exists many elaborate works concerning the cosmology of the bigravity theory, see e.g. [39][40][41]. In the case of minimally coupled matter fields and small graviton mass, the theory admits several interesting branches of solutions. Unfortunately, among them the self-accelerating branch is unstable due to the presence of three instantaneous modes, and a second branch of solutions admits an early time gradient instability [42]. Nevertheless, there exists attempts to overcome the gradient instability either by viable though finely tuned solutions in the case of a strongly interacting bimetric theory with m ≫ H 0 [43,44], or by demanding that M g ≫ M f as was proposed in [49]. See also other recent works concerning the phenomenology of bimetric gravity [45][46][47][48]50].

Dark matter context
On a quite different but equally fascinating topic is the phenomenology of dark matter at galactic scales. The evidence for dark matter is through the measurement of the rotation curves of spiral galaxies which turn out to be approximately flat, contrary to the Newtonian prediction based on ordinary baryonic matter [51,52]. The standard explanation is that the disk of galaxies is embedded into the quasi-spherical potential of a huge halo of dark matter, and that this dark matter is the same as the cold dark matter (CDM) which is evidenced at large cosmological scales notably with the fit of Λ-CDM to the CMB anisotropies [53,54]. Unfortunately this explanation faces severe challenges when compared to observations at galactic scales [55,56]. There are predictions of the Λ-CDM model that are not observed, like the phase-space correlation of galaxy satellites and the generic formation of dark matter cusps in the central regions of galaxies. Even worse, there are observations which are not predicted by Λ-CDM, such as the tight correlation between the mass discrepancy (luminous vs. dynamical mass) which measures the presence of dark matter and the involved scale of acceleration, the famous baryonic Tully-Fisher (BTF) relation for spiral galaxies [57][58][59], and its equivalent for elliptical galaxies, the Faber-Jackson relation [60].
Instead of additional, non-visible mass, Milgrom [61][62][63] proposed an amendment to the Newtonian laws of motion in order to account for the phenomenology of dark matter in galaxies, dubbed MOND for MOdified Newtonian Dynamics. According to the MOND hypothesis the change has a relevant influence on the motion only for very small accelerations, as they occur at astronomical scales, below the critical value a 0 ≃ 1.2 × 10 −10 m/s 2 measured for instance from the BTF relation. 1 A more elaborate version of MOND is the Bekenstein-Milgrom [64] modification of the Poisson equation of Newtonian gravity. All the challenges met by Λ-CDM at galactic scales are then solved -sometimes with incredible successby the MOND formula [55,56]. Unfortunately, MOND faces difficulties in explaining the dark matter distribution at the larger scales of galaxy clusters [65][66][67][68][69]. It is also severely constrained by observations in the solar system [70,71].
Reconciling Λ-CDM at cosmological scales and MOND at galactic scales into a single relativistic theory is a great challenge. There has been extensive works on this. One approach is to modify gravity by invoking new gravitational degrees of freedom without the presence of dark matter [72][73][74][75][76][77][78][79][80][81]. Primary examples are the tensor-vector-scalar (TeVeS) theory [73,74] and generalized Einstein-AEther theories [75,76]. Another approach is a new form of dark matterà la MOND, called dipolar dark matter (DDM). It is based on the dielectric analogy of MOND [82] -a remarkable property of MOND with possible far-reaching implications. Indeed the MOND equation represents exactly the gravitational analogue (in the non-relativistic limit) of the Gauss equation of electrostatics modified by polarisation effects in non-linear dielectric media. Some early realizations of this analogy were proposed in [83][84][85] and shown to also reproduce the cosmological model Λ-CDM. The best way to interpret this property is by a mechanism of gravitational polarisation, involving two different species of particles coupled to different gravitational potentials. This was the motivation for pushing forward a more sophisticated model in the context of a bimetric extension of GR [86,87]. In this model the two species of dark matter particles that couple to the two metrics are linked through an internal vector field. The model [87] is able to recover MOND at galactic scales and the cosmological model Λ-CDM at large scales. Furthermore the post-Newtonian parameters (PPN) in the solar system are the same as those of GR. Unfortunately, 1 A striking observation is the numerical coincidence that a0 ∼ √ Λ [55,56]. even if this model delivers a very interesting phenomenology, it is plagued by ghost instabilities in the gravitational sector already at the linear order of perturbations and within the decoupling limit [88]. A more promising model was then proposed in [88], with the hope to cure these instabilities in the gravitational sector.
In the present work we will study in detail the theoretical and phenomenological consistency of this recently proposed model for dipolar dark matter in the context of bimetric gravity [88]. We will first work out the covariant field equations in section 2 and analyze the linear field equations in section 3. As next, we will investigate in section 4 the mechanism of gravitational polarisation in the non relativistic approximation, and show how it successfully recovers the MOND phenomenology on galactic scales. Finally we will pay attention in section 5 to the matter sector and investigate the interactions in the decoupling limit. The outcome of section 5 is that a ghost instability is still present in the (dark) matter sector of the model. The paper ends with some short conclusions in section 6.

The covariant theory
A new relativistic model for dipolar dark matter has been recently proposed in [88]. Basically this model is defined by combining the specific dark matter sector of a previous model [87] with gravity in the form of ghost-free bimetric extensions of GR [13]. The dark matter in this model consists of two types of particles respectively coupled to the two dynamical metrics g µν and f µν of bigravity. In addition, these two sectors are linked together via an additional internal field in the form of an Abelian U(1) vector field A µ . This vector field is coupled to the effective composite metric g eff of bigravity which is built out of the two metrics g and f , and which has been shown to be allowed in the matter Lagrangian in the effective field theory sense at least up to the strong coupling scale [15,17]. The total matter-plus-gravity Lagrangian reads 2 Here R g and R f are the Ricci scalars of the two metrics, and the scalar energy densities of the ordinary matter modelled simply by pressureless baryons, 3 and the two species of pressureless dark matter particles are denoted by ρ bar , ρ g and ρ f respectively. In addition j µ g , j µ f stand for the mass currents of the dark matter, defined by where J µ g = r g ρ g u µ g and J µ f = r f ρ f u µ f are the corresponding conserved dark matter currents associated with the respective metrics g and f , thus satisfying ∇ g µ J µ g = 0 and ∇ f µ J µ f = 0. Here ρ g and ρ f denote the scalar densities and u µ g , u µ f are the four velocities normalized to g µν u µ g u ν g = −1 and f µν u µ f u ν f = −1. We also introduced two constants r g and r f specifying the ratios between the charge and the inertial mass of the dark matter particles.

JCAP12(2015)026
The metric independent matter degrees of freedom are the coordinate densities ρ * g = √ −gρ g u 0 They are conserved in the ordinary sense, ∂ µ J * µ g = 0 and ∂ µ J * µ f = 0. They relate to the classical currents J µ g , J µ f or j µ g , j µ f by Notice that the coupling term √ −g eff A µ (j µ g − j µ f ) in the action (2.1) is actually independent of any metric, neither g nor f nor g eff . We shall study in detail in section 5 the implications of the term √ −g eff A µ (j µ g − j µ f ) for the decoupling limit of the theory. The vector field A µ is generated by the dark matter currents and plays the role of a "graviphoton". The presence of this internal field is necessary to stabilize the dipolar medium and will yield the wanted mechanism of gravitational polarisation, as we shall see in section 4. It has a non-canonical kinetic term W(X ), where W is a function which is left unspecified at this stage, and Here λ is a constant and the field strength is constructed with the effective composite metric g eff µν given by The model (2.1) is defined by several constant parameters, the coupling constants M g , M f and M eff , the mass of the graviton m, the charge to mass ratios r g and r f of dark matter, the constant λ associated with the vector field, and the arbitrary constants α and β entering the effective metric (2.4). Parts of these constants, as well as the precise form of the function W(X ), will be determined in section 4 when we demand that the model reproduces the required physics of dark matter at galactic scales. In particular λ will be related to the MOND acceleration scale a 0 , and the model will finally depend on a 0 and the mass of the graviton m. Additionally, there will be still some remaining freedom in the parameters α and β, and in the coupling constants M g and M f .
In the model (2.1) the ghost-free potential interactions between the two metrics g and f take the particular form of the square root of the determinant of the effective metric [15,17,19] g eff where α and β are arbitrary constants, and we have defined where the square root matrix is defined by Interchanging the two metrics g and f does not change the form of the metric (2.4)-(2.5) except for a redefinition of the parameters α and β. Notice that G eff µν can be proved to be automatically symmetric [16,89], and in [88] it was shown, that it corresponds to the composite metric considered in the previous model [87]. The square root of the determinant of g eff µν is given in either forms respectively associated with the g or f metrics, by

JCAP12(2015)026
More explicitly, introducing the elementary symmetric polynomials e n (X) and e n (Y ) of the square root matrices X or Y , 4 e 0 (X) = 1 , we can write, still in symmetric guise, In this form we see that (2.6) corresponds to the right form of the acceptable potential interactions between the metrics g and f . We can first vary the Lagrangian (2.1) with respect to the two metrics g and f , which yields the following two covariant Einstein field equations [90] where G µν g and G µν f are the Einstein tensors for the g and f metrics. The tensors U µν (n) and V µν (n) are defined from the following sums of matrices [8] 5 by raising or lowering indices with their respective metrics, thus U µν and V µν (n) are indeed automatically symmetric. 4 As usual we denote the traces of rank-2 tensors as X µ µ = [X], X µ ν X ν µ = [X 2 ], and so on. Recall in particular that √ −g en(X) = √ −f e4−n(Y ), with e4(X) = det(X) and e4(Y ) = det(Y ). 5 Here we adopt a slight change of notation with respect to [8]. Our notation is related to the one of [8] by

JCAP12(2015)026
In the right sides of (2.9) the stress-energy tensors T µν bar and T µν g are defined with respect to the metric g, T µν f is defined with respect to f and T µν g eff with respect to g eff . Thus, for pressureless baryonic and dark matter fluids we have T µν bar = ρ bar u µ bar u ν bar , T µν g = ρ g u µ g u ν g and T µν f = ρ f u µ f u ν f in terms of the corresponding scalar densities and normalized four velocities, while T µν g eff represents in fact the stress-energy tensor of the vector field A µ , i.e.
where W X is the derivative of W with respect to its argument X defined by (2.3).
The equation of motion for the baryons follows a geodesic law a bar µ ≡ u ν bar ∇ g ν u bar µ = 0, whereas the equations of motion for the dark matter particles are non-geodesic, where the accelerations a g µ ≡ u ν g ∇ g ν u g µ and a f µ ≡ u ν f ∇ f ν u f µ are defined with their respective metrics. On the other hand, the equation for the vector field yields with ∇ g eff ν denoting the covariant derivative associated with g eff . This equation is obviously compatible with the conservation of the currents, since by (2.2) we have also ∇ g eff µ j µ g = 0 and ∇ g eff µ j µ f = 0. The equation (2.13) can be also written as Combining the equations of motion (2.12) with (2.14) we obtain the conservation law Alternatively, such global conservation law can be obtained from the scalarity of the total matter action under general diffeomorphisms.

Linear perturbations
In this section we will be interested in computing the linear perturbations of our model, which will be extensively used in the following section for the study of the polarisation mechanism. Following the ideas of [87] we will assume that the two fluids of dark matter particles are slightly displaced from the equilibrium configuration by displacement vectors y µ g and y µ f . This makes the dark matter medium to act as an analogue of a relativistic plasma in electromagnetism. The dark matter currents will be slightly displaced from the equilibrium current j µ 0 = ρ 0 u µ 0 as well, with j µ 0 satisfying ∇ g eff µ j µ 0 = 0. Furthermore, we will perturb the two metrics around a general background in the following way 6

1a)
f µν = (f µν + ℓ µν ) 2 , (3.1b) 6 By which we mean, for instance, gµν = (ḡµρ + hµρ)ḡ ρσ (ḡσν + hσν ). with main interest in background solutions whereḡ µν =f µν . In the latter case the potential interactions take the form where k µν = αh µν + βℓ µν , and traces are defined with the common background metricḡ =f . Working at first order in the displacement of the dark matter particles and assuming that their gradients are of the same order as the metric perturbations, the dark matter currents can then be expanded as [87] with the projector operator perpendicular to the four velocity of the equilibrium u µ 0 given as ⊥ µν = g µν eff +u µ 0 u ν 0 , and the remainder O(2) indicating that the expansion is at first order. The desired plasma-like solution is achieved after plugging our ansatz (3.3) into the equations of the vector field (2.13), which results in where we have defined the projected relative displacement as ξ µ ≡⊥ µ ν (y ν g − y ν f ). Note that ξ µ is necessarily a space-like vector. In section 4 we shall see that the spatial components ξ i of this vector, which can be called a dipole moment, define in the non-relativistic limit the polarisation field of the dark matter medium as P i = ρ * 0 ξ i , where ρ * 0 is the coordinate density associated with ρ 0 . But for the moment we only need to notice that ξ µ is a first order quantity, therefore the field strength F µν is itself a first order quantity, and hence the stress-energy tensor (2.11) with respect to g eff is already of second order, 7 For this reason, in the case of our desired plasma-like solution the contributions of the selfinteractions of the internal vector field will be at least third order in perturbation at the level of the action. For more detail see the comprehensive derivations in [87]. Let us first expand the Lagrangian (2.1) to first order in perturbation around a common backgroundḡ =f . We assume that there is no matter in the background, so the matter interactions do not contribute at that order. Using (3.2) we obtain (modulo a total divergence) where Gḡ µν is the Einstein tensor for the background metricḡ, and Lḡ is the background value of the Lagrangian. In order for the backgroundḡ µν in our interest to be a solution of the theory, we have to impose that linear perturbations vanish. We find

JCAP12(2015)026
This is only possible if With this choice we guarantee that we are expanding around the backgroundḡ µν =f µν that is a solution to the background equations of motion. As next we shall compute the action to second order in perturbations. Let us emphasize again that thanks to the plasma-like ansatz (3.3) the self-interaction of the internal vector field does not contribute to this order, see (3.5). Our quadratic Lagrangian above the background g µν simply becomes whereĒ ρσ µν is the Lichnerowicz operator on the backgroundḡ given in the general case by with Cḡ µρσν denoting the Weyl curvature of the background metric. Of course this expression is to be simplified using the background equations (3.7). Also bear in mind that αM 2 f = βM 2 g for having a common backgroundḡ =f . The field equations obtained by varying with respect to h µν and ℓ µν yield where we recall that k µν = αh µν + βℓ µν and k = [k] =ḡ µν k µν . Of course those equations can also be recovered directly from the general Einstein field equations (2.9).
In the present paper we shall mostly make use of that combination of the two Einstein field equations (3.11) which corresponds to the propagation of a massless spin-2 field. We subtract the two equations from each other so as to cancel the mass term, resulting in However we shall also use the contribution of the mass term in the linearised Bianchi identities associated with (3.11). Thus, last but not least we can act on (3.11) with ∇ḡ ν , yielding In deriving this relation we use the result that, for instance, ∇ḡ ν T µν g = ρ g a µ g = O(2) which comes from the equations of motion (2.12) and the fact that there is no matter in the background. We shall see in the next section how important is this relation for the present model to recover the looked-for phenomenology of dark matter at galactic scales.
-9 -We next consider the Newtonian or non-relativistic (NR) limit of our model. In the previous section we derived the field equations at linear order around a vacuum background metric g, which necessarily obeys the equations (3.7). Thus, the background is a de Sitter solution, with cosmological constant given by the graviton's mass, Λ ∼ m 2 . When computing the Newtonian limit applied to describe the physics of the local universe, e.g. the solar system or a galaxy at low redshift, we shall neglect the cosmological constant and shall approximate the de Sitter metricḡ by a flat Minkowski background η.
The best way to implement (and justify) this approximation is to perform an expansion on small scales, say r → 0. When we solve for the gravitational field in the solar system or a galaxy embedded in the de Sitter background, we get terms like 1−2M r −1 −Λr 2 /3 (de Sitter-Schwarzschild solution), where M is the mass of the Sun or the galaxy. On small scales the term 2M r −1 will dominate over Λr 2 /3, so we can ignore the influence of the cosmological constant if the size of the system is ≪ (6M/Λ) 1/3 . In the present case we shall require that the size of our system includes the very weak field region far from the system where the MOND formula applies. In that case the relevant scale is r 0 = M/a 0 , where a 0 ≃ 1.2 × 10 −10 m/s 2 is the MOND acceleration. So our approximation will make sense provided that r 0 ≪ (6M/Λ) 1/3 . For a galaxy with mass M ∼ 10 11 M ⊙ the MOND transition radius is r 0 ∼ 10 kpc. With the value of the cosmological constant Λ ∼ (3 Gpc) −2 we find (6M/Λ) 1/3 ∼ 700 kpc, so we are legitimate to neglect the influence of the cosmological constant. For the solar system, r 0 ∼ 0.04 pc while (6M ⊙ /Λ) 1/3 ∼ 150 pc and the approximation is even better.
Actually we shall make the approximationḡ ≃ η only on the particular massless spin-2 combination of the two metrics given by (3.12), namely where the Lichnerowicz operator (3.10) is now the flat one, and also in the constraint equation (3.13). We shall not need to consider the other, independent (massive) combination of the two metrics. In the NR limit when c → ∞ we parametrize the two metric perturbations h µν and ℓ µν by single Newtonian potentials U g and U f such that We shall especially pay attention to the potential U g since it represents the ordinary Newtonian potential felt by ordinary baryonic matter. Similarly we parametrize the internal vector field A µ by a single Coulomb scalar potential φ such that

JCAP12(2015)026
It is now straightforward to show that (4.1) reduces in the NR limit to a single scalar equation for a combination of the Newtonian potentials U g and U f , 8 where ∆ = ∇ 2 is the ordinary Laplacian, and ρ * bar , ρ * g and ρ * f are the ordinary Newtonian (coordinate) densities, satisfying usual continuity equations such as ∂ t ρ * g + ∇ · (ρ * g v g ) = 0. Notably, all non-linear corrections O(2) in (4.1) are negligible. In principle, we should a priori allow some parametrized post-Newtonian (PPN) coefficients in the spatial components ij of the metrics (4.2), say γ g and γ f . However, when plugging (4.2) into (4.1) we determine from the ij components of (4.1) that in fact γ g = 1 and γ f = 1. Similarly we find that the equations (2.13) governing the internal vector field reduce to a single Coulomb type equation, Notice that the constants α and β in the effective metric g eff to which is coupled the internal field intervene in the NR limit only in the expression Next, the equations of the baryons which are geodesic in the metric g reduce to dv bar dt = ∇U g , (4.7) while the equations of motion (2.12) of the dark matter particles become dv g dt = ∇U g + r g ∇φ , At this stage the important point is to recall that the two sectors associated with the metrics g and f do not evolve independently but are linked together by the constraint equation (3.13). We now show that this constraint provides a mechanism of gravitational polarisation in the NR limit and yields the MOND equation for the potential U g felt by baryons. In flat space-time this constraint is where k µν = αh µν + βℓ µν and k = η ρσ k ρσ . Plugging (4.2) into (4.9) we readily obtain a constraint on the gravitational forces felt by the DM particles, namely ∇ αU g + βU f = 0 .

JCAP12(2015)026
Note that the constraint (4.10) comes from a combination between the 00 and ij components of the metrics (4.2). Using (4.10) we express the equations of motion (4.8) solely in terms of the potential U g ruling the baryons, As we see, the gravitational to inertial mass ratio m g /m i of the f particles when measured with respect to the g metric is m g /m i = −α/β. We look for explicit solutions of (4.11) in the form of plasma-like oscillations around some equilibrium solution. We immediately see from (4.11) the possibility of an equilibrium (i.e. for which dv g /dt = dv f /dt = 0) when we have the following relation between constants, The equilibrium holds when ∇U g + r g ∇φ = 0, i.e. when the Coulomb force annihilates the gravitational force. To describe in the proper way the two DM fluids near or at equilibrium we use the NR limits of the relations (3.3), with an appropriate choice of the equilibrium configuration. From (4.11) we note that (α + β)v 0 = αv g + βv f is constant, hence we define the equilibrium in such a way that the two displacement vectors y g and y f with respect to that equilibrium obey αy g + βy f = 0, and in particular we choose the equilibrium velocity to be v 0 = 0. We shall now define the relative displacement or dipole moment vector by ξ = y g − y f . 9 Hence (3.3) imply where ρ * 0 is the common density of the two DM fluids at equilibrium in the absence of external perturbations (far from any external mass), and P = ρ * 0 ξ is the polarisation. The velocity fields of the two fluids are where dξ/dt = ∂ t ξ + v 0 · ∇ξ denotes the convective derivative, which reduces here to the ordinary derivative since v 0 = 0. By inserting (4.13) into (4.5) we can solve for the polarisation resulting in Furthermore, by combining (4.14) and (4.11) and making use of the previous solution for the polarisation (4.15), we readily arrive at a simple harmonic oscillator describing plasma like oscillations around equilibrium, namely (4.16) 9 Here yg and y f denote the spatial components of the displacement four vectors ⊥ µ ν y ν g and ⊥ µ ν y ν f considered in section 3, see

JCAP12(2015)026
In this case the plasma frequency is given by (4.17) Next, consider the equation for the Newtonian potential U g in the ordinary sector. Combining (4.10) with (4.4) we readily obtain To recover the usual Newtonian limit we must impose, in geometrical units, This condition being satisfied, the right-hand side of (4.18) can be rewritten with the help of (4.13) and the relation (4.12). We obtain an ordinary Poisson equation but modified by polarisation effects, We still have to show that the polarisation will be aligned with the gravitational field ∇U g , i.e. we grasp a mechanism of gravitational polarisation. This is a consequence of the constitutive relation (4.15) taken at the equilibrium point, neglecting plasma like oscillations. At this point we have ∇U g + r g ∇φ = 0. Thus, Finally (4.20) with (4.21) takes exactly the form of the Bekenstein-Milgrom equation [64].
To recover the correct deep MOND regime we must impose that when X → 0 where a 0 is the MOND acceleration scale. This is easily achieved with the choice together with fixing the constants M eff and λ to the values We thus recovered gravitational polarisation and the MOND equation (in agreement with the dielectric analogy of MOND [82]) as a natural consequence of bimetric gravity since it is made possible by the constraint equation (4.9) linking together the two metrics of bigravity. 10 The polarisation mechanism is possible only if we can annihilate the gravitational JCAP12(2015)026 force by some internal force, here chosen to be a vector field, and therefore assume a coupling between the two species of DM particles living in the g and f sectors.
Let us recapitulate the various constraints we have found on the parameters in the original action (2.1). These are given by (4.12) for having a polarisation process, (4.19) for recovering the Poisson equation and (4.24) for having the correct deep MOND regime. In addition, we recall (3.8) which was imposed in order to be able to expand the two metrics around the same background. Strictly speaking the relation (3.8) is not necessary for the present calculation because we neglected the influence of the background, and (3.8) may be relaxed for some applications. Finally, we note that the result can be simplified by absorbing the remaining constant r g together with the sum α + β into the following redefinitions of the vector field and the effective metric: A µ → r g A µ and g eff µν → (α + β) −2 g eff µν . When this is done, and after redefining also the graviton mass m, we see that we can always choose r g = 1 and α + β = 1 without loss of generality.
Finally the fully reduced form of the action (2.1) reads where we have moved for convenience the mass currents J µ g = ρ g u µ g and J µ f = α β ρ f u µ f to the g and f sectors, where the effective metric g eff µν of bigravity is given by (2.4) but with α + β = 1, and where the kinetic term of the vector field is defined by In addition the coupling constants M 2 g and M 2 f are constrained by (4.19). We can still further impose (3.8) if we insist that the two metrics g and f can be expanded around a common backgroundḡ =f . 11 Finally the theory depends also on the graviton's mass m, hopefully to be related to the observed cosmological constant Λ, and the MOND acceleration scale a 0 . Note that in such unified approach between the cosmological constant and MOND, it is natural to expect that a 0 and Λ should have comparable orders of magnitudes, i.e. a 0 ∼ √ Λ which happens to be in very good agreement with observations.

The decoupling limit
In this section we would like to focus on the matter sector and address the question whether or not the interactions between the matter fields reintroduces the BD ghost [12]. For this we derive the decoupling limit and pay special attention to the helicity-0 mode. Before doing 11 A simple choice satisfying all these requirement is α = β = 1 2 and M 2 g = M 2 f = 1 16π , i.e. each coupling constant takes half the GR value. With this choice the plasma frequency (4.17) becomes in agreement with the finding of the previous model [87]. so, we first restore the broken gauge invariance by introducing the Stueckelberg fields in the f metric (with a, b = 0, 1, 2, 3) where the Stueckelberg fields can be further decomposed into their helicity-0 π and helicity-1 A a counterparts, As next we can take the decoupling limit by sending M g , M f → ∞, while keeping the scale Λ 3 3 = M g m 2 constant. For our purpose, it is enough to keep track of the contributions of the helicity-0 mode π to the matter interactions. The BD ghost is hidden behind the higher derivative terms of the helicity-0 interactions after using all of the covariant equations of motion and constraints. We will therefore neglect the contribution of the helicity-1 field and assume g µν = η µν and where Π µν stands for Π µν ≡ ∂ µ ∂ ν π/Λ 3 3 . In the remaining of this section all indices are raised and lowered with respect to the Minkowski metric η µν . The effective metric in the decoupling limit corresponds to g eff µν →g eff In [88] the required criteria for ghost freedom in the kinetic and potential terms were investigated in detail and hence the new model was constructed in such a way that it fulfils these criteria. In the decoupling limit, the contribution to the equation of motion for the helicity-0 field coming from the potential term is at most second order in derivative for the allowed potentials. Therefore, we can concentrate on the problematic terms coming from the matter interactions. Let us for instance consider the coupling between j μ f and A µ , thus 12 which we can also write as Let us first compute the equation of motion with respect to the vector field, which is simply given by which can be also expressed as For simplicity we will drop the constants λ, r f and M eff and assume that W(X ) = X in the following.

JCAP12(2015)026
As next we can compute the contribution of the matter to the equation of motion with respect to the helicity-0 field, where we made use of Applying one of the derivatives we have The equation of motion for the vector field (5.8) can also be written as which we can use to express Using the equation of motion for the vector field and the conservation equation, the equation for the helicity-0 becomes

JCAP12(2015)026
which we can rewrite as where we have introduced Using the fact that the Christoffel symbol with respect to thef metric is given by 20) and further taking into account the following relation, we see immediately that R ν µσ = 0. Similarly, using the definition of the Christoffel symbol with respect to theg eff metric, 22) and the relation (α + β)δ µ ρ − βΠ µ ρ g ρσ eff (α + β)δ ν σ − βΠ ν σ = η µν , (5.23) we can show that alsoR ν µσ = 0. Thus, the contribution of the matter part to the equation of motion for the helicity-0 mode simplifies to (5.24) Using (5.21) this can be equally written as where j µ * = √ −g eff j μ f is the coordinate current and represents the matter degrees of freedom independent of the metric, while F κµ is also independent of the metric. We can now also apply the remaining derivative in front, The contribution of the matter interactions to the equation of motion for the π field contains higher derivative terms which we are not able to remove by invoking the covariant equations -17 -

JCAP12(2015)026
of motion for the vector field and the dark matter particle. Hence, this reflects the presence of a ghostly degree of freedom in the decoupling limit through the matter interaction term. The scale Λ 3 is not the cut-off scale of this theory anymore, but rather given by the scale of the mass of the introduced ghost m BD . For the precise plasma-like background solution that we considered here, the contributions of the matter interactions and the internal vector field will be at least third order in perturbation at the level of the action. Hence the ghost would enter at non-linear order in perturbations, and does not show up in our linear perturbation analysis of the gravitational polarisation in section 4. This indicates that the mass of the ghost m BD could be large.

Conclusions
In this work we followed the philosophy of combining the different approaches used to tackle the cosmological constant and dark matter problems. The origin of this dark sector constitutes one of the most challenging puzzle of contemporary physics. Here we want essentially to consider them on the same footing. The standard model of cosmology Λ-CDM, despites many observational successes, fails to explain the observed tiny value of the cosmological constant in the presence of large quantum corrections using standard quantum field theory techniques. This has initiated the development of IR modifications of GR, like massive gravity and bigravity. On the other hand, the model Λ-CDM does not account for many observations of dark matter at galactic scales, being unable to explain without fine-tuning the tight correlations between the dark and luminous matter in galaxy halos. Rather, the dark matter phenomenology at galactic scales is in good agreement with MOND [61][62][63].
In the present work we aimed at addressing these two motivations under the same umbrella using a common framework, namely the one of bigravity. An important clue in this respect is the fact that the MOND acceleration scale is of the order of the cosmological constant, a 0 ∼ √ Λ. An additional hope was to be able to promote the MOND formula into a decent relativistic theory in the context of bigravity. For this purpose we followed tightly the same ingredients as in the model proposed in [88], with two species of dark matter particles coupled to the two metrics of bigravity respectively, and linked via an internal vector field. This paper was dedicated to explore the theoretical and phenomenological consistency of this model and verify that it is capable to recover the MOND phenomenology on galactic scales. We first worked out the covariant field equations with a special emphasis on the contributions of the vector field and the dark matter particles. Because of the interaction term between the dark matter particles and the vector field, the stress-energy tensors are not conserved separately, but rather a combination of them, giving rise to a global conservation law. The divergence of the stress-energy tensor with respect to the effective metric g eff to which is coupled the vector field [15,17], is given by the interaction of the dark matter particles with the vector field, which has important consequences for the polarisation mechanism but also for our decoupling limit analysis.
As next we computed the linear field equations around a de Sitter background where the mass term shall play the role of the cosmological constant on large scales. On small scales where the post-Newtonian limit applies, the de Sitter background can be approximated by a flat Minkowski background. Considering a small perturbation of the Minkowski metric and computing the Newtonian limit, we were able to show that the polarisation mechanism works successfully and recovers the MOND phenomenology on galactic scales. We find that this polarisation mechanism is a natural consequence of bimetric gravity since it is made possible -18 -

JCAP12(2015)026
by the constraint equation (linearized Bianchi identity) linking together the two metrics of bigravity. In addition it strongly relies on the internal vector field generated by the coupling between the two species of DM particles living in the two metrics, as it is alive due to an appropriate cancelation between the gravitational force and the internal vector force.
We then studied the decoupling limit of the matter sector in details. We found that precisely because of this coupling to the vector field (whose importance is probably priceless for the polarisation to work) a ghost is reintroduced in the dark matter sector. For the particular solution we have for the dark matter the ghost is higher order in the perturbative expansion. The exact mass of the ghost shall be studied in a future work.