Combining cosmological and local bounds on bimetric theory

Ghost-free bimetric theory describes two nonlinearly interacting spin-2 fields, one massive and one massless, thus extending general relativity. We confront bimetric theory with observations of Supernovae type 1a, Baryon Acoustic Oscillations and the Cosmic Microwave Background in a statistical analysis, utilising the recently proposed physical parametrisation. This directly constrains the physical parameters of the theory, such as the mass of the spin-2 field and its coupling to matter. We find that all models under consideration are in agreement with the data. Next, we compare these results to bounds from local tests of gravity. Our analysis reveals that all two- and three-parameter models are observationally consistent with both cosmological and local tests of gravity. The minimal bimetric model (only $\beta_1$) is ruled out by our combined analysis.


Introduction
In our current era of precision cosmology, cosmological tests of gravity with remarkable accuracy become available. On the one hand, the direct detection of gravitational waves by the LIGO/Virgo collaboration [1] has confirmed some of the most fundamental predictions of general relativitygravitational waves [2]. Moreover, the joint detection of gravitational waves with their electromagnetic counterpart [3] show that gravitational waves indeed travel at the speed of light to large precision [4]. On the other hand, measurements of the cosmic microwave background [5,6] and local observables [7][8][9][10] challenge the standard model of cosmology. In fact, the discrepancy between the inferred value for the Hubble parameter today (H 0 ) from early and late-time measurements represents a substantial tension [11]. From a theoretical perspective, the cosmological constant problems remain unsolved [12,13]. In addition, concerns about the quantum consistency of a positive cosmological constant have been raised, in the context of quantum breaking [14][15][16][17] and the swampland 1 [19][20][21][22]. Therefore, it is worth investigating theories beyond the standard model of cosmology, which modify the gravitational theory or replace the cosmological constant by some dark energy [23]. To modify gravity, one has to break one of the axioms of Lovelock's theorem [24]. A particularly simple possibility is to add new degrees of freedom to the gravitational sector. The resulting scalar-tensor [25,26] and vector-tensor [27,28] theories are already tightly constrained by gravitational wave observations and forced into their simplest forms [29][30][31][32]. Also massive gravity [33,34] is tightly constrained observationally [35]. Contrarily, bimetric theory [36] remains mostly unconstrained by the determination of the propagation speed of gravitational waves [37,38]. In the present paper, we will further investigate the observational viability of the latter theory.
Historically, the question of whether the graviton can have a finite mass goes back to Fierz and Pauli, who presented a linear theory of a massive graviton in [39]. However, this theory failed even basic solar-system tests of gravity due to the van Dam-Veltman-Zakharov (vDVZ) discontinuity [40,41]. Shortly after, Vainshtein argued that nonlinear terms might cure the discontinuity and render a nonlinear theory of a massive graviton observationally viable [42]. However, Boulware and Deser argued that a fully nonlinear theory of a massive graviton will inevitably introduce an additional degree of freedom (d.o.f.) with negative norm 2 . Rather recently, de Rham, Gabadadze and Tolley identified a loop-hole in the argument, which lead to the construction of a fully nonlinear and ghostfree theory of massive gravity [33,34,46] 3 . The construction of a nonlinear theory for a massive graviton requires the introduction of a reference metric, which is kept fixed. A natural extension is to promote this fixed reference metric to be dynamical, which leads to ghost-free bimetric theory [36].
While massive gravity describes a single massive graviton, bimetric theory contains a massless and a massive spin-2 field [36]. This yields distinct phenomenological features. The mass m FP of a single massive graviton is tightly constrained. While gravitational wave observations provide the upper bound of m FP 10 −22 eV [32], the most stringent constraints come from the solar system with m FP 10 −33 eV [35]. Moreover, massive gravity does not give rise to viable homogenous and isotropic cosmological solutions [53][54][55][56]. Reviews on theoretical and phenomenological aspects of massive gravity can be found in [35,57,58].
To consistently combine the various observational and theoretical constraints on the parameter space, a new parametrisation of bimetric solutions was developed in [105]. The idea is to formulate solutions in terms of the following physical parameters: mass of the spin-2 field, its coupling constant to ordinary matter and the effective asymptotic cosmological constant. The procedure automatically yields theoretical constraints that ensure a viable cosmic expansion history, i.e. real-valued, nonsingular and devoid of the Higuchi ghost.
In this paper, we extend the analysis of [105] to test bimetric cosmology against data from Supernovae type 1a (SN1a), Baryon Acoustic Oscillations (BAOs) and the Cosmic Microwave Background (CMB). While bimetric theory has been tested against these data sets beforehand [60,64,67,68,106,107], here we directly constrain the physical parameters. We restrict our analysis to all theoretically viable bimetric models with up to three free interaction parameters 4 . In addition, we allow for non-zero spatial curvature in the statistical analysis.
Next, we compare the cosmological bounds to bounds from local tests of gravity. We use tests of the Newtonian gravitational potential from laboratory to planetary scale as summarised in [108][109][110][111]. The Yukawa parametrisation, where the gravitational potential is written as a superposition of the usual 1/r-term and a Yukawa-term, is appropriate for our purposes. We use the existing bounds to constrain the physical parameters of bimetric theory. We add further bounds from tests of the scalar curvature, i.e. measurements of the deflection and time delay of light [91,112,113].
The aforementioned frameworks are appropriate only in regimes without Vainshtein screening. In screened spacetime regions these bounds are not trustworthy, because deviations from general relativity (GR) are suppressed. However, only a subregion of the full bimetric parameter space supports Vainshtein screening. In this paper, we relate the bounds of [90,93], which ensure a working Vainshtein mechanism, to the physical parameters and identify the region of the physical parameter space that supports screening. We will see that the Vainshtein mechanism does not affect the observational bounds from local tests of gravity, such that these are directly applicable. Summarising, we consistently combine for the first time these observational and theoretical constraints from cosmology and local tests of gravity with each other. This paper is organised as follows. We start with a brief summary of bimetric theory, the cosmological solutions and the physical parameters in section 2. We proceed with the statistical analysis using cosmological data in section 3. In section 4, we discuss the bounds from local tests of gravity and the Vainshtein screening, which we confront with the cosmological constraints. We summarise our analysis and conclude in section 5.

Review of bimetric theory
In this section we briefly discuss bimetric theory, its mass spectrum, cosmological solutions and the physical parametrisation. For a thorough introduction to the field we refer to [114].
The construction of a nonlinear mass term involves two metric tensors, which we call g µν and f µν . The action of ghost-free bimetric theory for these metrics reads [33,34,36,46,115] where R g and R f are the Ricci scalars of g µν and f µν , resp. The parameter m g is the Planck mass of g µν and α is the Planck mass of f µν normalised to m g . The ghost-free bimetric potential U is defined in terms of the square-root matrix S = g −1 f and given by 5 Here, β n are real parameters with mass dimension two in our parametrisation. The functions e n are the elementary symmetric polynomials [115]. Matter fields are collectively denoted by Φ, which minimally couple 6 to the metric g µν via the matter Lagrangian L m . Hence, g µν is the physical metric that defines the geometry in which the matter fields propagate. Varying the action with respect to g µν and f µν yields the modified Einstein equations where G g,f µν are the Einstein tensors of g µν and f µν , resp. The terms coming from the variation of the bimetric potential are given by where the matrices Y λ (n)ν are given in, e.g. [115]. Since matter couples to the physical metric g µν , the matter stress-energy tensor is given by If the matter sector is invariant under diffeomorphisms, the stress-energy tensor of matter is conserved, where ∇ µ is the covariant derivative compatible with g µν . The Bianchi identity, ∇ µ G g µν = 0, then yields the so-called Bianchi constraint, The Bianchi identity in the f µν -sector yields another Bianchi constraint, which however coincides with eq. (2.7) by diffeomorphism invariance.

