Post-Newtonian parameter $γ$ and the deflection of light in ghost-free massive bimetric gravity
Manuel Hohmann




manuel.hohmann@ut.ee
Laboratory of Theoretical Physics, Institute of Physics, University of Tartu, W. Ostwaldi 1, 50411 Tartu, Estonia


We consider the parametrized post-Newtonian (PPN) limit of ghost-free massive bimetric gravity with two mutually non-interacting matter sectors coupled to the two metrics. Making use of a gauge-invariant differential decomposition of the metric perturbations, we solve the field equations up to the linear PPN order for a static, point-like mass source. From the result we derive the PPN parameter γ for spherically symmetric systems, which describes the gravitational deflection of light by visible matter. By a comparison to its value measured in the solar system we obtain bounds on the parameters of the theory. We further discuss the deflection of light by dark matter and find an agreement with the observed light deflection by galaxies. We finally speculate about a possible explanation for the observed distribution of dark matter in galactic mergers such as Abell 520 and Abell 3827.

Post-Newtonian parameter γ and the deflection of light
in ghost-free massive bimetric gravity
    Manuel Hohmann
    December 30, 2023
=============================================================================================





§ MOTIVATION

Current observations, together with their interpretation according to the standard ΛCDM model of cosmology, indicate that visible matter constitutes only 5% of the total matter content of the universe, while the remaining constituents are given by 26% dark matter and 69% dark energy <cit.>. These dark components of the universe arise purely from phenomenology: dark energy offers a potential explanation for the observed accelerating expansion of the universe <cit.>, while dark matter potentially explains the rotation curves of galaxies <cit.>, the formation of large scale structures in the early universe <cit.> or the lensing and peculiar motion in galaxy clusters <cit.>. Note that all of these observations are based on the gravitational influence of the dark universe on visible matter. Despite significant effort, no direct, non-gravitational interaction of visible matter and dark components has been observed. This raises the question whether any such non-gravitational interaction exists, or whether the coupling between the dark and visible sectors is purely gravitational, or even gravity itself constitutes at least part of the dark sector.

In this article we discuss a model which naturally features a purely gravitational coupling prescription for dark matter, while at the same time being able to accommodate dark energy. This model is based on the idea that the geometry of spacetime, which mediates the gravitational interaction, is described not by a single metric as in general relativity, but by two separate metrics. Further, each matter field couples to only one of these metrics, and there is no direct, non-gravitational coupling between matter fields associated to different metrics. These assumptions imply the existence of two different matter sectors, whose mutual interaction is mediated only by an interaction between the two metrics, so that they appear mutually dark. However attractive in its phenomenology, this idea also leads to potential theoretical issues, since a coupling between them requires at least one of the corresponding gravitons to be massive <cit.>, and such massive gravity theories generally suffer from the existence of a ghost instability <cit.>.

While it has been believed for several decades that the aforementioned ghost instability completely excludes any theories with gravitationally interacting massive spin 2 particles, it has been discovered recently that this is not the case, and that a particular, narrow class of theories avoids the ghost <cit.>; see <cit.> for a number of reviews. The most simple class of such theories indeed features two metric tensors, and allows for two separate, mutually non-gravitationally non-interacting classes of matter fields, each of which couples exclusively to one metric, and which interact with each other only through an interaction between the two metrics <cit.>. The interpretation of one of them as corresponding to dark matter, as we discussed above, has also been studied <cit.>, possibly involving an additional “graviphoton” vector field and reproducing modified Newtonian dynamics (MOND) on galactic scales <cit.>. Note, however, that in contrast to the latter we do not introduce a graviphoton, or aim to model dark matter as a gravitational effect as in MOND. We further remark that in massive gravity theories also the massive graviton is a potential dark matter candidate, an idea which has only recently been considered <cit.>.

Besides providing possible explanations for the observed dark sector of the universe, any viable theory of gravity must of course also pass tests in the solar system. An important tool for testing metric gravity theories with high-precision data from solar system experiments is the parametrized post-Newtonian (PPN) formalism <cit.>. The main idea of the PPN formalism is to express the metric tensor as a perturbation around a flat background, and then expand the perturbation in terms of certain integrals over the gravitating matter distribution. The coefficients of these potentials in the metric perturbation are characteristic for a given gravity theory, and can directly be linked to observable quantities. In this article we focus on a particular PPN parameter, conventionally denoted γ, which has been measured to high precision in numerous solar system experiments <cit.>, in particular through very long baseline interferometry <cit.>, via the Shapiro delay of radio signals <cit.> and using combined observations of the motion of bodies in the solar system <cit.>. To present all observations are in full agreement with the general relativity value γ = 1 <cit.>.

While all of the aforementioned experiments observed the gravitational interaction within the visible sector, it should be noted that the PPN parameter γ has also been determined through the deflection of visible light by galaxies, whose total gravitating matter content contains a significant contribution from dark matter <cit.>. Also these observations, although less precise, are in agreement with the general relativity value γ = 1. It is needless to say that understanding the light deflection by dark matter is essential for a correct interpretation of observations where the dark matter distribution is reconstructed from lensing under the assumption that dark matter deflects light in the same way as visible matter. This is particularly important in the case of galactic mergers, such as the so-called “Bullet Cluster” 1E0657-558 <cit.>, the “Train Wreck Cluster” Abell 520 <cit.>, MACS J0025.4-1222 <cit.> or Abell 3827 <cit.>, where visible and dark matter appear clearly separated from each other. However, it is not a priori clear that this assumption is valid in a theory in which dark and visible matter couple differently to gravity.

As mentioned above, the PPN formalism in its standard form requires a single dynamical metric for the description of gravity. In order to discuss theories with multiple metric tensors, we need an extension of this standard PPN formalism. A possible extension, which features massless and massive gravity modes, but includes only one type of gravitating source matter, has been introduced and applied in <cit.>. A complementary extension to multiple metrics and a corresponding number of matter sectors, but including only massless gravity modes, has been developed and applied in <cit.>. In this article we choose to make use of the latter, and to extend it to also allow calculating the PPN parameter γ for both dark and visible matter in massive gravity. This is the simplest possible extension, and the first step towards a fully general extension of the formalism to massive gravity theories; the latter would allow for a calculation of all PPN parameters.

We remark that the perturbative expansion of the metric in a weak field limit, which is an important ingredient to the PPN formalism, is not always valid in the context of bimetric gravity due to the Vainshtein mechanism <cit.>, and that a full, non-linear treatment is required in order to determine the gravitational dynamics close to the source mass. This non-linear mechanism typically suppresses all deviations from general relativity within a given radius around the source mass, called the Vainshtein radius. A perturbative treatment is valid only outside this radius. We will not discuss the Vainshtein mechanism in this article, and restrict ourselves to the case of theories in which the Vainshtein radius is sufficiently small so that the perturbative treatment is valid on solar system scales and above.

The outline of this article is as follows. In section <ref> we briefly review the action and field equations of ghost-free massive bimetric gravity. We then perform a perturbative expansion of these field equations in section <ref>, using an adapted version of the PPN formalism. We further simplify the obtained equations using gauge-invariant perturbation theory in section <ref>. This will yield us a set of equations, which we will solve for a static, point-like mass source in section <ref>, and thus determine the effective Newtonian gravitational constant and PPN parameter γ. We will connect our result to observations, in particular of the deflection of light, in section <ref>. We end with a conclusion in section <ref>. A few lengthy calculations are displayed in the appendix. In appendix <ref> we derive the linearized interaction potential connecting the two metrics. In appendix <ref> we list derivatives of the Yukawa potential. We show how to check our solution of the field equations in appendix <ref>.



§ ACTION AND FIELD EQUATIONS

In this section we start our discussion of the post-Newtonian limit of bimetric gravity with a brief review of its action and gravitational field equations, which are derived by variation with respect to the two metrics. We then trace-reverse the field equations, as this will be more convenient when we construct their solution. These trace-reversed field equations will be the main ingredient for our calculation.

The starting point for our derivation is the action functional

    S = _Md^4x[m_g^2/2√(- g)R^g + m_f^2/2√(- f)R^f - m^4√(- g)∑_n = 0^4β_ne_n(√(g^-1f))
    
    + √(- g)ℒ_m^g(g,Φ^g) + √(- f)ℒ_m^f(f,Φ^f)]

for two metric tensors g_μν, f_μν and two sets of matter fields Φ^g,f, each of which couples to only one metric tensor, and between which there is no direct, non-gravitational interaction. One may thus interpret, e.g., Φ^g as visible matter, constituted by the standard model fields and governed by the standard matter Lagrangian ℒ^g, and Φ^f as dark matter, constituted by a distinct set of fields with possibly different structure of the Lagrangian ℒ^f. However, we will not make any assumptions on the constituting fields of the two matter types here, or on their Lagrangians, as these will not be relevant for our calculation.

Note that since there are two metrics, there is no canonical prescription for raising or lowering tensor indices. Therefore, we will not raise or lower indices automatically, but provide definitions for all tensor fields with fixed index positions. In the action (<ref>) this applies to the Ricci scalars, each of which is defined solely through its corresponding metric, such that they are related to the Ricci tensors by

    R^g = g^μνR^g_μν ,    R^f = f^μνR^f_μν .

Note further the appearance of the (1,1) tensor field g^-1f, which we assume to have a square root A such that

    A^μ_σA^σ_ν = g^μσf_σν .

This is certainly the case in a sufficiently small neighborhood of the flat proportional background metrics g_μν = η_μν, f_μν = c^2η_μν, which we will henceforth consider. The functions e_0, …, e_4 in the action are the matrix invariants

    e_k(A) = A^μ_1_[μ_1⋯ A^μ_k_μ_k] = 1/k!(4 - k)!ϵ^μ_1 ⋯μ_k λ_1 ⋯λ_4 - kϵ_ν_1 ⋯ν_k λ_1 ⋯λ_4 - kA^ν_1_μ_1⋯ A^ν_k_μ_k