Vacuum solutions and mass eigenstates
Before moving to cosmology, we briefly review an important class of solutions in bimetric theory [36,123]. Let both metric tensors be proportional as f µν = c 2 g µν with conformal factor c. This ansatz is a solution to the modified Einstein equations only in vacuum, i.e. for T µν = 0. The Bianchi constraint (2.7) implies that c is a constant. Upon this ansatz, the modified Einstein equations (2.3) simplify to where the effective cosmological constants are given in terms of the bimetric parameters as (2.9) Subtracting both equations in (2.8) and using that both metrics are conformally related yields Λ g = Λ f ≡ Λ, or in terms of the bimetric parameters This is a polynomial in c of degree 4. Hence, there are up to four real-valued solutions for c in terms of the bimetric parameters α and β n , each corresponds to a vacuum solution. Bimetric theory has a well-defined mass spectrum around proportional solutions. Let us perturb the metrics around the backgroundḡ µν as where δg µν m g and δf µν αm g are small fluctuations around the proportional background f µν = c 2ḡ µν . The mass eigenstates are linear superpositions of the metric fluctuations. The spectrum contains a massless mode, δG µν , and a massive mode, δM µν , with Fierz-Pauli mass 7 (2.12) We can express the metric fluctuations in terms of the mass eigenstates as Here, we definedᾱ = αc that parametrises the mixing of the massless and the massive mode in the metric fluctuations. From here we can already identify the GR and massive gravity limit of bimetric theory [80,92]. The GR-limit isᾱ → 0 because then the fluctuation of g µν is aligned with the massless mode while the fluctuation of f µν is aligned with the massive mode. In the opposite limitᾱ → ∞ bimetric theory approaches massive gravity as then δg µν is aligned with the massive mode while δf µν coincides with the massless mode.

FLRW solutions
Assuming homogeneity and isotropy according to the cosmological principle, both metric can be cast into the Friedmann-Lemaître-Robertson-Walker (FLRW) form [59][60][61], (2.14) where a and b are the scale factors of g µν and f µν , respectively, while X is the lapse of f µν . These metric functions depend on timte t only. Further, k is the spatial curvature, which must be common to both metrics [124]. In this coordinate system, k < 0, k = 0, k > 0 describes a closed, flat, or open universe, resp. Let us define the Hubble rates and the scale factor ratio as where dot denotes derivative with respect to cosmic time t.
On the bidiagonal FLRW ansatz (2.14), the Bianchi constraint (2.7) simplifies to This equation has two branches of solution. We use the so-called dynamical branch solution with X =ḃ/ȧ for the remainder of this paper 8 . Then the Hubble rates are related as H = yH f . We take the matter source to be a perfect fluid with stress-energy tensor with energy density ρ, pressure p and u µ the 4-velocity of the fluid. We split energy density and pressure into a non-relativistic and a relativistic part as ρ = ρ m + ρ r and p = p m + p r , resp. The subscript m refers to non-relativistic matter, while the subscript r refers to radiation. The conservation of energy-momentum (2.6) leads to the continuity equatioṅ with w i = p i /ρ i the equation of state. Therefore, the energy density of matter and radiation evolve with the scale factor a of the physical metric as in standard cosmology. To be explicit, the continuity equation is solved by (2.20) Here we defined the energy density of the dynamical dark energy induced by the bimetric potential as ρ de m 2 g = β 0 + 3β 1 y + 3β 2 y 2 + β 3 y 3 (2.21) with the time-dependent function y. The Friedmann equation of f µν is sourced by potential energy only with energy density Interpreting the effect of spatial curvature as an energy source, we can define the energy density Then, the Friedmann eq. (2.20) of g µν can be written as It remains to determine the time-evolution of ρ de , which in general cannot be written as a polynomial in 1/a. Subtracting the the Friedmann equations in (2.20) yields This equation determines the time evolution of y in terms of the time evolution of the matter energy densities. It is a quartic polynomial in y, hence there are up to four real-valued solutions with different time evolutions of y. Following [60,67], we can distinguish three classes of solutions by inspecting the early-time behavior, for which the energy densities classically diverge: • Finite branch. The first possible early-time behavior is that y approaches zero so as to compensate the diverging value of the energy density in the linear piece of the polynomial (2.25). Hence, during early times the scale factor ratio can be approximated as y ∼ m 2 g β 1 /(α 2 ρ). Existence of the finite branch requires β 1 > 0. The scale factor ratio monotonically increases in time. In the asymptotic future, the energy densities vanish and y approaches a constant value c, which corresponds to the lowest-lying, strictly positive root of eq. (2.10) [105]. For positive matter and radiation energy densities, this branch of solution satisfies the cosmological stability bound, a generalization of the Higuchi bound to FLRW spacetime [79,126].
• Infinite branch. The alternative is that y is large as the energy densities are large. Then, at early times, the evolution can be approximated as y 3 ∼ −ρ/(m 2 g β 3 ). This reveals that the infinite branch exists only if β 3 < 0. The scale factor ratio monotonically decreases in time. In the asymptotic future, when the energy densities vanish, y approaches the highest lying, strictly positive root of eq. (2.10) [105]. This branch, however, inevitably violates the cosmological stability bound and is therefore not a physical solution [79,126].
Due to the high degree of the polynomial (2.25), there are further exotic branches of solutions, which however do not give rise to a viable expansion history [79]. Therefore, we pick out the finite branch solution in section 3. The scale factor ratio then satisfies 0 < y < c [79,105].
Let us pause for a moment and summarise. The bimetric potential contributes a dynamical dark energy to the Friedmann equation in the form of ρ de , which depends on the scale factor ratio y. Since y approaches a constant value c in the asymptotic future, ρ de approximates the effect of a cosmological constant at late times. We have identified two possible time evolutions for y yielding a consistent expansion history. Out of these, only the finite branch solution satisfies the cosmological stability bound. Therefore, y decreases when looking back in time, implying that also ρ de decreases. In other words, the (self-)interactions of the massive and massless spin-2 field are dynamically increasing as time proceeds. Hence, the (self-)interaction energy is necessarily phantom or negative [69,79].

Physical parametrisation
The solutions presented thus far are parametrised in terms of the Planck mass ratio α and the interaction parameters β n . Vacuum solutions are labeled by c that are the roots of eq. (2.10). The bimetric action and hence the equations of motion are invariant under a rescaling of these parameters implying that one of them is redundant.
In [105] a parametrisation that is invariant under this rescaling is proposed and worked out. It relies on the following physical parameters: the mass m FP of the spin-2 field, its coupling to matter α, and the effective cosmological constant Λ. The hence baptised physical parametrisation allows to unambiguously combine various theoretical and observational bounds on the parameter space of bimetric theory.
Let us give more detail on the physical parametrisation and on the already inferred theoretical consistency bounds. Since there are up to four real-valued roots c of eq. (2.10), the relation between the physical and interaction parameters is not unique. However, a vacuum solution and the cosmic expansion history are well-defined only under the following conditions: • Consistent de Sitter vacuum. Without loss of generality we restrict ourselves to roots with c > 0 9 . We impose that the mass of the spin-2 field is positive, m FP > 0, and we restrict ourselves to de Sitter vacua, Λ > 0. For a massive spin-2 field propagating in de Sitter space, unitarity forbids the mass to be arbitrarily small. Instead, the Higuchi bound [127] 3m 2 FP > 2Λ (2.26) has to be satisfied. If the bound is violated, the helicity-0 mode of the massive spin-2 field has a negative norm (Higuchi ghost).
• Consistent asymptotic future. The cosmic expansion has to approach a consistent de Sitter vacuum in the asymptotic future. Since only the finite branch solution to the Friedmann equations is physical, the scale factor ratio satisfies 0 < y < c. Hence, the lowest-lying strictly positive root of eq. (2.10) as the unique physical vacuum solution. Identifying the asymptotic future of the cosmic expansion history with the physical vacuum solution implies further, model-dependent conditions on the bimetric parameters.
• Consistent cosmic expansion history. The consistency of the early universe for a cosmic expansion history on the finite branch requires [67] β 1 > 0 , (2.27) as can be seen from eq. (2.20). Otherwise H 2 ≤ 0 during early times.
These restrictions imply a unique relation between the physical and interaction parameters. The details for all bimetric models with up to three non-vanishing interaction parameters are worked out in Ref. [105]. For completeness, we summarise the explicit expressions in table 4. Using the physical parametrisation, in section 3 we compute bounds from cosmological observations.