of this square root, while their coefficients β_0, …, β_4 are constant, dimensionless parameters to the action. The remaining parameters are the Planck masses m_g, m_f and the interaction mass m, all of which are of mass dimension. Any choice of the constant parameters determines a particular action, and hence a particular theory.

By variation with respect to the metric tensors we obtain the gravitational field equations

    m_g^2(R^g_μν - 1/2g_μνR^g) + m^4V^g_μν = T^g_μν ,   
    m_f^2(R^f_μν - 1/2f_μνR^f) + m^4V^f_μν = T^f_μν .

Here we have defined the energy-momentum tensors as usual through

    T^g_μν = -2/√(- g)δ(√(- g)ℒ_m^g(g,Φ_g))/δ g^μν ,   
    T^f_μν = -2/√(- f)δ(√(- f)ℒ_m^f(f,Φ_f))/δ f^μν .

The potential terms V^g,f_μν are given by

    V^g_μν = g_μρ∑_n = 0^3(-1)^nβ_nY_n^ρ_ν(A) ,   
    V^f_μν = f_μρ∑_n = 0^3(-1)^nβ_4 - nY_n^ρ_ν(A^-1) ,

where the functions Y_0, …, Y_3 are defined as

    Y_n(A) = ∑_k = 0^n(-1)^ke_k(A)A^n - k

and analogously for Y_n(A^-1). We remark that the action (<ref>), and hence also the field equations (<ref>), are fully symmetric under a simultaneous exchange g_μν↔ f_μν and β_n ↔β_4 - n.

While we could work directly with the field equations (<ref>), it turns out to be more convenient to use the trace-reversed equations instead, which read

    m_g^2R^g_μν + m^4V̅^g_μν = T̅^g_μν ,   
    m_f^2R^f_μν + m^4V̅^f_μν = T̅^f_μν .

Here we have trace-reversed each term with its corresponding metric, i.e., we have applied the definitions


    V̅^g_μν   = V^g_μν - 1/2g_μνg^ρσV^g_ρσ ,    T̅^g_μν   = T^g_μν - 1/2g_μνg^ρσT^g_ρσ ,
    V̅^f_μν   = V^f_μν - 1/2f_μνf^ρσV^f_ρσ ,    T̅^f_μν   = T^f_μν - 1/2f_μνf^ρσT^f_ρσ .


The field equations (<ref>) are the equations we will be working with during the remainder of this article. In order to calculate the post-Newtonian limit, we will need a perturbative expansion of these equations. This will be done in the next section.



§ POST-NEWTONIAN APPROXIMATION

We now come to a perturbative expansion of the field equations (<ref>) displayed in the previous section. For this purpose, we first briefly review the notion of velocity orders in section <ref>, and label the relevant components of the metric and energy-momentum tensors. We then discuss the metric ansatz, essentially following the construction developed in <cit.>, in section <ref>. Finally, we apply these constructions to the field equations under consideration. We derive and solve the field equations at the zeroth velocity order in section <ref>, and derive the second order equations in section <ref>.



 §.§ Expansion in velocity orders

A central ingredient of the PPN formalism is the assumption that the gravitating source matter is constituted by a perfect fluid. Since there are two different types of matter Φ^g,f in the theory we consider, which interact only gravitationally, we apply this assumption to each of them. Their energy-momentum tensors therefore take the form


    T^g μν   = (ρ^g + ρ^gΠ^g + p^g)u^g μu^g ν + p^gg^μν ,
    
    T^f μν   = (ρ^f + ρ^fΠ^f + p^f)u^f μu^f ν + p^ff^μν ,


with rest energy densities ρ^g,f, specific internal energies Π^g,f, pressures p^g,f and four-velocities u^g,f μ. Note that the four-velocities are normalized with their corresponding metrics,

    u^g μu^g νg_μν = u^f μu^f νf_μν = -1 .

We further assume that the source matter is slow-moving within our chosen frame of reference, so that the velocity components satisfy

    v^g,f i = u^g,f i/u^g,f 0≪ 1 .

We then assign orders of magnitude 𝒪(n) ∝ |v⃗|^n to all dynamical quantities. For the matter variables we assign ρ^g,f∼Π^g,f∼𝒪(2) and p^g,f∼𝒪(4), based on their values for the matter constituting the solar system. For the metrics we assume a small perturbation around a flat, proportional background solution, where we expand the perturbation in velocity orders in the form


    g_μν = η_μν + h_μν   = η_μν + h^(1)_μν + h^(2)_μν + h^(3)_μν + h^(4)_μν + 𝒪(5) ,
    
    c^-2f_μν = η_μν + e_μν   = η_μν + e^(1)_μν + e^(2)_μν + e^(3)_μν + e^(4)_μν + 𝒪(5)


with constant c > 0. Not all of these components are relevant for the post-Newtonian limit, while others vanish due to symmetries and conservation laws. The only non-vanishing components that are relevant for our calculation of the PPN parameter γ in this article are

    h^(2)_00 ,    h^(2)_ij ,    e^(2)_00 ,    e^(2)_ij .

Further, we only consider quasi-static solutions, so that changes of the metric are induced only by the motion of the source matter. We therefore assign another velocity order ∂_0 to any time derivative. We finally assume that the source matter is located in a bounded region and that the metrics are asymptotically flat, so that the metric perturbations and their derivatives vanish at infinity.

Since we are interested only in the second order metric perturbations (<ref>), and hence need to solve the field equations only up to the second velocity order, it is also sufficient to expand the energy-momentum tensors (<ref>) to the second velocity order. The only relevant components are

    T^g(2) 00 = ρ^g ,    T^f(2) 00 = ρ^f/c^2 ,    T^g(2) ij = T^f(2) ij = 0 .

In order to use them in the field equations (<ref>), we need to lower their indices with their corresponding metrics, which yields

    T^g(2)_00 = ρ^g ,    T^f(2)_00 = c^2ρ^f ,    T^g(2)_ij = T^f(2)_ij = 0 .

Finally, we also need to trace-reverse these terms with their corresponding metrics, from which we obtain

    T̅^g(2)_00 = 1/2ρ^g ,   T̅^f(2)_00 = c^2/2ρ^f ,   T̅^g(2)_ij = 1/2ρ^gδ_ij ,   T̅^f(2)_ij = c^2/2ρ^fδ_ij .

These terms will enter the trace-reversed field equations (<ref>).



 §.§ Metric ansatz and PPN parameters

Another important ingredient of the PPN formalism is an ansatz for the metric in terms of potentials, which are integrals over the source matter distribution. Their coefficients in the metric are observable quantities which allow a characterization of the gravity theory under examination. For single metric theories there is a standard form for this PPN metric ansatz <cit.>. Here we use a generalization to multimetric theories developed in <cit.>. Our ansatz for the second order metric perturbation reads


    h^(2)_00   = -α^ggχ^g - α^gfχ^f ,
    
    h^(2)_ij   = 2θ^ggχ^g_,ij + 2θ^gfχ^f_,ij - [(γ^gg + θ^gg)χ^g + (γ^gf + θ^gf)χ^f]δ_ij ,
    
    e^(2)_00   = -α^fgχ^g - α^ffχ^f ,
    
    e^(2)_ij   = 2θ^fgχ^g_,ij + 2θ^ffχ^f_,ij - [(γ^fg + θ^fg)χ^g + (γ^ff + θ^ff)χ^f]δ_ij ,


where = ∂^i∂_i and indices are raised and lowered with the flat metric η_μν. The PPN potentials we have introduced here are second order derivatives of the superpotentials

    χ^g(t, x⃗) = -∫ρ^g(t, x⃗') |x⃗ - x⃗'| d^3x' ,   χ^f(t, x⃗) = -c^3∫ρ^f(t, x⃗') |x⃗ - x⃗'| d^3x' .

The definition of χ^f contains a factor c^3, which originates from the volume element of the spatial part of the unperturbed contribution c^2η_μν of the metric f_μν. In the PPN metric we further have twelve PPN parameters α^g,f g,f, γ^g,f g,f, θ^g,f g,f. Note that these are not independent, as we have not yet fixed a gauge for the metric. The gauge freedom allows us to apply a diffeomorphism generated by a vector field ξ^μ, provided that it preserves the perturbation ansatz (<ref>). This means that the vector field ξ^μ must be of the same order as the metric perturbations. Recall that under a diffeomorphism the metrics change according to

    δ_ξg_μν = (ℒ_ξg)_μν = 2g_σ(μ∇^g_ν)ξ^σ ,   δ_ξf_μν = (ℒ_ξf)_μν = 2f_σ(μ∇^f_ν)ξ^σ ,

where ℒ denotes the Lie derivative. At the linear perturbation level, which is sufficient for our calculation here, this yields the transformation of the metric perturbations

    δ_ξh_μν = δ_ξe_μν = 2η_σ(μ∂_ν)ξ^σ .

Further demanding consistency with the PPN metric ansatz (<ref>) we find that the only allowed and relevant vector field is of second velocity order and can be written as <cit.>

    ξ_0 = 0 ,   ξ_i = λ^gχ^g_,i + λ^fχ^f_,i ,

with two constants λ^g,f, and where we have defined ξ_μ = η_μνξ^ν. Under a diffeomorphism generated by this vector field the metric perturbations change by

    δ_ξh^(2)_00 = δ_ξe^(2)_00 = 0 ,   δ_ξh^(2)_ij = δ_ξe^(2)_ij = 2λ^gχ^g_,ij + 2λ^fχ^f_,ij .

By choosing λ^g = -θ^gg and λ^f = -θ^ff we can always eliminate these two PPN parameters from the metric ansatz. In the remainder of our calculation we will adopt this gauge, in which θ^gg = θ^ff = 0, as this turns out to be compatible with the standard PPN gauge for single metric theories <cit.>.



 §.§ Background solution