Bounds from cosmological observations
We now compare bimetric theory to cosmological observables. We restrict ourselves to background observables that constrain the Hubble rate at different redshifts. In particular, we use Supernovae type 1a, Baryon Acoustic Oscillations and the first peak of the oscillations in the Cosmic Microwave Background, which we discuss in greater detail in section 3.2.

Parametrisation used for model fitting
We have already explained how to get the time-evolution of the scale factor ratio y, the time evolution of ρ de and hence the time evolution of H. To ease these computations, it is convenient to introduce the energy density parameters The Friedmann equation (2.24) can then be written in the compact form Evaluating this equation at present times yields the following relation between the parameters as We choose to express the parameter Ω m0 in terms of the other parameters in our statistical analysis. In order to compute Ω de , we introduce the following constant parameters: The relations between B n and the physical parametersᾱ, Ω FP , and Ω Λ are presented in appendix A as derived in [105]. Note that B n is rescaling-invariant and hence observable. Further, in analogy to the definition ofᾱ, we define the rescaled scale factor ratio asȳ = αy. To be explicit, the energy density parameter of the induced dynamical dark energy is given as where E = H/H 0 is the normalised Hubble rate, as follows from eq. (2.21). The scale factor ratioȳ is determined by the quartic polynomial (2.25), which in terms of the energy density parameters reads As last ingredient, we need to determine Ω de0 as initial condition. In order to do so, we determine the value of the scale factor ratio at present timeȳ 0 using the f µν -Friedmann equation (2.20) evaluated today: This polynomial has up to three real-valued roots, out of which we pick the one that satisfies 0 < y 0 <ᾱ in order to ensure that we are on the finite branch. Plugging the result into eq. (3.5) evaluated today, determines Ω de0 , which in turn determines Ω m0 via eq. (3.3).

Cosmological data
We now present the cosmological observables that are used in our statistical analysis. We use three different types of observations: Supernovae type 1a (SN1a), Baryon Acoustic Oscillations (BAOs) and Cosmic Microwave Background (CMB). Before discussing them in detail, let us first introduce some useful quantities to enlighten the notation. We start with the definition of the following integral: where E(z) = H(z)/H 0 . We use this quantity to write the various distance indicators, such as the luminosity distance d L and the angular diameter distance d A : Another important quantity is the comoving sound horizon, defined as . (3.10) Here Ω b0 and Ω γ0 are the baryons and photons density parameters. Ω γ0 is fixed by the temperature T CMB of the CMB [129]: Here, h = H 0 /100 km s −1 Mpc −1 is the normalised Hubble rate today. The photons density parameter is related to the total radiation density parameter via N eff is the effective number of neutrino species, that we fix to the standard value N eff = 3.046 for our analysis. We now proceed introducing separately the different observables.

Supernovae type 1a
Supernovae of type 1a serve as standard candles whose luminosities can be inferred independently of their redshift. This allows to build up the local distance ladder. We use the catalogue of the Joint Light-curve Analysis (JLA) that contains 740 SN1a events as reported in Ref. [130]. In order to compare the luminosity distance to the observed apparent magnitude m obs , we have to callibrate the supernovae as [130] where m is related to the luminosity distance as in Eq. (3.13). Here, X 1 and C are the light-curve parameters of the supernova. For the nuisance parameters α, β and ∆ M we use the best-fit values α = 0.140 ± 0.006 and β = 3.139 ± 0.072 as obtained in Ref. [130]. The parameter ∆ M is a correction to the absolute magnitude M of the supernova that depends on the stellar mass M * of the host galaxy of the supernova as Marginalising over M, the log-likelihood for the SN-data is given by where we have defined in terms of y = m − 5 log 10 D L . The covariance matrix C SN is given in Ref. [130]. The errors on α, β and ∆ M are added in quadrature.

CMB
CMB observations constrain background cosmology through the measurement of two different distance ratios [131]. The first is the ratio between the angular diameter distance d A (z * ) and the sound horizon r s (z * ) at decoupling epoch z * . This ratio is measured through the location of the first peak of the CMB power spectrum: The second is the angular distance divided by the Hubble horizon at the decoupling epoch, d A (z * )H(z * )/c. This last information is usually implemented in the shift parameter, defined as: which differs from d A (z * )H(z * )/c by a factor of √ 1 + z * . This parameter is obtained assuming which correspond to neglect all contributions (other than Ω m0 ) to the evolution history at the time of decoupling H(z * ). We use that the early-universe is not altered in bimetric theory compared to general relativity as our working assumption. This is justified because on the finite branch the energy density arising from the bimetric potential does not contribute to the Hubble rate at large redshifts. For our statistical analysis we follow the standard procedure [132], implementing CMB data in the three distance priors (l A , R, Ω b0 h 2 ). We use the latest Planck values [129] In our analysis R(z * ) and l A (z * ) are computed from eq. (3.19) and eq. (3.18), where the decoupling epoch z * is obtained via the fitting function [133]: Defining the CMB data vector as the likelihood for the CMB is written as: where the covariance matrix C CMB can be found in [129].

BAO
Barionic Acoustic Oscillations (BAOs) are another important tool to constrain background cosmology. We follow a procedure similar to [107], using measurements at 10 different redshifts z = 0. 106 The other relevant quantities in our study are the effective distance measure d V and the redshiftweighted comoving distance d M : The likelihood for the statistical analysis is derived using the covariance matrices taken from the original references [134][135][136][137][138].

Statistical analysis and results
Using the aforementioned data sets, we perform a statistical analysis via MCMC sampling for all models with up to three non-vanishing interaction and vacuum energy parameters β n . The ΛCDMmodel corresponds to the β 0 -model in this sense. We classify models according to their number of free parameters β n : models with one, two or three free parameters β n are referred to as one-, twoor three-parameter models, respectively. In terms of the physical parameters, that means that one, two or all three of the parametersᾱ, Ω FP and Ω Λ are independent. The best-fit values and the 1-dimensional (1σ) error intervals for all free and derived parameters are reported in tables 1 and 2.
In order to estimate the goodness-of-fit for each model, we compute the Bayesian Information Criterion (BIC) [139], where χ 2 = −2 ln L is evaluated at maximum likelihood, k is the number of model parameters and N is the number of data points. For our combined data set SN+CMB+BAO we have N = 740+3+10 = 753. The number of model parameters is given by k = 3+k withk the number of independent physical parameters as will become clear from the next paragraph. We will use the BIC of the ΛCDM-model as reference. The difference ∆BIC allows to estimate the preference of one model over another as [107,140]: strong support (∆BIC < −12), favorable (∆BIC < −6), inconclusive (∆BIC < 6), disfavored (∆BIC < 12), strongly disfavored (∆BIC ≥ 12). The corresponding ∆BIC are reported in tables 1 and 2 as well. Independent parameters and priors. We use the theoretical bounds of [105] that ensure a consistent cosmology, which we summarise in table 4 Reference model and consistency of data sets. Let us start with a discussion of the results for the ΛCDM-model as reference. The 1-and 2-dimensional posterior distributions for the parameters H 0 , Ω m0 and Ω Λ = Ω de0 are presented in the left panel of fig. 1. While supernovae do not constrain H 0 due to the degeneracy with the absolute magnitude, also CMB and BAOs do not yield robust constraints individually. However, this geometric degeneracy is broken when combining all three sets of data. We then obtain the value H 0 = (68.8 ± 1.2)km/s/Mpc, which is in good agreement with current constraints [5,130,141]. Our value is slightly but not significantly larger, which can be understood as an artifact of using the distance priors [129]. Combining the data sets stabilises spatial curvature to the value Ω k0 = 0.002 ± 0.003 such that the universe is spatially flat within 68% c.l. Also the other parameters as reported in table 1 are in good agreement with current constraints. We conclude that our analysis is robust. At the best-fit point we have χ 2 = 694.7 leading to BIC = 721.2.

β 1 -model
As discussed before, the β 1 -model is the only bimetric one-parameter model that can possibly give rise to a viable expansion history. The posterior distributions are presented in the right panel of fig. 1. Also this model is subject to the geometric degeneracy such that both H 0 and Ω k0 are unconstrained by CMB and BAO individually. Combining all three data sets stabilises these parameters. Data then favors a slightly spatially curved universe with Ω k0 = −0.007 +0.004 −0.003 . The Hubble rate today is constrained to the remarkable high value of H 0 = 73.1 +1.6 −1.4 km/s/Mpc, which is in perfect agreement with local determinations of H 0 [142]. We emphasise that we obtained this result without local prior  on H 0 . The value of the spin-2 mass is constrained to be small with m FP = (1.82 ± 0.02) × 10 −32 eV, which is close to the Higuchi bound. However, the cosmological early-and late-time measurements are incompatible within this model, as can be seen already from the right panel of fig. 1. Quantitatively, ∆BIC = +31 indicates that this model is indeed strongly disfavored by cosmological background data. Therefore, this model does not solve the H 0 -tension, despite its large favored value. We conclude that the bimetric β 1 -model is statistically ruled out. Nonetheless, we summarise the obtained best-fit values for the other parameters in table 1.

Two-parameter models
As already discussed, consistency of the early universe on the finite branch requires β 1 > 0. This means that there are four possibilities for the two-parameter models: the β 0 β 1 model and the β 1 β n models with n = 2, 3, 4. From our analysis it emerges that the results for all the β 1 β n models are very similar and statistically compatible. For this reason, we explicitly show only the results for the β 1 β 3 model and we refer to it as β 1 β n for the rest of this section. In fig. 2 we plot the 1-and 2dimensional marginalised posterior distribution for these models, where we can see that the data sets are compatible. This is true also for SN data, that are not shown explicitly in this figure. We conclude that the data sets can be combined in the statistical analysis.
In table 1 we summarise the results of the combined analysis. The main difference between the β 0 β 1 model and the β 0 β n models are the values ofᾱ and m FP . In both cases,ᾱ is constrained only from above with a maximum value of 0.2 for the β 0 β 1 model and 0.016 for the β 1 β n models. In the β 0 β 1 model the value of m FP is constrained to be (1.33 ± 0.02) × 10 −32 eV, which lies close to the Higuchi bound. In the β 1 β n models, however, the value of the Fierz-Pauli mass is constrained only weakly from below with m FP > 4.24 × 10 −31 eV. We will discuss the constrains onᾱ and m FP in more detail in section 4.3, where we will compare the results of the statistical analysis with local tests of gravity. Despite these differences, the other parameters for the β 0 β 1 and β 1 β n models are very similar with each other and compatible with ΛCDM. The only exception is spatial curvature, that seems to be slightly positive for the β 0 β 1 model. However, this model still prefers an almost flat universe with Ω k0 = 0.002 +0.005 −0.001 , which remains very close to 0 at 1σ and compatible with the ΛCDM value. Moreover, the value of χ 2 in both models coincides with ΛCDM. This makes the β 0 β 1 and β 1 β n models slightly disfavored with a ∆BIC = +6.6, since these models have more free parameters with respect to ΛCDM.

Three-parameter models
We finish by discussing models with three free parameters. We will see that the parameter constraints are similar across these models. Data sets are compatible for all the three parameter models, as we can see in figs. 3 to 5. Therefore, combining the data sets for these models is trustworthy 11 .
• Models with β 4 = 0. We first discuss the three-parameter models without vacuum energy in the f -sector. We present the 1-and 2-dimensional marginalised posterior distribution in the left and right panel of fig. 3 and the left panel of fig. 4 for the β 0 β 1 β 2 -, β 0 β 1 β 3 -and β 0 β 1 β 4 -model, respectively. The parameter constraints are reported in table 2 and agree with the ΛCDM values within 1σ and are all quite similar across these models. The only exception is that the β 0 β 1 β 2 -model seems to prefer a slightly curved universe at 1σ. However, the value of the curvature is very small and compatible with ΛCDM at 1σ and therefore we do not consider it as a statistically relevant difference. The bimetric parametersᾱ and m FP are constrained only weakly. We will discuss more about the constraints on these parameters in theᾱ − m FP plane when comparing with local tests in section 4.3. We find that these models are consistent with the current data because χ 2 = 694.7 as in the ΛCDM-model. However, the BIC is higher due to the two additional free parameters indicating that these models are statistically disfavored.
• Models with β 0 = 0. The result for this models is very similar to models with β 4 = 0. The only exception is thatᾱ is more constrained of about one order of magnitude, as it can be seen in table 2.

Bounds from local tests of gravity
Having identified those regions of the parameter space that are in agreement with current observations in background cosmology, we now turn to local tests of gravity. In section 4.1 we give a brief review on local tests of gravity and present the existing constraints on the parametersᾱ and m FP . However, the existing bounds were derived without taking into account the Vainshtein screening mechanism. In regions where Vainshtein screening is active, the bounds onᾱ and m FP are not trustworthy. In section 4.2 we identify the regions of the parameter space that give rise to Vainshtein screening in spherically symmetric systems. We finally confront cosmological with local bounds in section 4.3 in parameter regions without Vainshtein screening.