Recall from section <ref> that we have expanded the metrics around a flat Minkowski background, as usual in the PPN formalism. For the PPN formalism to be applicable in this form, it is necessary that this background is a solution of the field equations at the zeroth velocity order. Since both the Ricci tensors R^g,f(0)_μν and the energy-momentum tensors T̅^g,f(0)_μν at the zeroth velocity order vanish, these simply reduce to

    m^4V̅^g(0)_μν = 0 ,    m^4V̅^f(0)_μν = 0 ,

where we assume m > 0. In order to determine the potential terms V̅^g,f(0)_μν, and also the second velocity order in the next section, it is useful to first linearize the potential in the metric perturbations. Since this is a rather lengthy procedure, we have deferred it to appendix <ref>. Here we make use of the result (<ref>), from which we read off the zeroth order contribution


    V̅^g(0)_μν   = -(β̃_0 + 3β̃_1 + 3β̃_2 + β̃_3)η_μν = -(β_0 + 3cβ_1 + 3c^2β_2 + c^3β_3)η_μν ,
    V̅^f(0)_μν   = -(β̃_1 + 3β̃_2 + 3β̃_3 + β̃_4)c^-2η_μν = -(β_1 + 3cβ_2 + 3c^2β_3 + c^3β_4)c^-1η_μν ,


where we used the abbreviations β̃_k = c^kβ_k. We require that these equations, which are polynomial in c, possess at least one common positive solution c > 0. Note that a particular fixed c solves both equations if and only if the parameters in the action (<ref>) satisfy

    β_0 = -3cβ_1 - 3c^2β_2 - c^3β_3 ,   β_4 = -c^-3β_1 - 3c^-2β_2 - 3c^-1β_3 .

The condition that the background equations are solved by proportional flat metrics therefore completely determines the two parameters β_0 and β_4 in the action in terms of a new free parameter c > 0. In the following we will therefore replace β_0 and β_4, and hence β̃_0 and β̃_4, using

    β̃_0 = -3β̃_1 - 3β̃_2 - β̃_3 ,   β̃_4 = -β̃_1 - 3β̃_2 - 3β̃_3 ,

and keep β_1, β_2, β_3 and c as free parameters of the class of theories we discuss.



 §.§ Second order field equations

For the remainder of our calculation and in order to determine the second order metric perturbations we will need to expand the field equations (<ref>) to the second velocity order. The only relevant components are given by


    m_g^2R^g(2)_00 + m^4V̅^g(2)_00   = T̅^g(2)_00 ,    
    m_g^2R^g(2)_ij + m^4V̅^g(2)_ij   = T̅^g(2)_ij ,
    
    m_f^2R^f(2)_00 + m^4V̅^f(2)_00   = T̅^f(2)_00 ,    
    m_f^2R^f(2)_ij + m^4V̅^f(2)_ij   = T̅^f(2)_ij .


We have already calculated the necessary components (<ref>) of the energy-momentum tensor at the second velocity order. The components of the Ricci tensor are easily obtained and yield the standard textbook result <cit.>


    R^g(2)_00   = -1/2 h^(2)_00 ,    
    R^g(2)_ij   = -1/2( h^(2)_ij - h^(2)_00,ij + h^(2)_kk,ij - h^(2)_ik,jk - h^(2)_jk,ik) ,
    
    R^f(2)_00   = -1/2 e^(2)_00 ,    
    R^f(2)_ij   = -1/2( e^(2)_ij - e^(2)_00,ij + e^(2)_kk,ij - e^(2)_ik,jk - e^(2)_jk,ik) .


Finally, we also need the second velocity order contribution from the potential terms. Using the result (<ref>) derived in appendix (<ref>) one finds the components


    V̅^g(2)_00   = 1/4β̃(3h^(2)_00 - 3e^(2)_00 - h^(2)_ii + e^(2)_ii) ,
    V̅^g(2)_ij   = 1/4β̃[2h^(2)_ij - 2e^(2)_ij + (h^(2)_kk - e^(2)_kk - h^(2)_00 + e^(2)_00)δ_ij] ,
    V̅^f(2)_00   = 1/4c^2β̃(3e^(2)_00 - 3h^(2)_00 - e^(2)_ii + h^(2)_ii) ,
    V̅^f(2)_ij   = 1/4c^2β̃[2e^(2)_ij - 2h^(2)_ij + (e^(2)_kk - h^(2)_kk - e^(2)_00 + h^(2)_00)δ_ij] ,


where we introduced the abbreviation

    β̃ = β̃_1 + 2β̃_2 + β̃_3 .

These are the field equations we will be using during the remainder of this article. However, directly working with these equations poses two difficulties. First, the field equations possess a gauge freedom, as discussed in section <ref>, and so the solution will be unique only after gauge fixing. Second, the equations turn out to be involved and cumbersome to solve due to the mixing of tensor components. Both of these difficulties can be solved straightforwardly by performing a gauge-invariant differential decomposition of the metric perturbations. We will detail this formalism in the next section.



§ GAUGE-INVARIANT DIFFERENTIAL DECOMPOSITION

In the previous section we have performed an expansion of the gravitational field equations (<ref>) into velocity orders and obtained the second order equations (<ref>). Instead of solving them directly for the metric perturbations (<ref>), we will first bring them into a significantly simpler form in this section. For this purpose, we employ the formalism of gauge-invariant perturbations, which is well-known from cosmology <cit.>. We apply this procedure in several steps. First, we decompose the metric perturbations into gauge-invariant potentials in section <ref>. In section <ref> we further decompose these potentials into velocity orders as required by the PPN formalism. Using the expressions obtained, we then decompose the Ricci tensors (<ref>) in section <ref>, the potentials (<ref>) in section <ref> and the energy-momentum tensors (<ref>) in section <ref>. This will finally yield us a full decomposition of the field equations (<ref>) in section <ref>.



 §.§ Decomposition of the metrics

We start with a differential decomposition of the metric perturbations. Using the split into time and space components, we introduce the decomposition


    h_00   = -2ϕ^g,    
    h_0i   = ∂_iB^g + B^g_i,    
    h_ij   = -2ψ^gδ_ij + 2_ijE^g + 4∂_(iE^g_j) + 2E^g_ij ,
    
    e_00   = -2ϕ^f,    
    e_0i   = ∂_iB^f + B^f_i,    
    e_ij   = -2ψ^fδ_ij + 2_ijE^f + 4∂_(iE^f_j) + 2E^f_ij


into four scalars ϕ^g,f, ψ^g,f, B^g,f, E^g,f, two divergence-free vectors B^g,f_i, E^g,f_i and one trace-free, divergence-free tensor E^g,f_ij. Here _ij denotes the trace-free second derivative _ij = ∂_i∂_j - 1/3δ_ij. From these quantities we further derive the potentials

    I_1^g,f = ϕ^g,f + ∂_0B^g,f - ∂_0^2E^g,f ,   
    I_3^g,f = B^g,f ,   
    I_4^g,f = E^g,f ,
    
    I_2^g,f = ψ^g,f + 1/3 E^g,f ,   
    I^g,f_i = B^g,f_i - 2∂_0E^g,f_i ,   
    I'^g,f_i = E^g,f_i ,   
    I^g,f_ij = E^g,f_ij .

The advantage of using these potentials becomes apparent when we consider gauge transformations of the metric, i.e., diffeomorphisms generated by a vector field ξ^μ which preserve the perturbation ansatz (<ref>) as discussed in section <ref>. Here we introduce a differential decomposition for ξ_μ = η_μνξ^ν of the form

    ξ_0 = X ,   ξ_i = ∂_iX' + X_i

into two scalars X, X' and one divergence-free vector X_i. One now easily computes from the decomposition (<ref>) the transformations

    δ_ξϕ^g,f = -∂_0X ,   δ_ξψ^g,f = -1/3 X' ,   δ_ξB^g,f = ∂_0X' + X ,   δ_ξE^g,f = X' ,
    δ_ξB^g,f_i = ∂_0X_i ,   δ_ξE^g,f_i = 1/2X_i ,   δ_ξE^g,f_ij = 0 .

The potentials (<ref>) hence transform as

    δ_ξI_1^g,f = δ_ξI_2^g,f = 0 ,   δ_ξI_3^g,f = ∂_0X' + X ,   δ_ξI_4^g,f = X' ,
    δ_ξI^g,f_i = 0 ,   δ_ξI'^g,f_i = 1/2X_i ,   δ_ξI^g,f_ij = 0 .

Finally, defining the linearly related potentials

    I_1^± = I_1^g ± I_1^f ,   
    I_2^± = I_2^g ± I_2^f ,   
    I_3^± = I_3^g ± I_3^f ,   
    I_4^± = I_4^g ± I_4^f ,
    
    I^±_i = I^g_i ± I^f_i ,   
    I'^±_i = I'^g_i ± I^f_i ,   
    I^±_ij = I^g_ij± I^f_ij ,

we see that the six scalar potentials I_1^±, I_2^±, I_3^-, I_4^-, the three vectors I^±_i, I'^-_i and the two tensors I^±_ij are invariant under gauge transformations, while the remaining two scalars I_3^+, I_4^+ and the vector I'^+_i are pure gauge degrees of freedom corresponding to the two scalars and the vector constituting the diffeomorphisms. The only physical degrees of freedom are the gauge invariant potentials. Since the gravitational field equations are derived from a diffeomorphism invariant action, we can fully express them in terms of these gauge invariants. In the following we will do so by introducing a suitable differential decomposition of the Ricci tensors, potentials and energy-momentum tensors.



 §.§ Gauge invariant potentials and velocity orders

Recall from section <ref> that we have decomposed the metric perturbations into velocity orders and that only the components (<ref>) are relevant for the calculation we present in this article. We now apply the decomposition into velocity orders to the gauge invariant potentials above in order to determine which of them will be relevant for our calculation. A comparison of the relevant components (<ref>) with the differential decomposition (<ref>) shows that only the quantities

    ϕ^g,f(2) ,   ψ^g,f(2) ,    E^g,f(2) ,    E^g,f(2)_i ,    E^g,f(2)_ij