Local tests of gravity
In this section we briefly summarise local tests of gravity. For details and reviews we refer to Ref. [108][109][110][111]113]. Extra degrees of freedom (such as the massive spin-2 field of bimetric theory) modify the Newtonian gravitational potential as compared to GR. This is usually attributed to a fifth force. If the force mediator is massive, it contributes a Yukawa potential to the gravitational potential felt by a massive test body. Following standard notation in the literature, the gravitational potential for a static and spherically-symmetric configuration around a massive source can be parametrised as with m P the Planck mass. This is referred to as Yukawa parametrisation. The first term corresponds to the usual potential arising from massless gravitons. The second term arises from a massive force mediator with Compton wavelength λ and coupling to matter ξ. While in the limit ξ 1 the fifth force is suppressed, it is most prominent on length scales r ∼ λ. As we will see in section 4.2, these parameters are related to the bimetric parameters as In fig. 6 we present the constraints from local tests of gravity on the bimetric parametersᾱ and m FP as indicated by the black lines. The gray-shaded region is excluded by 95% c.l. To compute the constraints, we follow the procedure as described in [108][109][110]. The laboratory constraints correspond to the Irvine [143] 12 , Eöt-Wash [145] and Stanford [146] experiments. The geophysical constraints correspond to the lake [147,148] and tower [149][150][151] experiments. The constraints from measurements ��� �� (α) Figure 6. Bounds on the bimetric parametersᾱ and mFP from local tests of gravity. The black lines result from tests of the Newtonian gravitational potential as reported in [108-110, 145, 146, 155]. Tests of the scalar curvature yield the red line [91,112]. The gray-shaded region is excluded by 95% c.l. Note that the Vainshtein mechanism is not taken into account in the derivation of these bounds.
with the LAGEOS satellite are reported in [152][153][154]. The constraint from Lunar-Laser-Ranging (LLR) results from lunar precession [109]. The planetary constraints are reported in [155]. The aforementioned tests of the Newtonian potential provide only weak constraints on modifications if the mass of the fifth force mediator is small. In this region of the parameter space, tests of the scalar curvature, i.e. the deflection and delay of light induced by matter, provide powerful constraints. We use the gravitational slip parameter γ, which is unity in GR, γ = 1. Within the context of bimetric theory, the parameter has been calculated to be [91] γ = 3 + 2ᾱ 2 e −mFPr 3 + 4ᾱ 2 e −mFPr . The most precise value has been obtained within the solar system. The measurement of the time delay of radar signals sent between the Earth and the Cassini spacecraft on its way to Saturn yields γ − 1 = (2.1 ± 2.3) × 10 −5 [112]. The radio signals were passing by the sun at a distance of 1.6 solar radii, i.e. r 7.44 × 10 −3 AU 1.31 × 10 −15 eV −1 . Using these values defines the red line in fig. 6. The region above that line is excluded by 95% c.l. For small spin-2 masses, the coupling to matter is constraint to beᾱ 1.7 × 10 −3 . A collection of bounds on γ from time delay and deflection of light form galaxy halo to solar system scales can be found in [113], which however are less constraining than the bound used here. Summarising, the most stringent constraint comes from LLR withᾱ < 8.9 × 10 −6 at 95% c.l. for m FP 6.5 × 10 −15 eV. On the other hand, for m FP 10 −2 eV the bounds onᾱ are weak. This, however, is not the end of the story in the framework of bimetric theory due to the Vainshtein screening mechanism as we discuss in the next section.

Vainshtein screening
Local tests of gravity put tight constraints on modifications of the gravitational potential. Historically, the linear theory of Fierz and Pauli [39] was discarded because it did not pass basic solar system tests. This is due to the so-called vDVZ discontinuity [40,41]: the helicity-0 mode of the massive graviton couples to the trace of the energy-momentum tensor even in the limit m FP → 0. Vainshtein argued that including nonlinear interaction terms for the massive spin-2 field restores GR in that limit thus curing the problem [42]. This is indeed the case as explicitly demonstrated in [156][157][158]. For a review on the Vainshtein mechanism we refer to [159].
Bimetric theory incorporates the Vainshtein mechanism as explicitly demonstrated for static and spherically symmetric systems [89,90,93,160]. We leave the technical details to appendix B, where we repeat and extend the analysis by relaxing the assumption of asymptotic flatness to allow for a non-vanishing cosmological constant. Here we limit ourselves to discussing the general features. The radius below which nonlinearities have to be taken into account is referred to as Vainshtein radius r V = (r S /m 2 FP ) 1/3 with r S the Schwarzschild radius of the matter source. For larger radii, r r V , the linear approximation is valid. The resulting gravitational potential for small and large radii can be written as [89,93] where m 2 P = (1 +ᾱ 2 )m 2 g is the unscreened Planck mass. This justifies our identification (4.2). On the other hand, since r V depends on the spin-2 mass and the mass of the central source, the constraints onᾱ and m FP summarised in fig. 6 are not trustworthy.
However, Vainshtein screening exists only in a subregion of the parameter space [90,93], as we review in appendix B. In this section, we will identify this subregion in terms ofᾱ and m FP . The Vainshtein mechanism relies on the following three necessary conditions: • Consistent screening regime. The existence of a nonlinear solution that restores GR in spacetime regions close to the matter source requires Since the denominator is positive if m 2 FP > 0, this bound implies that β 3 must be strictly positive 13 , β 3 > 0. That means, any bimetric model with β 3 = 0 does not have a working Vainshtein mechanism (for some related caveats, see [93]).
• Consistent asymptotics. The nonlinear equations must give rise to a solution that matches the linearised solution for r r V . This requires to be satisfied. Combined with eq. (4.5) it follows that β 2 must be strictly negative, β 2 < 0. Bimetric models with β 2 = 0 do not give rise to Vainshtein screening.
• Consistent Vainshtein-Yukawa solution. The nonlinear solution realising Vainshtein screening has to be smoothly connected to the linear solution realising the Yukawa-type fifth force without branch cuts. The existence of such a well-defined Vainshtein-Yukawa solution imposes another restriction on the parameters. The explicit expression is lengthy so that we shift its presentation to eq. (B.28).
Following [105], we express these conditions in terms of the physical parametersᾱ, m FP , and Λ in this section. The relations between the interaction and physical parameters are collected in appendix A. It suffices to study models with all interaction parameters β 1 , β 2 and β 3 being free as only those models can possibly give rise to a viable cosmic expansion history while incorporating the local Vainshtein mechanism. Cosmologically viable one-and two-parameter models are not able to incorporate Vainshtein screening. The β 1 β 2 β 3 -model is the only three-parameter model with a possibly viable background cosmology and a working Vainshtein mechanism. In fig. 7 we collect all the aforementioned theoretical bounds on the physical parameter space of that model. The expressions are too lengthy to display them here explicitly. The blue-and red-shaded regions are excluded as these do not give rise to a viable expansion history. The purple-shaded region violates the bound in eq. (4.5). In the greenshaded region the bound in eq. (4.6) is violated. The Vainshtein-Yukawa solution does not exist in the yellow-shaded region, which represents the most stringent bound. To summarise, we expand the most-stringent bound on the physical parameters for m 2 FP Λ to find If this approximate bound is satisfied, the β 1 β 2 β 3 -model incorporates the local Vainshtein mechanism. As can be seen from fig. 7, the parameter region giving rise to a viable cosmic expansion history is only slightly larger. Before moving on, let us point out the following caveats. The equivalence principle is violated in two-body systems due to nonlinearities [161]. The assumption of spherical symmetry entering the derivation of the gravitational potential might not be appropriate to describe, e.g. the earth-moon system. Further, Vainshtein screening is weaker or even completely absent in systems with cylindrical or planar symmetry, respectively, as demonstrated in the context of galileons [162]. This has to be taken into account for gravity tests without spherical symmetry.

Local vs. cosmological bounds
We are finally in the position to compare the constraints from local tests of gravity, shown in fig. 6, to the parameter regions that give rise to Vainshtein screening. Further, we compare the local constraints  to the constraints that we inferred from background cosmology in section 3.3. Since only the β 1 β 2 β 3model incorporates Vainshtein screening, the local bounds are directly applicable to all the other models considered in this paper. The discussed observational bounds in this section always refer to 95% confidence level. Let us start discussing the β 1 -model. As mentioned before, the cosmological data sets are not compatible, strongly disfavoring this model. If we nonetheless combine the data sets, the spin-2 mass is constrained to m FP = (1.82 ± 0.02) × 10 −32 eV. The coupling to matter is fixed to the value  Figure 11. The combined bounds from local and cosmological tests for the β1β2β4-and β1β3β4-model. The region above the gray line is excluded by local tests of gravity at 95% c.l. The red-shaded region is excluded by theoretical constraints. The region above the blue line is excluded by our bounds from cosmology at 95% c.l. ᾱ = 1/ √ 3 0.58, which is incompatible with the observational bound from Cassiniᾱ 1.7 × 10 −3 in that mass regime. Since this model does not give rise to Vainshtein screening, we conclude that the β 1 -model is completely ruled out by the combined analysis of cosmological and local data.
Next let us discuss the two-parameter models, which do not give rise to Vainshtein screening as well. This implies that the bounds from local tests of gravity are indeed trustworthy. In fig. 8 we compare the bounds from local and cosmological tests for the two parameter models β 0 β 1 and β 1 β 3 . The blue line defines the boundary of the region that is consistent with cosmological data. The gray-shaded region is excluded by local tests of gravity. For the β 0 β 1 -model, cosmological data constraints the spin-2 mass to m FP = (1.33 ± 0.02) × 10 −32 eV while Cassini constrains the coupling to matter toᾱ 1.7 × 10 −3 , which is slightly more stringent than the cosmological bound. For the β 1 β n -model, the physical parameters are correlated as can be seen in the right panel of fig. 8. Indeed, they are approximately related asᾱ 2 m 2 FP 12 n Λ for n = {2, 3, 4} (see appendix A.1) [105]. As before, Cassini provides the most stringent constraint onᾱ. Due to the parameter correlation, this implies a more stringent constraint on the spin-2 mass than inferred from cosmology of m FP 1.85 × 10 −30 eV at 95% c.l. Therefore, all these two-parameter models are in agreement with observational constraints from cosmology and local tests of gravity, but driven to their GR-limits.
In figs. 9 to 11 we combine theoretical, cosmological and local bounds in single plots for all the three parameter models. In all these plots, the regions above the blue lines are excluded at 95% c.l. by our cosmological bounds inferred in section 3.2. The red-shaded regions in the plots are excluded by theoretical constraints that ensure a consistent cosmic expansion history, which we used as hard priors in the cosmological analysis. The region excluded by local bounds at 95% c.l. (if the Vainshtein mechanism is not taken into account) is shown as a gray-shaded region. For large spin-2 masses, the bounds from cosmology are the most stringent ones. For smaller spin-2 masses, the most stringent bound is provided by tests of the scalar curvature with the Cassini spacecraft. Summarising, without taking Vainshtein screening into account, all three-parameter models are consistent with both local tests of gravity and background cosmology. However, since the coupling of the massive spin-2 field to matter is forced to be small, i.e.ᾱ 1.7 × 10 −3 , all these models are essentially forced into their GR-limits. Although local tests of gravity allow for large couplings to matter if the spin-2 mass is large, m FP 10 −2 eV, cosmological data excludes this parameter region. Therefore, deviations from GR are substantially suppressed for these models on all scales.
Let us discuss further about the cosmological bounds onᾱ and m FP . From inspecting figs. 9 to 11 we notice that the 95% contour line falls into two regimes. The first is the small spin-2 mass regime. Here, below some critical value that depends on the model, the observational constraints becomes independent of m FP . This is particularly evident for the models with β 0 = 0. In this regime, the coupling to matterᾱ is bounded from above as reported in tables 1 and 2. The second regime is above the critical value for m FP , where the contour line can be approximated as: We give the explicit value of the coefficients c 1 and c 2 in table 3, as the result of a polynomial fit. This behavior can be seen in figs. 9 to 11, that however show only a portion of theᾱ − m FP plane scanned by the cosmological analysis 14 . Indeed, the approximation (4.8) is valid up to the Planck scale, i.e. beyond the zoomed-in regions in the plots 15 . For the β 1 β 2 β 3 -model, we need to take into account the Vainshtein mechanism when including local tests of gravity. The parameter region, which supports Vainshtein screening, is depicted as a hatched area in the right panel of fig. 10. As we have seen before, the region incorporating Vainshtein screening is almost exactly complementary to the region that is excluded by our theoretical bounds. In particular, the entire parameter region that is consistent with cosmological data gives rise to Vainshtein screening. From the right panel of fig. 10 we see that there is a small overlap between the region excluded by local tests and the region where the Vainshtein mechanism is active. In this interesting region, which is characterized by a potentially large value ofᾱ, local tests are not Model c 1 c 2 β 0 β 1 β 4 0.97 ± 0.01 23.6 ± 0.4 β 0 β 1 β 3 1.117 ± 0.006 32.9 ± 0.2 β 0 β 1 β 2 0.964 ± 0.006 24.6 ± 0.2 β 1 β 2 β 3 1.047 ± 0.007 35.5 ± 0.2 β 1 β 2 β 4 1.021 ± 0.006 32.7 ± 0.2 β 1 β 3 β 4 0.959 ± 0.008 32.7 ± 0.3 Table 3. The value of numerical coefficients for the fitting formula (4.8), obtained with a polynomial fit. We also show 1σ errors in the determination of the coefficients.
trustworthy. Unfortunately this small region seems to be excluded by cosmological data at 95% c.l. in our statistical analysis. However, this is probably caused by prior domination. Indeed, this is a very small region of the full parameter space explored by the statistical analysis and it is close to the hard prior given by the theoretical constraints. To further investigate this issue, one should perform another statistical analysis focusing on this small region of parameter space.
For the models with β 0 = 0, the 95% contour lies very close to the theoretical priors except for a few wiggles. We interpret the wiggles as an artifact of the finite size of the Markov chains and conclude that the used data sets are not constraining these parameters. Instead, the statistically favored regions almost exactly coincide with the parameter regions allowed by our theoretical priors. To stabilise the value ofᾱ and m FP , we would need to include other data sets, which however is beyond the scope of the present analysis.
On the other hand, models with β 4 = 0 are subject to rather weak theoretical priors. Nonetheless, the 95% lines show the same tendency as for the other models. Hence, the data sets constrain these parameters substantially beyond the prior knowledge, but is not able to stabilise their values. To understand why the upper-right region in theᾱ−m FP -plane is excluded also for these models, we note the following. Expanding the expression for β 0 in the limit m 2 FP Λ, we find that β 0 − 12 nᾱ 2 m 2 FP for the β 0 β 1 β n -model with n = {2, 3, 4}. This allows to conclude that data disfavors large negative values of β 0 . This means that, in order to be cosmologically viable, ρ de is allowed to change its sign only at sufficiently large redshifts.