at the second velocity order will be relevant. Using the relations (<ref>), while taking into account that time derivatives are weighted with an additional velocity order 𝒪(1), then yields the relevant potentials

    I_1^g,f(2) = ϕ^g,f(2) ,   
    I_2^g,f(2) = ψ^g,f(2) + 1/3 E^g,f(2) ,   
    I_4^g,f(2) = E^g,f(2) ,
    
    I'^g,f(2)_i = E^g,f(2)_i ,   
    I^g,f(2)_ij = E^g,f(2)_ij .

Finally, transitioning to the linearly related potentials (<ref>) shows that the only relevant gauge-invariant potentials are the five scalars I_1^±(2), I_2^±(2), I_4^-(2), the vector I'^-(2)_i and the two tensors I^±(2)_ij, while the scalar I_4^+(2) and the vector I'^+(2)_i are pure gauge quantities. From the former we can now calculate the relevant components of the Ricci tensor and the potential.



 §.§ Decomposition of the Ricci tensors

We now perform a differential decomposition of the Ricci tensors, similar to the differential decomposition (<ref>) of the metric introduced above. Here we use the defining relations

    R^g,f_00 = K_1^g,f ,   
    R^g,f_0i = ∂_iK_3^g,f + K^g,f_i ,   
    R^g,f_ij = 1/3K_2^g,fδ_ij + _ijK_4^g,f + 2∂_(iK'^g,f_j) + K^g,f_ij .

From this definition and the second order field equations (<ref>) follows that the only relevant components for our calculation are given by

    K_1^g,f(2) =  I_1^g,f(2) ,   
    K_2^g,f(2) = 4 I_2^g,f(2) -  I_1^g,f(2) ,
    
    K_4^g,f(2) = I_2^g,f(2) - I_1^g,f(2) ,   
    K'^g,f(2)_i = 0 ,   
    K^g,f(2)_ij = - I^g,f(2)_ij .

Comparing these expressions with the gauge transformations (<ref>) we see that they contain only gauge-invariant potentials, as expected from the fact that they originate from a diffeomorphism invariant action.



 §.§ Decomposition of the potentials

For the trace-reversed potentials we proceed in full analogy to the decomposition (<ref>) of the Ricci tensors. Here we use the decomposition

    V̅^g,f_00 = U_1^g,f ,   V̅^g,f_0i = ∂_iU_3^g,f + U^g,f_i ,   V̅^g,f_ij = 1/3U_2^g,fδ_ij + _ijU_4^g,f + 2∂_(iU'^g,f_j) + U^g,f_ij .

From the second order field equations (<ref>) we read off that the relevant components are given by


    U_1^g(2) = -c^2U_1^f(2)   = -1/2β̃(3I_1^-(2) - 3I_2^-(2) +  I_4^-(2)) ,
    
    U_2^g(2) = -c^2U_2^f(2)   = 1/2β̃(3I_1^-(2) - 15I_2^-(2) + 5 I_4^-(2)) ,
    
    U_4^g(2) = -c^2U_4^f(2)   = β̃I_4^-(2) ,
    
    U'^g(2)_i = -c^2U'^f(2)_i    = β̃I'^-(2)_i ,
    
    U^g(2)_ij = -c^2U^f(2)_ij   = β̃I^-(2)_ij .


Again we see that these depend only on gauge invariant potentials, as expected.



 §.§ Decomposition of the energy-momentum tensors

We finally also need to perform a differential decomposition of the energy-momentum tensors. Following the same prescription as for the Ricci tensors and the potentials we define

    T̅^g,f_00 = Q_1^g,f ,   T̅^g,f_0i = ∂_iQ_3^g,f + Q^g,f_i ,   T̅^g,f_ij = 1/3Q_2^g,fδ_ij + _ijQ_4^g,f + 2∂_(iQ'^g,f_j) + Q^g,f_ij .

For the expressions (<ref>) for the second order trace-reversed energy-momentum tensors of the perfect fluid then follow the relevant components

    Q_1^g(2) = 1/2ρ^g ,   
    Q_1^f(2) = 1/2c^2ρ^f ,   
    Q_2^g(2) = 3/2ρ^g ,   
    Q_2^f(2) = 3/2c^2ρ^f ,
    
    Q_4^g,f(2) = 0 ,   
    Q'^g,f(2)_i = 0 ,   
    Q^g,f(2)_ij = 0 .

These are all expressions we need for the field equations (<ref>) at the second velocity order.



 §.§ Decomposition of the field equations

We now have all expressions at hand which are necessary to perform a differential decomposition of the second order field equations (<ref>) and to fully express them in terms of gauge-invariant quantities. It is an important feature of the differential decomposition that it is unique and bijective under the boundary conditions mentioned in section <ref>, which imply that all metric perturbations and their derivatives vanish at infinity. It thus follows that the field equations (<ref>) are equivalent to the decomposed field equations


    m_g^2K_1^g(2) + m^2U_1^g(2)   = Q_1^g(2) ,    
    m_f^2K_1^f(2) + m^2U_1^f(2)   = Q_1^f(2) ,
    
    m_g^2K_2^g(2) + m^2U_2^g(2)   = Q_2^g(2) ,    
    m_f^2K_2^f(2) + m^2U_2^f(2)   = Q_2^f(2) ,
    
    m_g^2K_4^g(2) + m^2U_4^g(2)   = Q_4^g(2) ,    
    m_f^2K_4^f(2) + m^2U_4^f(2)   = Q_4^f(2) ,
    
    m_g^2K'^g(2)_i + m^2U'^g(2)_i    = Q'^g(2)_i ,    
    m_f^2K'^f(2)_i + m^2U'^f(2)_i    = Q'^f(2)_i ,
    
    m_g^2K^g(2)_ij + m^2U^g(2)_ij   = Q^g(2)_ij ,    
    m_f^2K^f(2)_ij + m^2U^f(2)_ij   = Q^f(2)_ij .


We can now insert the expressions for the differential components of the Ricci tensors, the potentials and the energy-momentum tensors which we derived above. We start with the trace-free, divergence-free tensor equations (<ref>). Inserting the components K^g,f(2)_ij, U^g,f(2)_ij, Q^g,f(2)_ij yields the equations


    -m_g^2/2( I^+(2)_ij +  I^-(2)_ij) + m^4β̃ I^-(2)_ij   = 0 ,
    
    -m_f^2/2( I^+(2)_ij -  I^-(2)_ij) - m^4β̃/c^2I^-(2)_ij   = 0 .


Note that together with the boundary conditions they yield the trivial solution I^±(2)_ij = 0. We then continue with the divergence-free vector equations (<ref>). Using the expressions for K'^g,f(2)_i, U'^g,f(2)_i, Q'^g,f(2)_i we obtain

    m^4β̃I'^-(2)_i = 0 ,    -m^4β̃/c^2I'^-(2)_i = 0 .

These equations are equivalent as a consequence of the Bianchi identities, which follow from the diffeomorphism invariance of the action (<ref>). Also these equations yield a trivial solution I'^-(2)_i = 0. We are thus left with the scalar equations (<ref>), (<ref>) and (<ref>), which take the form


    1/2ρ^g    = m_g^2/2( I_1^+(2) +  I_1^-(2)) - m^4β/2(3I_1^-(2) - 3I_2^-(2) +  I_4^-(2)) ,
    c^2/2ρ^f    = m_f^2/2( I_1^+(2) -  I_1^-(2)) + m^4β/2c^2(3I_1^-(2) - 3I_2^-(2) +  I_4^-(2)) ,
    3/2ρ^g    = m_g^2/2(4 I_2^+(2) + 4 I_2^-(2) -  I_1^+(2) -  I_1^-(2)) + m^4β/2(3I_1^-(2) - 15I_2^-(2) + 5 I_4^-(2)) ,
    3c^2/2ρ^f    = m_f^2/2(4 I_2^+(2) - 4 I_2^-(2) -  I_1^+(2) +  I_1^-(2)) - m^4β/2c^2(3I_1^-(2) - 15I_2^-(2) + 5 I_4^-(2)) ,
    
    0    = m_g^2/2(I_2^+(2) + I_2^-(2) - I_1^+(2) - I_1^-(2)) + m^4β I_4^-(2) ,
    
    0    = m_f^2/2(I_2^+(2) - I_2^-(2) - I_1^+(2) + I_1^-(2)) - m^4β/c^2I_4^-(2) .


Note that also these equations are not independent, but are related to each other as a consequence of the Bianchi identities. Indeed, one easily checks that

    K_2^g,f(2) - 3K_1^g,f(2) - 4 K_4^g,f(2) = 0 ,   
    Q_2^g,f(2) - 3Q_1^g,f(2) - 4 Q_4^g,f(2) = 0 ,
    
    U_2^g(2) - 3U_1^g(2) - 4 U_4^g(2) = 6m^4β̃( I_1^-(2) - 2 I_2^-(2)) = -c^2(U_2^f(2) - 3U_1^f(2) - 4 U_4^f(2)) ,

which shows that the corresponding linear combinations of the scalar equations become identical. Symbolically, this can be written as

    (<ref>) - 3(<ref>) - 4(<ref>) = -c^2[(<ref>) - 3(<ref>) - 4(<ref>)] .

Hence, one of the equations (<ref>) is redundant and can be omitted. The remaining five equations then determine the five gauge-invariant scalar potentials I_1^±(2), I_2^±(2), I_4^-(2). We will solve these equations in the following section for the special case of a static point mass source.



§ STATIC SPHERICALLY SYMMETRIC SOLUTION

Using the gauge invariant field equations (<ref>) derived in the preceding section we are now in the position to construct an explicit solution. The starting point will be a point mass, which we discuss in section <ref>. We will then determine a solution for the gauge invariant potentials in section <ref>. From these we will derive the metric components in section <ref>, reversing the procedure detailed in section <ref>. By comparison with the metric ansatz (<ref>) we read off the PPN parameters in section <ref>. We finally discuss a few limiting cases in section <ref>.



 §.§ Point-mass source

The matter source we consider for our solution is a static point mass located at the origin of our coordinate system, which is constituted by masses M^g and M^f with respect to the two matter sectors. Invoking the interpretation of the matter sectors as visible and dark matter, this would correspond to a source containing both visible and dark matter, unless one of the masses vanishes. This choice is the most general one, and includes the physically relevant case of a galaxy with a dark matter component, as we discuss later in section <ref>. A source of this type is characterized by the matter variables

    ρ^g = M^gδ(x⃗) ,   ρ^f = M^fδ(x⃗)/c^3 ,   Π^g,f = 0 ,    p^g,f = 0 ,    v^g,f_i = 0 ,

where we have normalized the delta function in ρ^f with the spatial volume element c^3 of the unperturbed metric f^(0)_μν = c^2η_μν. Note that this factor cancels the volume element in the corresponding superpotential (<ref>). Using isotropic spherical coordinates, the superpotentials thus read

    χ^g,f = -M^g,fr .

For later convenience we also list the second order derivatives of the superpotentials, which take the form

    χ^g,f_,ij = M^g,f(x_ix_j/r^3 - δ_ij/r) ,   χ^g,f = -2M^g,f/r .

These will be used when we read off the PPN parameters in section <ref>.



 §.§ Gauge-invariant potentials

We will now determine the gauge-invariant potentials I_1^±(2), I_2^±(2), I_4^-(2) by solving the scalar part (<ref>) of the field equations at the second velocity order, where we assume the matter source given by the point mass introduced above. It will turn out to be convenient to use rescaled mass units

    m̃_g = m_g ,   m̃_f = cm_f   M̃^g = M^g ,   M̃^f = cM^f .

We then start with the purely algebraic equations (<ref>) and (<ref>). Using the definitions above these take the form


    0    = m̃_g^2/2(I_2^+(2) + I_2^-(2) - I_1^+(2) - I_1^-(2)) + m^4β̃I_4^-(2) ,
    
    0    = m̃_f^2/2(I_2^+(2) - I_2^-(2) - I_1^+(2) + I_1^-(2)) - m^4β̃I_4^-(2) .


Here we choose to solve these equations for the potentials I_2^±(2). The solutions are given by

    I_2^±(2) = I_1^±(1) - (1/m̃_g^2∓1/m̃_f^2)m^4β̃I_4^-(2) .

We can use this relation to eliminate I_2^±(2) from the remaining equations. Using the linear combination (<ref>), which together with the boundary conditions yields

    I_1^-(2) = 2I_2^-(2) ,

we can then solve for I_4^-(2) and obtain the solution

    I_4^-(2) = 1/2μ^2I_1^-(2) ,

where we have defined the mass parameter

    μ = m^2√(β̃(1/m̃_f^2 + 1/m̃_g^2)) .

We now take a suitable linear combination of the scalar equations (<ref>) and (<ref>), so that the terms involving I_1^+(2) cancel. Eliminating I_2^-(2) and I_4^-(2) with the relations (<ref>) and (<ref>) we obtain

    I_1^-(2) - μ^2I_1^-(2) = 2/3(M̃^g/m̃_g^2 - M̃^f/m̃_f^2)δ(x⃗) ,

which is a screened Poisson equation for I_1^-(2). The solution is given by

    I_1^-(2) = -(M̃^g/m̃_g^2 - M̃^f/m̃_f^2)e^-μ r/6π r .

From the relations (<ref>) and (<ref>) then immediately follows

    I_2^-(2) = -(M̃^g/m̃_g^2 - M̃^f/m̃_f^2)e^-μ r/12π r ,   
    I_4^-(2) = -(M̃^g/m̃_g^2 - M̃^f/m̃_f^2)e^-μ r/12πμ^2r .

Inserting these results into the scalar equation (<ref>), we obtain an equation for I_1^+(2), which is most conveniently expressed as

    (I_1^+(2) + m̃_g^2 - m̃_f^2/m̃_g^2 + m̃_f^2I_1^-(2)) = M̃^g + M̃^f/m̃_g^2 + m̃_f^2δ(x⃗) .

This is an ordinary Poisson equation and one immediately reads off the solution

    I_1^+(2) = -M̃^g + M̃^f/m̃_g^2 + m̃_f^21/4π r - m̃_g^2 - m̃_f^2/m̃_g^2 + m̃_f^2I_1^-(2) = -M̃^g + M̃^f/m̃_g^2 + m̃_f^21/4π r + m̃_g^2 - m̃_f^2/m̃_g^2 + m̃_f^2(M̃^g/m̃_g^2 - M̃^f/m̃_f^2)e^-μ r/6π r .

Finally, making use of the relation (<ref>) yields

    I_2^+(2) = -M̃^g + M̃^f/m̃_g^2 + m̃_f^21/4π r + m̃_g^2 - m̃_f^2/m̃_g^2 + m̃_f^2(M̃^g/m̃_g^2 - M̃^f/m̃_f^2)e^-μ r/12π r .

This completes the solution of the field equations in terms of gauge-invariant potentials.



 §.§ Metric components

Before we calculate the metric components from the solution for the gauge-invariant potentials, it is convenient to introduce the abbreviations

    ℐ_M = M̃^g + M̃^f/8π(m̃_g^2 + m̃_f^2) ,   ℐ_± = -1/24πμ^2(M̃^g/m̃_g^2±M̃^f/m̃_f^2) ,   𝒟 = -m̃_g^2 - m̃_f^2/m̃_g^2 + m̃_f^2

for a few frequently occurring constants. Further, we use the shorthand notation

    𝒴_μ(r) = e^-μ r/r ,   𝒴_0(r) = 1/r

for the Yukawa and Coulomb potentials. Using these abbreviations, the solution derived above takes the simple form

    I_1^-(2) = 4μ^2ℐ_-𝒴_μ ,   
    I_2^-(2) = 2μ^2ℐ_-𝒴_μ ,   
    I_4^-(2) = 2ℐ_-𝒴_μ ,
    
    I_1^+(2) = -2ℐ_M𝒴_0 + 4μ^2𝒟ℐ_-𝒴_μ ,   
    I_2^+(2) = -2ℐ_M𝒴_0 + 2μ^2𝒟ℐ_-𝒴_μ .

In order to separate these potentials into the potentials for the individual metrics, we further need to fix the pure gauge potential I_4^+(2). A convenient choice, which turns out to be compatible with the standard PPN gauge, is given by

    I_4^+(2) = 2ℐ_+𝒴_μ .

Together with the relations (<ref>) we then obtain the potentials


    I_1^g,f(2)   = -ℐ_M𝒴_0 + 2μ^2(𝒟± 1)ℐ_-𝒴_μ ,
    
    I_2^g,f(2)   = -ℐ_M𝒴_0 + μ^2(𝒟± 1)ℐ_-𝒴_μ ,
    
    I_4^g,f(2)   = (ℐ_+ ±ℐ_-)𝒴_μ .


Now using the relations (<ref>) we obtain the quantities


    ϕ^g,f(2)   = -ℐ_M𝒴_0 + 2μ^2(𝒟± 1)ℐ_-𝒴_μ ,
    ψ^g,f(2)   = -ℐ_M𝒴_0 + μ^2(𝒟± 1)ℐ_-𝒴_μ - 1/3(ℐ_+ ±ℐ_-)𝒴_μ ,
    Ẽ^g,f(2)   = (ℐ_+ ±ℐ_-)𝒴_μ ,


and finally using their definition (<ref>) yields the components of the metric perturbations


    h^(2)_00   = 2ℐ_M𝒴_0 - 4μ^2(𝒟 + 1)ℐ_-𝒴_μ ,
    
    e^(2)_00   = 2ℐ_M𝒴_0 - 4μ^2(𝒟 - 1)ℐ_-𝒴_μ ,
    
    h^(2)_ij   = 2[ℐ_M𝒴_0 - μ^2(𝒟 + 1)ℐ_-𝒴_μ]δ_ij + 2(ℐ_+ + ℐ_-)∂_i∂_j𝒴_μ ,
    
    e^(2)_ij   = 2[ℐ_M𝒴_0 - μ^2(𝒟 - 1)ℐ_-𝒴_μ]δ_ij + 2(ℐ_+ - ℐ_-)∂_i∂_j𝒴_μ .


For later use we now insert the constants (<ref>) and the Yukawa and Coulomb potentials (<ref>). Note that second order derivatives of these potentials contain also delta functions, which must be taken into account for deriving further quantities from the metric perturbations. We have listed the relevant formulas in appendix <ref>. Using these formulas we obtain


    h^(2)_00   = M̃^g + M̃^f/4π(m̃_g^2 + m̃_f^2)r + m̃_f^2M̃^g - m̃_g^2M̃^f/3πm̃_g^2(m̃_g^2 + m̃_f^2)re^-μ r ,
    
    e^(2)_00   = M̃^g + M̃^f/4π(m̃_g^2 + m̃_f^2)r - m̃_f^2M̃^g - m̃_g^2M̃^f/3πm̃_f^2(m̃_g^2 + m̃_f^2)re^-μ r ,
    
    h^(2)_ij   = [M̃^g + M̃^f/4π(m̃_g^2 + m̃_f^2)r - 2m̃_g^2m̃_f^2M̃^f - m̃_g^4M̃^f - 3m̃_f^4M̃^g/18πm̃_g^2m̃_f^2(m̃_g^2 + m̃_f^2)re^-μ r - 2M̃^f/9m̃_f^2μ^2δ(x⃗)]δ_ij
       =+ [μ r(μ r + 3) + 3]M̃^f/6πm̃_f^2μ^2r^5e^-μ r(x_ix_j - 1/3r^2δ_ij) ,
    
    e^(2)_ij   = [M̃^g + M̃^f/4π(m̃_g^2 + m̃_f^2)r - 2m̃_g^2m̃_f^2M̃^g - m̃_f^4M̃^g - 3m̃_g^4M̃^f/18πm̃_g^2m̃_f^2(m̃_g^2 + m̃_f^2)re^-μ r - 2M̃^g/9m̃_g^2μ^2δ(x⃗)]δ_ij
       =+ [μ r(μ r + 3) + 3]M̃^g/6πm̃_g^2μ^2r^5e^-μ r(x_ix_j - 1/3r^2δ_ij) .


Note that the off-diagonal contribution of h^(2)_ij depends only on the mass M̃^f, while the off-diagonal contribution of e^(2)_ij contains only the mass M̃^g. This is a consequence of our gauge choice (<ref>), and the reason for making this choice.



 §.§ PPN parameters

We can now read off the PPN parameters by comparing the solution (<ref>) to the PPN metric ansatz (<ref>). Since we have used rescaled mass units (<ref>), it is convenient to replace the superpotentials (<ref>), which for a point mass source take the form (<ref>), by the correspondingly rescaled superpotentials

    χ̃^g = χ^g = -M^gr = -M̃^gr ,   χ̃^f = cχ^f = -cM^fr = -M̃^fr .

We thus use the modified PPN metric ansatz


    h^(2)_00   = 2α̃^ggM̃^g + α̃^gfM̃^f/r ,    
    h^(2)_ij   = 2γ̃^ggM̃^g + γ̃^gfM̃^f/rδ_ij + 2θ̃^ggM̃^g + θ̃^gfM̃^f/r^3x_ix_j ,
    
    e^(2)_00   = 2α̃^fgM̃^g + α̃^ffM̃^f/r ,    
    e^(2)_ij   = 2γ̃^fgM̃^g + γ̃^ffM̃^f/rδ_ij + 2θ̃^fgM̃^g + θ̃^ffM̃^f/r^3x_ix_j .


Note that the observable parameters α^gg = α̃^gg, γ^gg = γ̃^gg and θ^gg = θ̃^gg, which govern the gravitational interaction within the visible matter sector, are unaffected by this rescaling, and that only the PPN parameters involving the dark sector receive constant factors. We then read off the PPN parameters


    α̃^gg   = 3m̃_g^2 + 4m̃_f^2e^-μ r/24πm̃_g^2(m̃_f^2 + m̃_g^2) ,    α̃^gf   = 3 - 4e^-μ r/24π(m̃_f^2 + m̃_g^2) ,
    α̃^ff   = 3m̃_f^2 + 4m̃_g^2e^-μ r/24πm̃_f^2(m̃_f^2 + m̃_g^2) ,    α̃^fg   = 3 - 4e^-μ r/24π(m̃_f^2 + m̃_g^2) ,
    γ̃^gg   = 3m̃_g^2 + 2m̃_f^2e^-μ r/24πm̃_g^2(m̃_f^2 + m̃_g^2) ,    γ̃^gf   = 9m̃_f^2 + 2(m̃_g^2 - 2m̃_f^2)e^-μ r/72πm̃_f^2(m̃_f^2 + m̃_g^2) - μ r(μ r + 3) + 3/36πm̃_f^2μ^2r^2e^-μ r ,
    γ̃^ff   = 3m̃_f^2 + 4m̃_g^2e^-μ r/24πm̃_f^2(m̃_f^2 + m̃_g^2) ,    γ̃^fg   = 9m̃_g^2 + 2(m̃_f^2 - 2m̃_g^2)e^-μ r/72πm̃_g^2(m̃_f^2 + m̃_g^2) - μ r(μ r + 3) + 3/36πm̃_g^2μ^2r^2e^-μ r ,
    θ̃^gg   = 0 ,    θ̃^gf   = μ r(μ r + 3) + 3/12πm̃_f^2μ^2r^2e^-μ r ,
    θ̃^ff   = 0 ,    θ̃^fg   = μ r(μ r + 3) + 3/12πm̃_g^2μ^2r^2e^-μ r .


We find that the gauge condition θ̃^gg = θ̃^ff = 0, which we have introduced in section <ref>, is satisfied, due to our choice (<ref>). From these parameters we can in particular derive the observable quantities

    G_eff = α̃^gg = 3m̃_g^2 + 4m̃_f^2e^-μ r/24πm̃_g^2(m̃_f^2 + m̃_g^2) ,   γ = γ̃^gg/α̃^gg = 3m̃_g^2 + 2m̃_f^2e^-μ r/3m̃_g^2 + 4m̃_f^2e^-μ r ,

which are the effective Newtonian constant and the usual PPN parameter γ. Both quantities depend on the distance r between the mass source and the location where the gravitational field is probed, in contrast to general relativity, where both quantities are constant. It is further remarkable that γ depends only on the ratio m̃_f/m̃_g of the two Planck masses and the graviton mass μ, and that this result essentially resembles the observable parameters of scalar-tensor theory with a general potential <cit.>, or the more general Horndeski class of theories <cit.>, which depend on the Brans-Dicke parameter ω and the scalar field mass.



 §.§ Limiting cases

We finally discuss a few interesting limiting cases for the mass parameters m̃_g,f and μ and their consequences for the PPN parameters. These are in particular:



  * It is well known that in the limit m̃_f → 0, while keeping the parameters m and β_k in the interaction potential fixed, one obtains the general relativity limit for the visible sector <cit.>. Note that in this limit we also have μ→∞. The PPN parameters (<ref>) then take the form


    α̃^gg = α̃^gf = α̃^fg = α̃^ff   = 1/8πm̃_g^2 ,
    γ̃^gg = γ̃^gf = γ̃^fg = γ̃^ff   = 1/8πm̃_g^2 ,
    θ̃^gg = θ̃^gf = θ̃^fg = θ̃^ff   = 0 ,


while the observable parameters (<ref>) are given by

    G_eff = 1/8πm̃_g^2 ,   γ = 1 ,

as usual in general relativity.


  * For equal Planck mass parameters m̃_g = m̃_f one obtains the PPN parameters


    α̃^gg = α̃^ff   = 3 + 4e^-μ r/48πm̃_g^2 ,    α̃^gf = α̃^fg   = 3 - 4e^-μ r/48πm̃_g^2 ,
    γ̃^gg = γ̃^ff   = 3 + 2e^-μ r/48πm̃_g^2 ,    γ̃^gf = γ̃^fg   = 3 - 2e^-μ r/48πm̃_g^2 - (μ r + 1)e^-μ r/12πm̃_g^2μ^2r^2 ,
    θ̃^gg = θ̃^ff   = 0 ,    θ̃^gf = θ̃^fg   = [μ r(μ r + 3) + 3]e^-μ r/12πm̃_g^2μ^2r^2


and the observable parameters

    G_eff = 3 + 4e^-μ r/48πm̃_g^2 ,   γ = 3 + 2e^-μ r/3 + 4e^-μ r .

We remark that this result is similar to the PPN parameter γ in higher-order gravity, except for an additional scalar contribution and a different sign due to the massive graviton being a ghost in the latter class of theories <cit.>. Note that the effective Planck mass for the visible sector,

    m_Pl^2 = lim_r →∞1/8π G_eff ,

is given by m_Pl^2 = 2m̃_g^2.


  * In the limit μ→∞ of a highly massive graviton we find the PPN parameters


    α̃^gg = α̃^gf = α̃^fg = α̃^ff   = 1/8π(m̃_g^2 + m̃_f^2) ,
    γ̃^gg = γ̃^gf = γ̃^fg = γ̃^ff   = 1/8π(m̃_g^2 + m̃_f^2) ,
    θ̃^gg = θ̃^gf = θ̃^fg = θ̃^ff   = 0 ,


from which follow the observable parameters

    G_eff = 1/8π(m̃_g^2 + m̃_f^2) ,   γ = 1 .

In this case the effective Planck mass (<ref>) turns out to be m_Pl^2 = m̃_g^2 + m̃_f^2.


This concludes our discussion of the post-Newtonian limit of ghost-free bimetric gravity for a static point mass. The PPN parameters we have obtained now allow us to discuss observable effects, and in particular the deflection of light by both dark and visible matter. This will be done in the following section.



§ CONFRONTATION WITH OBSERVATIONS

In the previous section we obtained both a general result and a number of limiting cases for the effective gravitational constant and the PPN parameter γ, as well as additional PPN parameters which govern effects involving a second, dark type of matter. We can now compare our results with observations, in particular of the deflection of light. We will restrict ourselves to visible matter in section <ref> and derive bounds on the parameters of ghost-free massive bimetric gravity from solar system experiments. In section <ref> we will discuss the deflection of visible light by dark matter and its consistency with observations of lensing effects by galaxies. We will further speculate on a possible explanation for the lensing effects observed in the vicinity of galactic mergers, in particular Abell 520 and Abell 3827.



 §.§ Solar system consistency

We have remarked in section <ref> that our result (<ref>) for the effective Newtonian constant G_eff and the PPN parameter γ has essentially the same form as the corresponding result for scalar-tensor gravity with a general potential <cit.>, or the more general Horndeski class of theories <cit.>. Hence the experimental constraints on the parameters of these theories derived from measurements of γ can directly be translated to constraints on the parameters of ghost-free massive bimetric gravity, and in particular to the ratio m̃_f/m̃_g of the Planck masses and the graviton mass μ. An important obstacle that must be taken into account is the fact that γ is not constant, but depends exponentially on the distance r between the gravitating mass source and the observer. This restricts the possible experimental tests of γ to those for which such an interaction distance can be defined. The most precise observation of γ which satisfies this condition is the measurement of the Shapiro time delay of radio signals between Earth and the Cassini spacecraft on its way to Saturn, from which a value γ - 1 = (2.1 ± 2.3) · 10^-5 was obtained <cit.>. These were passing by the sun at a distance of 1.6 solar radii, so that we define the interaction distance r_0 ≈ 7.44 · 10^-3AU. Following the same procedure as detailed in <cit.>, we find that the area of the parameter space shown in figure <ref> is excluded at 2σ confidence level. Note, however, that the assumption of a constant interaction distance for this experiment is only an approximation, and that more accurate results are obtained from a thorough treatment of light propagation in the solar gravitational field <cit.>.





 §.§ Light deflection by dark matter

The full set (<ref>) of PPN parameters, which we derived in section <ref>, allows us to discuss also the gravitational interaction of dark matter. We first consider the parameter α̃^gf, which can be interpreted as an effective Newtonian constant for the gravitational influence of dark matter Φ^f on visible matter Φ^g. For short distances, μ r < ln(4/3), we see that α̃^gf becomes negative, so that the gravitational interaction between dark and visible matter becomes repulsive; however, taking into account the bounds shown in figure <ref>, we see that this is possible only on scales significantly smaller than the solar system, and hence does not play any role for the observed dark matter concentrations. On the scales of galaxies or even galactic clusters we can safely assume μ r ≫ 1, and thus use the PPN parameters obtained in the limit μ→∞ in section <ref>. In this limit, the gravitational effects on both test masses and light become indistinguishable between visible and dark matter sources. In particular, it follows that the deflection of visible light by dark matter is likewise governed by a PPN parameter

    γ̅ = γ̃^gf/α̃^gf→ 1

in the limit of large scales. This agrees with observations of the deflection of light by galaxies, which contain significant amounts of dark matter in addition to the visible mass <cit.>.

Our result plays an important role in particular for the observed light deflection by galactic mergers, such as most prominently the so-called “Bullet Cluster” 1E0657-558 <cit.> or more recently MACS J0025.4-1222 <cit.>. Measurements of the mass distribution in these and other mergers using weak lensing together with x-ray imaging show that the gas component of the merger, which is heated by the collision and which constitutes the major amount of visible matter, is not at the same location as the dominant gravitating matter contribution, and that the motion of the latter is largely unaffected by the collision. This leads to the conclusion that their dark matter content is non-interacting, so that the dark matter components of the colliding objects pass through each other <cit.>. However, observations of the so-called “Train Wreck Cluster” Abell 520 <cit.> or Abell 3827 <cit.> show a more differentiated picture. While also Abell 520 shows evidence for dark matter components which have passed through each other unaffectedly, one has further identified another dark mass concentration in the central region, which is difficult to explain if dark matter is non-interacting. Similar stress on the non-interacting dark matter model is put by an observed separation between stellar and dark matter in Abell 3827. A possible explanation for these observations is to assume that dark matter also possesses a component which interacts non-gravitationally <cit.>.

The bimetric class of theories we studied in this article allows for an interesting tentative model for the aforementioned observations, which hint towards the existence of both interacting and non-interacting dark matter components. Invoking the interpretation of the matter sector Φ^f as dark matter, as suggested in <cit.>, and further assuming that Φ^f contains an interacting component, would suggest that the central dark matter concentration in Abell 520 and the separated dark matter concentration in Abell 3827 result from a collision of these interacting components, while any dark matter constituted by massive gravitons, as suggested in <cit.>, would pass the merger unaffectedly, and could thus account for the dark matter concentrations away from the center of Abell 520 or the unaffected dark matter halos in Abell 3827. Future extensions of our work presented here will be necessary in order to quantitatively assert the viability of such models.



§ CONCLUSION

We have considered the post-Newtonian limit of ghost-free massive bimetric gravity with two mutually non-interacting matter sectors. From the assumption that the vacuum field equations are solved by two flat metrics proportional to the Minkowski metric, we have derived restrictions on the parameters in the action. For this restricted class of theories we have derived the field equations up to the second velocity order by making use of a suitable extension of the PPN formalism to multiple metrics. We have solved these equations for a point-like mass source using a gauge-invariant differential decomposition of the metric perturbations. From this solution we have read off the effective gravitational constant G_eff and the PPN parameter γ for the visible matter sector. By comparing our result to the observed value determined by the Cassini tracking experiment we have derived combined bounds on two parameters of the theory, namely on the mass of the massive graviton and on the ratio of the Planck masses occurring in the bimetric action.

We have further discussed the interpretation of the additional matter sector as a possible constituent of dark matter. From our experimental bounds we then concluded that on scales significantly larger than the solar system, and hence in particular on the observationally relevant scales of galaxies and clusters, the gravitational effects caused by visible and dark matter become indistinguishable from each other. It thus follows that dark matter should deflect light in the same way as visible matter does, in agreement with measurements of the PPN parameter γ through the lensing effect of galaxies, which contain a significant dark matter component. Another possible experimental test of this result could be performed by searching for possible (non-)correlations between the ratio of dark to visible matter of a galaxy and its light deflection. Such an analysis would be most effective with data of higher precision than available to date <cit.>.

On a more speculative note, we have considered that besides the second matter sector also massive gravitons could contribute to the observed dark matter content of the universe. The assumption that the former contains non-gravitational self-interactions, while the latter interacts only gravitationally, then provides a tentative explanation for the observed separation of apparently different dark matter components in galactic mergers such as Abell 520 and Abell 3827. The question arises whether such different dark matter constituents could be distinguished also in other processes besides galactic mergers, for example, by their light deflection properties. An extension of our work presented here to the light deflection caused by massive graviton concentrations might answer this question.

There are also other possibilities to further extend the theoretical analysis we presented in this article. While we have studied only linear perturbations of flat vacuum solutions, considering also the quadratic perturbation order would allow us to calculate the PPN parameter β, and thus open the possibility for additional tests using solar system observations. This would ultimately lead to a full generalization of the formalism developed in <cit.> to massive gravity theories. Further, one may also include cosmological corrections to the PPN formalism along the lines of <cit.>, and thus relax the condition of a flat background. Finally, one may consider more general theories with N > 2 metric tensors and a corresponding number of matter sectors <cit.>, or involving an effective metric <cit.>, both of which allow for ghost-free matter coupling prescriptions <cit.>. We intend to study these generalizations in future research.


The author is happy to thank the members of the Laboratory of Theoretical Physics at the University of Tartu for fruitful discussions. He gratefully acknowledges the full financial support of the Estonian Research Council through the Startup Research Grant PUT790 and the European Regional Development Fund through the Center of Excellence TK133 “The Dark Side of the Universe”.




§ LINEARIZATION OF THE POTENTIAL

In this appendix we show how to obtain the linearized potentials, which enter the gravitational field equations at the zeroth and second velocity order as shown in sections <ref> and <ref>. The starting point for our derivation is a linear perturbation ansatz for the metrics, which we write in the form

    g_μν = η_μν + h_μν ,  
    f_μν = c^2(η_μν + e_μν) .

Up to the linear perturbation order, we can then write their inverses as

    g^μν = η^μν - η^μρη^νσh_ρσ + 𝒪(h^2) ,   
    f^μν = 1/c^2(η^μν - η^μρη^νσe_ρσ) + 𝒪(e^2) .

For their product we find

    g^μρf_ρν = c^2(δ^μ_ν - D^μ_ν) + 𝒪({h,e}^2) ,

where we introduced the perturbation tensor

    D^μ_ν = η^μρ(h_ρν - e_ρν) .

Since the matrix g^μρf_ρν is given as a perturbation of the Kronecker symbol δ^μ_ν, we can find its square root A^μ_ν as defined in (<ref>) using a series expansion analogously to the well-known Taylor series

    √(1 + x) = 1 + x/2 + 𝒪(x^2) .

This series expansion yields

    A^μ_ν = c(δ^μ_ν - 1/2D^μ_ν) + 𝒪({h,e}^2) .

For later use we also need to expand powers of A into linear perturbations. These are given by

    (A^k)^μ_ν = c^k(δ^μ_ν - k/2D^μ_ν) + 𝒪({h,e}^2) .

The matrix invariants e_k(A) defined by (<ref>) then take the form


    e_0(A)    = 1 ,
    
    e_1(A)    = A^μ_μ = c(4 - 1/2D^μ_μ) + 𝒪({h,e}^2) ,
    
    e_2(A)    = 1/2(A^μ_μA^ν_ν - A^μ_νA^ν_μ) = c^2(6 - 3/2D^μ_μ) + 𝒪({h,e}^2) ,
    
    e_3(A)    = 1/6(A^μ_μA^ν_νA^ρ_ρ - 3A^μ_νA^ν_μA^ρ_ρ + 2A^μ_νA^ν_ρA^ρ_μ) = c^3(4 - 3/2D^μ_μ) + 𝒪({h,e}^2) ,
    
    e_4(A)    = 1/24(A^μ_μA^ν_νA^ρ_ρA^σ_σ - 6A^μ_μA^ν_νA^ρ_σA^σ_ρ + 3A^μ_νA^ν_μA^ρ_σA^σ_ρ + 8A^μ_νA^ν_ρA^ρ_μA^σ_σ
       =- 6A^μ_νA^ν_ρA^ρ_σA^σ_μ) = c^4(1 - 1/2D^μ_μ) + 𝒪({h,e}^2) .


For the matrices Y_n defined via (<ref>) we then find the expressions


    Y_0^μ_ν(A)    = δ^μ_ν + 𝒪({h,e}^2) ,
    
    Y_1^μ_ν(A)    = c(-3δ^μ_ν - 1/2D^μ_ν + 1/2D^ρ_ρδ^μ_ν) + 𝒪({h,e}^2) ,
    
    Y_2^μ_ν(A)    = c^2(3δ^μ_ν + D^μ_ν - D^ρ_ρδ^μ_ν) + 𝒪({h,e}^2) ,
    
    Y_3^μ_ν(A)    = c^3(-δ^μ_ν - 1/2D^μ_ν + 1/2D^ρ_ρδ^μ_ν) + 𝒪({h,e}^2) .


In order to obtain the corresponding expressions for A^-1 = √(f^-1g) instead of A, one simply replaces D by -D and c by c^-1. We can now calculate the potentials (<ref>). Using renormalized parameters β̃_k = c^kβ_k we obtain


    V^g_μν   = [(β̃_0 + 3β̃_1 + 3β̃_2 + β̃_3)η_μν - (1/2β̃_1 + β̃_2 + 1/2β̃_3)η_μνη^ρσ(h_ρσ - e_ρσ)
       =+ (β̃_0 + 7/2β̃_1 + 4β̃_2 + 3/2β̃_3)h_μν - (1/2β̃_1 + β̃_2 + 1/2β̃_3)e_μν] + 𝒪({h,e}^2) ,
    
    V^f_μν   = 1/c^2[(β̃_1 + 3β̃_2 + 3β̃_3 + β̃_4)η_μν + (1/2β̃_1 + β̃_2 + 1/2β̃_3)η_μνη^ρσ(h_ρσ - e_ρσ)
       =+ (3/2β̃_1 + 4β̃_2 + 7/2β̃_3 + β̃_4)e_μν - (1/2β̃_1 + β̃_2 + 1/2β̃_3)h_μν] + 𝒪({h,e}^2) .


Finally, we calculate the trace-reversed potentials (<ref>), which are given by


    V̅^g_μν   = [-(β̃_0 + 3β̃_1 + 3β̃_2 + β̃_3)η_μν + (1/4β̃_1 + 1/2β̃_2 + 1/4β̃_3)η_μνη^ρσ(h_ρσ - e_ρσ)
       =- (β̃_0 + 5/2β̃_1 + 2β̃_2 + 1/2β̃_3)h_μν - (1/2β̃_1 + β̃_2 + 1/2β̃_3)e_μν] + 𝒪({h,e}^2) ,
    V̅^f_μν   = 1/c^2[-(β̃_1 + 3β̃_2 + 3β̃_3 + β̃_4)η_μν - (1/4β̃_1 + 1/2β̃_2 + 1/4β̃_3)η_μνη^ρσ(h_ρσ - e_ρσ)
       =- (1/2β̃_1 + 2β̃_2 + 5/2β̃_3 + β̃_4)e_μν - (1/2β̃_1 + β̃_2 + 1/2β̃_3)h_μν] + 𝒪({h,e}^2) .


These expressions can now be used in the post-Newtonian field equations at the zeroth velocity order in section <ref> and at the second velocity order in section <ref>.



§ DERIVATIVES OF THE YUKAWA POTENTIAL

During our calculation we have frequently encountered (mostly second order) derivatives of the Yukawa potential, for which we introduced the shorthand notation

    𝒴_k(r) = e^-kr/r .

Taking into account the singularity at the origin, its second derivatives are given by

    ∂_i∂_j𝒴_k = {[kr(kr + 3) + 3]x_ix_j/r^5 - (kr + 1)δ_ij/r^3}e^-kr - 4π/3δ_ijδ(x⃗) ,

which is a straightforward generalization of the well-known formula for the Coulomb potential <cit.>. Taking the trace yields the standard formula

    𝒴_k = k^2e^-kr/r - 4πδ(x⃗) .

These formulas cover all expressions which appear in the final result for the Ricci tensor and the interaction potential. Note that during intermediate steps also fourth order derivatives of the Yukawa potential occur in derivatives of the metric perturbations. For completeness we also list the corresponding expressions. From the formula given above immediately follows

    ∂_i∂_j𝒴_k = k^2{[kr(kr + 3) + 3]x_ix_j/r^5 - (kr + 1)δ_ij/r^3}e^-kr - 4π k^2/3δ_ijδ(x⃗) - 4π∂_i∂_jδ(x⃗)

and thus

    𝒴_k = k^4e^-kr/r - 4π k^2δ(x⃗) - 4πδ(x⃗) .

These are all terms which occur during our calculation.



§ CHECKING THE FIELD EQUATIONS IN COMPONENTS

Since we have used a rather technical transformation of the field equations to gauge invariant potentials in section <ref> and the corresponding inverse transformation of their solution to metric components in section <ref>, it is appropriate to check the obtained result also using the field equations in their original component form as shown in section <ref>. While this is rather cumbersome using the explicit expressions (<ref>) and requires careful tracking of singular contributions from higher derivatives of Coulomb and Yukawa potentials, it becomes considerably simpler by using the abbreviations (<ref>), starting from the expressions (<ref>) and finally evaluating higher derivatives using the formulas shown in appendix <ref>.

From the expressions (<ref>) one easily reads off the traces


    h^(2)_ii   = 2[3ℐ_M𝒴_0 - 3μ^2(𝒟 + 1)ℐ_-𝒴_μ + (ℐ_+ + ℐ_-)𝒴_μ] ,
    
    e^(2)_ii   = 2[3ℐ_M𝒴_0 - 3μ^2(𝒟 - 1)ℐ_-𝒴_μ + (ℐ_+ - ℐ_-)𝒴_μ]


of the spatial components of the metric perturbations. Using the formulas (<ref>) for the potential at the second velocity order we then obtain


    V^g(2)_00 = -c^2V^f(2)_00   = -β̃ℐ_-(3μ^2𝒴_μ + 𝒴_μ) ,
    
    V^g(2)_ij = -c^2V^f(2)_ij   = -β̃ℐ_-(3μ^2𝒴_μδ_ij - 𝒴_μδ_ij - 2∂_i∂_j𝒴_μ) .


Further, we need to evaluate second order derivatives of the metric, which read


    h^(2)_00,ij   = 2ℐ_M∂_i∂_j𝒴_0 - 4μ^2(𝒟 + 1)ℐ_-∂_i∂_j𝒴_μ ,
    
    e^(2)_00,ij   = 2ℐ_M∂_i∂_j𝒴_0 - 4μ^2(𝒟 - 1)ℐ_-∂_i∂_j𝒴_μ ,
     h^(2)_00   = 2ℐ_M𝒴_0 - 4μ^2(𝒟 + 1)ℐ_-𝒴_μ ,
     e^(2)_00   = 2ℐ_M𝒴_0 - 4μ^2(𝒟 - 1)ℐ_-𝒴_μ ,
    
    h^(2)_kk,ij   = 2[3ℐ_M∂_i∂_j𝒴_0 - 3μ^2(𝒟 + 1)ℐ_-∂_i∂_j𝒴_μ + (ℐ_+ + ℐ_-)∂_i∂_j𝒴_μ] ,
    
    e^(2)_kk,ij   = 2[3ℐ_M∂_i∂_j𝒴_0 - 3μ^2(𝒟 - 1)ℐ_-∂_i∂_j𝒴_μ + (ℐ_+ - ℐ_-)∂_i∂_j𝒴_μ] ,
     h^(2)_ij   = 2[ℐ_M𝒴_0 - μ^2(𝒟 + 1)ℐ_-𝒴_μ]δ_ij + 2(ℐ_+ + ℐ_-)∂_i∂_j𝒴_μ ,
     e^(2)_ij   = 2[ℐ_M𝒴_0 - μ^2(𝒟 - 1)ℐ_-𝒴_μ]δ_ij + 2(ℐ_+ - ℐ_-)∂_i∂_j𝒴_μ ,
     h^(2)_ii   = 2[3ℐ_M𝒴_0 - 3μ^2(𝒟 + 1)ℐ_-𝒴_μ + (ℐ_+ + ℐ_-)𝒴_μ] ,
     e^(2)_ii   = 2[3ℐ_M𝒴_0 - 3μ^2(𝒟 - 1)ℐ_-𝒴_μ + (ℐ_+ - ℐ_-)𝒴_μ] ,
    
    h^(2)_ik,jk   = 2[ℐ_M∂_i∂_j𝒴_0 - μ^2(𝒟 + 1)ℐ_-∂_i∂_j𝒴_μ + (ℐ_+ + ℐ_-)∂_i∂_j𝒴_μ] ,
    
    e^(2)_ik,jk   = 2[ℐ_M∂_i∂_j𝒴_0 - μ^2(𝒟 - 1)ℐ_-∂_i∂_j𝒴_μ + (ℐ_+ - ℐ_-)∂_i∂_j𝒴_μ] .


Inserting these expressions into the formulas (<ref>) for the Ricci tensor then yields the components


    R^g,f(2)_00   = -ℐ_M𝒴_0 + 2μ^2(𝒟± 1)ℐ_-𝒴_μ ,
    
    R^g,f(2)_ij   = -ℐ_M𝒴_0δ_ij + μ^2(𝒟± 1)ℐ_-(𝒴_μδ_ij - ∂_i∂_j𝒴_μ) .


Inserting the expressions (<ref>) and (<ref>) into the second order field equations (<ref>), applying the definitions (<ref>) and using the relations for the Coulomb and Yukawa potentials listed in appendix <ref> finally yields


    m̃_g^2R^g(2)_00 + m^4V̅^g(2)_00   = -m̃_g^2(M̃^g + M̃^f)𝒴_0 - (m̃_g^2M̃^f - m̃_f^2M̃^g)(𝒴_μ - μ^2𝒴_μ)/8π(m̃_g^2 + m̃_f^2)
       = M̃^g/2δ(x⃗) = T̅^g(2)_00 ,
    m̃_f^2/c^2R^f(2)_00 + m^4V̅^f(2)_00   = -m̃_f^2(M̃^g + M̃^f)𝒴_0 + (m̃_g^2M̃^f - m̃_f^2M̃^g)(𝒴_μ - μ^2𝒴_μ)/8π c^2(m̃_g^2 + m̃_f^2)
       = M̃^f/2δ(x⃗) = T̅^f(2)_00 ,
    m̃_g^2R^g(2)_ij + m^4V̅^g(2)_ij   = -m̃_g^2(M̃^g + M̃^f)𝒴_0 - (m̃_g^2M̃^f - m̃_f^2M̃^g)(𝒴_μ - μ^2𝒴_μ)/8π(m̃_g^2 + m̃_f^2)δ_ij
       = M̃^g/2δ(x⃗)δ_ij = T̅^g(2)_ij ,
    m̃_f^2/c^2R^f(2)_ij + m^4V̅^f(2)_ij   = -m̃_f^2(M̃^g + M̃^f)𝒴_0 + (m̃_g^2M̃^f - m̃_f^2M̃^g)(𝒴_μ - μ^2𝒴_μ)/8π c^2(m̃_g^2 + m̃_f^2)δ_ij
       = M̃^f/2δ(x⃗)δ_ij = T̅^f(2)_ij .


This shows that the field equations are indeed satisfied.

apsrev4-1