Conclusions and outlook
Massive spin-2 fields might play a crucial role in our universe. Their (self-)interaction energy can drive the observed cosmic accelerated expansion at late times. Furthermore, a massive spin-2 field serves as ideal dark matter candidate as it interacts with ordinary matter only gravitationally [84][85][86]. Recently, another intriguing aspect has been explored within a completely different context. As argued in [163], massive spin-2 states could be essential in solving the black hole information paradoxon.
In this paper, we have further studied the observational viability of bimetric theory. We confronted all models with up to three free interaction parameters β n with measurements of SN1a, BAOs and CMB. We find that all two and three parameter models are in perfect agreement with these data sets. Compared to the ΛCDM model, bimetric models are statistically disfavored due to the larger number of free parameters. However, the self-accelerating solutions with β 0 = 0 appear more appealing from a theoretical perspective as these models do not suffer the cosmological constant problems [12,13].
The best-fit values for the cosmological parameters do not significantly differ from the best-fit values of the ΛCDM model. We have also studied spatial curvature within bimetric theory, and our results show that a spatially flat universe is always preferred.
The only exception from the previous summary is the β 1 -model, which represents the simplest bimetric model that gives rise to self-accelerating solutions on the finite branch. Our statistical analysis shows that the aforementioned data sets are inconsistent leading to a substantially higher χ 2 at the best fit point. We therefore conclude that the β 1 -model is statistically disfavored. This conclusion is further strengthened when simultaneously confronting this model with local tests of gravity. Cosmological data selects a spin-2 mass of m FP = (1.82 ± 0.02) × 10 −32 eV, for which local tests of gravity demand a coupling to matter ofᾱ 1.7 × 10 −3 at 95% c.l. Within the β 1 -model, this parameter is fixed to beᾱ = 3 −1/2 . Therefore, this model is in severe conflict with observations.
We also analysed the other bimetric models with up to three free parameters, for which we compared the bounds from background cosmology to the ones obtained from local tests of gravity. All models are consistent with both types of observations even if the Vainshtein mechanism is not active. However, the coupling of the massive spin-2 field to ordinary matter is forced to be small, α 10 −3 (see section 4.1 and tables 1 and 2 for the precise values). For large spin-2 masses, the cosmological observations force the coupling to be even smaller. This suppresses the deviations from GR on all scales for these models.
For the β 1 β 2 β 3 -model, the Vainshtein mechanism removes the constraint onᾱ and m FP implied by local tests for sufficiently small values of the spin-2 mass. In this region, whereᾱ can be large, the theory is not in its GR limit. As we mentioned in section 4.3, the fact that this small region is excluded by our cosmological analysis might be a numerical effect, and it would be interesting to further investigate this issue in another analysis. Further, more general bimetric models, i.e. those with four or all five parameters β n being free, might not be forced into their GR-limits. In any case, even if the restricted models studied in this paper turn out to be forced into their GR-limits, their theoretical motivation would persist: self-accelerating models with β 0 = 0 have an effective cosmological constant that is technical natural in the sense of t'Hooft [164].
The data sets used in this paper are only weekly constraining the physical parameter space ofᾱ and m FP . To stabilise their values and effectively constrain bimetric theories, a possibility would be taking into account further measurements. Cosmological perturbations around the FLRW background solution are expected to have constraining power as these probe different redshifts and scales. However, in bimetric theory linear scalar perturbations are plagued by gradient instabilities at early times [65,70,74,[77][78][79]. The instability sets in when the Hubble rate is of the order of the spin-2 mass, H m FP [69]. Nonlinear terms become as important as the linear term rendering perturbation theory invalid. Indeed, taking nonlinearities into account might stabilise the solution in analogy to the Vainshtein mechanism [69,81,82,128]. Structure formation and nonlinear cosmological perturbations are therefore main open issues within bimetric theory, which should be addressed in the future.

A Details on the physical parametrisation
In order to provide an overview over the relation between the physical and interaction parameters, we summarise the results of Ref. [105] in this appendix. We immediately work in terms of the cosmological parameters that we use for parameter fitting.
In the β 1 -model, there is only one free physical parameter, in terms of which the interaction parameter is given as while the mixing angle is fixed to beᾱ = 1/ √ 3.

A.1 Two-parameter models
Next, we turn to the two-parameter models with β 1 = 0. We choose to express the interaction parameters on terms of Ω FP and Ω Λ , while the mixing angleᾱ is a derived parameter. However, one can equivalently eliminate one of the physical parameters byᾱ. For the β 0 β 1 -model, the relation between the parameters is in terms of which the mixing angle is determined to beᾱ = Ω FP /Ω Λ − 1.

A.2 Three-parameter models
We continue by stating the parameter relations for the three-parameter models with β 1 = 0. For the β 0 β 1 β 2 -model, the relations read (A.10)

Model
Viable cosmology Vainshtein screening Table 4. This table summarises the theoretical constraints on all bimetric models with up to three nonvanishing interaction parameters. The cosmology constraints were derived in Ref. [105] and follow from requiring a consistent cosmic expansion history. We only state the most stringent constraints other than the Higuchi bound, 3ΩFP > 2ΩΛ. Only the β1β2β3-model has a working Vainshtein mechanism. We report the approximation of the most stringent bound.
In the β 0 β 1 β 3 -model, the physical and interaction parameters are related as Next, the β 0 β 1 β 4 -model gives rise the the following relations, Turning to models with β 0 = 0, the parameter relations for the β 1 β 2 β 3 -model read Next, the β 1 β 2 β 4 -model has the following parameter relations, (A.14) -26 -Finally, in the β 1 β 3 β 4 -model the interaction and physical parameters are related as

B Vainshtein screening in static, spherically symmetric systems
In this appendix, we will solve the bimetric equations of motion for the static, spherically symmetric, and bidiagonal ansatz. This allows to study the Vainshtein mechanism in these systems, which restores General Relativity on small length scales. However, a solution that incorporates the Vainshtein mechanism does not exist for every set of bimetric parameters. Instead, demanding existence of such a solution implies conditions on the bimetric parameters. The situation was studied in detail in Refs. [90,93] however with flat asymptotics. We aim at generalizing the analysis allowing for a non-zero cosmological constant [166] and compute the conditions for the Vainshtein mechanism to work in this more general setup. As we will see, our results reduce to the ones of Ref. [90,93] in the limit of vanishing cosmological constant. Assuming staticity and spherical symmetry, the metrics g µν and f µν can be written in the bidiagonal case as [93] ds 2 g = −e −νg dt 2 + e λg dr 2 + r 2 dΩ 2 , (B.1) The metric functions ν g,f , λ g,f , and µ depend on radius r only. The field µ can be thought of as a Stückelberg field that restores diffeomorphisms in the r-dimension.
First, let us introduce the following short-hand notation, To simplify the equation of motion, we assume that the gravitational fields and their derivatives are small, {λ g,f , ν g,f } 1 and {rλ g,f , rν g,f } 1, but we keep all nonlinearities in µ. This results in the following set of Einstein equation for g µν , With the same approximations, the Einstein equations of f µν can be arranged to Note that we do not present the φφ-components of the Einstein equations as these coincide with the θθ-components. The last ingredient is the Bianchi constraint, that simplifies to Even in this simplified approximation, the equations are not solvable analytically. To find analytic solutions, we will go to large radii where all metric functions and their derivatives have to be small because we assume spacetime to be asymptotically de Sitter while remaining far inside the de Sitter horizon 3/ √ Λ. The other limit is for small radii compared to the Compton wavelength of the massive spin-2 field. First, we will solve the equations of motion for r r S with r S the Schwarzschild radius of a localised source to be defined later. For large radii smaller than the de Sitter horizon, all metric functions, including µ, are assumed to be small such that both metrics approach de Sitter solution asymptotically. We expand the Einstein equations and the Bianchi constraint further in µ 1 and rµ 1. It is convenient to introduce the auxiliary functions, In terms of these functions, the linearised Bianchi constraint (B.5) reads, which is independent of the +-fields. We can find linear combinations of the Einstein equations for g µν and f µν such that the +-fields drop out of the expression, where we have used the Bianchi constrained already to eliminate ν − and ν − in terms of λ − and λ − . The equations reduce to two coupled differential equations for λ − and µ, while ν − is algebraically determined by the other two fields. These equations are solved by where C 1 is a constant of integration to be determined later. We already fixed the other constant of integration by requiring that the fields vanish in the limit r → ∞ in order to ensure that the metrics are asymptotically bidiagonal and de Sitter. Next, we turn to the other set of linearly independent equations that contain also the +-fields. These can be arranged to 0 = (ᾱ 2 − 1)(λ − + rλ − ) + (1 +ᾱ 2 )(λ + + rλ + ) − 2(1 +ᾱ 2 )Λr 2 , 0 = (ᾱ 2 − 1)λ − − (1 +ᾱ 2 )(λ + − rν + ) + 2(1 +ᾱ 2 )Λr 2 , Upon using eq. (B.10), these equations are solved by where C 2 is another constant of integration. Again, we have already used that the functions have to be bidiagonal and de Sitter in the limit r → ∞ to fix one of the constants of integration. Having solved the differential equations for all the fields, we can present the final solutions for the metric fields in the linearised limit. They are given by (1 +ᾱ 2 )r , (B.13) In the limit of a vanishing cosmological constant, Λ = 0, our linearised solutions reduce the the results obtained in Refs. [90,93]. Let us comment on the region of validity of these solutions. As we will see later, C 1 ∼ r S . Hence, all metric fields are small on scales r r S and inside the de Sitter horizon, r Λ −1/2 . This is not true however for the Stückelberg field µ, which is small only on scales r (r S /m 2 FP ) 1/3 . In general, this scale is larger than the Schwarzschild radius of the source. These solutions are hence not appropriate to describe the scales r S r (r S /m 2 FP ) 1/3 . In the next section, we seek solutions that are valid on these scales.

B.1.2 Compton regime
Next, we want to solve the equations on scales much smaller than the Compton wavelength of the massive spin-2 field, i.e. for r m −1 FP . Technically that means that we omit the terms m 2 FP r 2 · {λ g,f , ν g,f } 1 in the modified Einstein equations. Under this approximation we can integrate the tt-component of eq. (B.3) to where we defined the Schwarzschild radius of a source of mass M and radius R * . The solution is valid outside the compact object and we fixed the constants of integration by demanding regularity at the surface of the compact object and at the origin r = 0. Upon using the solution for λ g , the rr-component of eq. (B.3) simplifies to The tt-component of eq. (B.4) can easily be integrated to Here, we fixed the constant of integration by demanding regularity at the origin r = 0. With this result, the rr-component of eq. (B.3) simplifies to We have obtained expression for all the fields that show up in the Bianchi constraint in terms of µ.

(B.19)
In here, we defined the Vainshtein radius that we already encountered before. Our result coincides with the one obtained in Ref. [90]. Eq. (B.19) represents a seventh order algebraic polynomial for µ and as such has up to seven solutions. As argued in Ref. [93] for γ > 0, three solutions are real-valued for all r, while two solutions are real only up to a certain critical radius that depends on the other parameters. The last two solutions are complex-valued. The three everywhere real-valued solutions are shown in Fig. 12. In the opposite case γ < 0, two of the real-valued solutions become complex-valued and the remaining realvalued solutions are neither asymptotically bidiagonal nor restore GR inside the Vainshtein regime. The case γ = 0 is special because the polynomial is only of degree five. There is one solution that realises the Vainshtein mechanism for small radii, which however is not asymptotically bidiagonal and it is complex-valued for very small radii. Although there might be special situations in which the Vainshtein mechanism works even for this case [93,160], we exclude γ = 0 from our analysis and restrict to γ > 0.
We can infer further information about the general behavior of µ. For radii much larger than the Vainshtein radius, the term (r V /r) 3 on the right hand side of eq. (B.19) is small, which implies that the Stückelberg field is small as well, µ 1. This gives rise to three asymptotic values for µ as can be seen in Fig. 12. One of these is given by µ = 0 such that we can neglect nonlinearities in µ. This solution corresponds to the one obtained when linearising the field equations as in appendix B.1.1. For small radii, the term (r V /r) 3 on the right hand side is large. This allows to identify two classes of solutions. Either, µ on the right hand side compensates that large term. This is the case for µ = −1 or µ = ±1/ √ γ. Note that the second solution exists only for γ > 0. These are the three asymptotic values that show up in Fig. 12. The other possibility is that µ diverges as r → 0. These solutions coincide, if we choose the constants of integration to be (B.23)

B.2 Existence of Vainshtein-Yukawa branch
So far, we have studied analytic solutions on two different scales. For large radii, r r V , the gravitational potential is a combination of the 1/r-law and a Yukawa-type potential: as follows from eqs. (B.13) and (B.23). On smaller length scales, r S r r V the Stückelberg field µ is nonlinear. The nonlinearities are such that the gravitational potential is given by ν g = − r S r (B.25) Figure 13. Exclusion plots for a working Vainshtein mechanism for the β1β2β3-model. Left panel: In the blue-shaded region β < d1/d2 is implied while in the yellow-shaded region d2 < 0 holds. In the overlap, there is no everywhere real solution for the Stückelberg field µ that incorporates the Vainshtein mechanism at small radii while giving rise to the correct asymptotics. Right panel: In the green-shaded region µ 1 is no solution for large radii while in the purple-shaded region γ < 1. The yellow-shaded region shows the overlap of the bounds presented in the left panel. In the gray shaded region the Higuchi bound is violated. In summary, only outside the shaded regions, the bimetric model incorporates a working Vainshtein mechanism.
as follows from eq. (B.16) on a solution where µ is constant. This is summarised in eq. (4.4). For the screening mechanism to work it is necessary that the solution of the Stückelberg field realising these two regimes exists and is real-valued for every r without branch cuts.
This point was studied in detail in [90]. Firstly, the nonlinear equation (B.19) determining µ must give rise to a solution with µ 1 for r r V , which matches the linearised solution. This is the case only if β < √ γ (B.26) is satisfied. This bound is stated in eq. (4.6) and ensures a consistent asymptotic behavior of the nonlinear solutions. This is the condition ensuring consistent asymptotics as stated in eq. (4.6). Secondly, the linearised solution can be smoothly connected only to the nonlinear solution with µ = −1/ √ γ, which is required to be larger than the solution with µ = −1. Therefore, a consistent screening regime requires the parameters to satisfy γ > 1 (B.27) as stated in eq. (4.5). Thirdly, there must be a solution that is real-valued for every r and interpolates between the two aforementioned asymptotic limits without branch cuts. This is the case if the parameters satisfy the following bound: where d 1 = −1 + 6(1 +ᾱ 2 ) √ γ(1 + γ) − (13 + 12ᾱ 2 )γ , d 2 = 1 + 3ᾱ 2 − 6(1 +ᾱ 2 ) √ γ + 3(1 +ᾱ 2 )γ . (B.29) The bound in eq. (B.27) implies β 3 > 0. Combining eqs. (B.26) and (B.27) then implies β 2 < 0. A theoretically consistent expansion history requires β 1 > 0. Therefore, the β 1 β 2 β 3 -model is the only bimetric model with up to three free parameters that allows for Vainshtein screening and has a consistent background cosmology. In fig. 13 the bounds in eqs. (B.26) to (B.28) are presented in terms of the physical parametersᾱ and m 2 FP /Λ for that model. To express the interaction in terms of the physical parameters, we use the parameter relations reported in appendix A. In the left panel the violation of the bound (B.28) is shown. In the overlap of the blue-and yellow-shaded regions, the Vainshtein-Yukawa solution does not exist. In the right panel, we combine this bound with the other two. The shaded regions represent, where the bounds are violated and hence where the Vainshtein screening mechanism does not work. Only the region left white gives rise to screening. In fig. 7, the Vainshtein bounds are combined with the bounds ensuring a consistent background cosmology.