Stable cosmology in chameleon bigravity

The recently proposed chameleonic extension of bigravity theory, by including a scalar field dependence in the graviton potential, avoids several fine-tunings found to be necessary in usual massive bigravity. In particular it ensures that the Higuchi bound is satisfied at all scales, that no Vainshtein mechanism is needed to satisfy solar-system experiments, and that the strong coupling scale is always above the scale of cosmological interest all the way up to the early universe. This paper extends the previous work by presenting a stable example of cosmology in the chameleon bigravity model. We find a set of initial conditions and parameters such that the derived stability conditions on general flat Friedmann background are satisfied at all times. The evolution goes through radiation dominated, matter dominated, and de Sitter eras. We argue that the parameter space allowing for such a stable evolution may be large enough to encompass an observationally viable evolution. We also argue that our model satisfies all known constraints due to gravitational wave observations so far and thus can be considered as a unique testing ground of gravitational wave phenomenologies in bimetric theories of gravity.

The recently proposed chameleonic extension of bigravity theory, by including a scalar field dependence in the graviton potential, avoids several fine-tunings found to be necessary in usual massive bigravity. In particular it ensures that the Higuchi bound is satisfied at all scales, that no Vainshtein mechanism is needed to satisfy Solar System experiments, and that the strong coupling scale is always above the scale of cosmological interest all the way up to the early Universe. This paper extends the previous work by presenting a stable example of cosmology in the chameleon bigravity model. We find a set of initial conditions and parameters such that the derived stability conditions on general flat Friedmann background are satisfied at all times. The evolution goes through radiation-dominated, matter-dominated, and de Sitter eras. We argue that the parameter space allowing for such a stable evolution may be large enough to encompass an observationally viable evolution. We also argue that our model satisfies all known constraints due to gravitational wave observations so far and thus can be considered as a unique testing ground of gravitational wave phenomenologies in bimetric theories of gravity.

I. INTRODUCTION
Bimetric theories are an intensively studied class of massive gravity theories considered as an alternative to general relativity (GR). On one hand they predict new phenomena, such as the graviton oscillation [1,2]. On the other hand, bimetric theories contain both a massless and a massive spin-2 field. It has been nontrivial to construct a consistent theory of massive gravity. The first bimetric model free of the so-called Boulware-Deser ghost was proposed by Hassan and Rosen [3], based on the de Rham-Gabadadze-Tolley (dRGT) ghost-free massive gravity model [4].
The bigravity [3], although allowing for a stable cosmological evolution, still requires an important fine-tuning of its parameters in order to be consistent. On one hand, it has been shown that to accommodate a stable evolution, the mass parameter m (controlling the graviton potential terms) needs to be generically much larger than today's Hubble parameter, i.e. m H 0 [5,6]. This condition forbids the graviton mass to account for the accelerated expansion of the Universe today. On the other hand, one needs another fine-tuning for (i) the Vainshtein mechanism [7] to effectively screen extra forces on Solar System scales, for (ii) letting the theory be differentiable from GR by leaving nontrivial phenomenology, while (iii) satisfying the Higuchi bound m T > O(1)H 0 [8], where m T is the mass of the tensor modes (proportional but not equal to m). Finally, the strong coupling is encountered at a rather low scale Λ 3 = (M Pl m 2 ) 1/3 easily by going early enough in the history of the Universe, which makes the need for a (partial) UV completion all the more important.
In response to these practical issues, it has been recently proposed to add a new chameleonlike degree of freedom to the theory [9]. In this model, the constant coefficients appearing in the graviton potential are promoted to be general functions of the new scalar field φ, and matter is coupled to gravity through a φ-dependent effective metric. In this way, the effective graviton mass m T becomes environment dependent, so that m 2 T scales as the local energy density of matter ρ. This mechanism allows us to evade the need for the Vainshtein mechanism to screen the extra gravitational forces on Solar System scales, and lets the theory be viable against strong coupling, Higuchi bound, and instabilities up to the very early Universe. The scalar field also has a high enough mass to be possibly not detectable by fifth-force experiments [9]. A possible cosmological application of the chameleonic extension of bigravity theory has been studied in [10].
In this work we study further the model presented in Ref. [9]. Indeed, notwithstanding the arguments in favor of stability and wider applicability that were given, it is important to study the compatibility of the theory versus the observed cosmic evolution. First, we present a detailed study of the scaling solutions, including the conditions for stability under homogeneous perturbations. Second, we present the stability conditions derived from studying the action that is quadratic in inhomogeneous linear perturbations around a flat Friedmann-Lemaître-Robertson-Walker (FLRW) spacetime. Finally, we present a viable set of parameters and initial conditions that upon numerical integration leads to a stable cosmological evolution, including radiation-dominated, matter-dominated, and de Sitter phases.
The text is organized as follows. In Sec. II we review the chameleon bigravity model presented in Ref. [9], defining the action and the background equations obtained from its variation. In Sec. III we present the scaling solutions of the model, and their respective stability under homogeneous perturbations. In Sec. IV we discuss inhomogeneous linear perturbations of the model, and the derivation of the stability conditions in a general flat FLRW universe. In Sec. V we present the numerical integration, as well as the chosen parameters and initial conditions. Finally, we conclude in Sec. VI and briefly present future extensions of this work.

A. Action
The chameleon bigravity model is defined by the total action S tot = S EH + S m + S φ + S mat [9]. In this model, the usual ghost-free bimetric theory is supplemented by a scalar field φ, coupled to both metrics via the promotion of the coefficients found in the graviton potential into the functions β i (φ). The gravitational part of the action is given explicitly by where M g and M f stand for the respective bare Planck masses of the gravitational g and f sectors. We also define κ ≡ M 2 f /M 2 g for later convenience. Just as in the usual bigravity case the construction of the potentials U i relies on powers of the metric square root s α β ≡ ( g −1 f ) α β such that s α γ s γ β = g αδ f δβ . By defining T n ≡ Tr[s n ], we have The potentials U 0 and U 4 constitute the two cosmological constants of the metric sectors g and f , respectively. The terms β i (φ)U i also play the role of potentials for the field φ. Finally, to implement the chameleon mechanism, the matter sector is coupled nonminimally to the metric g µν , i.e.
where ψ stands for the different matter fields,g µν = A 2 (φ)g µν , and A(φ) is a universal coupling function. In order to simplify the treatment, we adopt the choice of general functions A(φ) and β i (φ), following Ref. [9]. We thus set with i ∈ {0, · · · , 4}. These choices are sufficient to obtain a scaling solution described in Sec. III. We will use these specific functions for our numerical work.

B. Background equations
In order to study cosmological backgrounds, we choose a flat FLRW ansatz for both metrics, i.e.
Under these assumptions, the computation of the metric square root s µ ν becomes much simpler. We further define the Hubble parameters associated with each gravitational sector, H ≡ȧ/a and H f ≡ (aξ)˙/(acξ 2 ), where the dot stands for a derivative with respect to the cosmic time t. On such a FLRW background, the equations of motion become the two Friedmann equations (with R and U defined below) as well as the two dynamical equations (with J defined below) and the equation of motion for the chameleon scalar field (with Q defined below). In these equations we have used and ρ and P are, respectively, the total energy density and pressure of the matter fields. By combining the Friedmann In order to represent perfect fluids in the latter analysis, one can choose, for instance, to use k-essence scalar fields, where X α ≡ − 1 2g µν ∂ µ ψ α ∂ ν ψ α is the canonical kinetic term for a scalar field ψ α . One can then identify pressure P α , energy density ρ α , and the sound speed squared c 2 s,α in the Jordan frame as

A. Scaling solutions
It is possible to find exact and approximate scaling solutions to Eqs. (8)- (12). We find that in radiation-and cosmological-constant-dominated eras there exist exact scaling solutions. In the matter-dominated era one can find an exact scaling solution only for β = 0. When 0 < β 1 this turns into an approximate scaling solution. For a radiation-dominated or de Sitter Universe, on the other hand, the exact scaling solutions persist for any value of β.
From the Friedmann equation (9) for f µν , we can show that both ξ = constant and c = constant in any scaling solution. Assuming a power law behavior of the scale factor, all terms in the Friedmann equations (8) and (9) should scale as t −2 . Then one can immediately see from the graviton potential terms that if ξ is constant, then where we have used the standard scaling of the scale factor a(t) ∼ t 2/n (with n = 4 for radiation domination and n = 3 for matter domination, here with β = 0) and introduced the e-folding number N e = ln (a(t)/a i ). Here, t i (> 0) is the initial time and a i = a(t = t i ). Denoting a derivative with respect to the e-folding time by a prime, one obtains In the case of an exponential increase of the scale factor, i.e. in a purely de Sitter or Λ-dominated universe, this last equation (18) can be extended with the value n = 0, since all background quantities (excepting the scale factor) can be taken as constant. Finally, we also have In a radiation-dominated universe (and in de Sitter) the scaling expressions presented above can be shown to satisfy all background equations trivially.
On the other hand, in a matter-dominated universe, once we adopt the choices in Eq. (6), we combine background equations to find the following condition including β: As c and κ are positive, this condition with β = 0 can be satisfied only if λ ≤ √ 3. Since we are interested in the regime λ β to have m 2 T ∝ ρ [9], the condition (20) implies that there is no exact scaling solution in a matter-dominated era unless β = 0. However, if β is not zero but small enough then the system with λ β exhibits an approximate scaling behavior. Therefore, we impose that β ≈ 0 to allow for an approximate scaling solution.

B. Stability under homogeneous perturbation of the scaling solutions
For practicality, the chameleon scalar field and the Hubble expansion rate are rendered dimensionless using mass parameters of the theory, i.e., The equations are then written in terms of ln h, ϕ, ξ, and c. Homogeneous perturbations of the fields are defined as where is a small expansion parameter, h 0 is the initial background value of h, andξ and c (0) are the constant values of ξ and c, respectively, for the scaling solutions. The background equations are then expanded to first order in . After using the zeroth order equations of motion to set, for instance, c 0 , κ, c 4 , and the initial amount of matter (either radiation or dust) in terms of c (0) and the other background variables, one can solve the linearized equations for the variables h (1) , ξ (1) , and c (1) in terms of ϕ (1) and its derivatives.
Upon making these replacements, one finds the dynamics is uniquely determined by a second-order equation for ϕ (1) . This can be written as during radiation domination (with general β), and during matter domination (with β = 0), where with i = r, m. To guarantee the stability during radiation and matter dominations, respectively, it is necessary and sufficient that Here, it is understood thatξ r andξ m in (27) and (28) are the constant values of ξ in radiation-and matter-dominated epochs, respectively.

IV. STABILITY CONDITIONS OF PERTURBATIONS
One can define the perturbations of the fields with respect to the spatially flat FLRW background as follows. In the Arnowitt-Deser-Misner (ADM) decomposition, the (perturbed) metrics are written as One can then decompose the lapses, shifts, and 3D metrics separately as where Φ,Φ, δN i , δÑ i , δγ ij , and δγ ij are the perturbations. In particular, we are free to choose N = 1 by the time reparametrization invariance, and we also have that N i =Ñ i = 0 in our particular background. One may use other equivalent definitions of the perturbations; for instance, as long as the background equations of motion are taken into account, any definitions that differ only at second order will be equivalent as far as the quadratic action is concerned. Finally, as perturbations are studied only linearly and on a spatially homogeneous and isotropic background, one can decompose the perturbations of the shifts and 3D metrics into SO(3) scalar, vector, and tensor representations, i.e.
The Laplacian is defined as ∆ ≡ δ kl ∂ k ∂ l , and we use the notation O (ij) ≡ 1 2 (O ij + O ji ) to denote symmetrization of the indices. The latin indices of partial derivatives and perturbations can be raised and lowered with δ ij and δ ij . The perturbations of the chameleon scalar field and matter fields are The full action is then expanded to second order in the linear perturbations just defined. In particular the perturbations to the metric square root can be computed along the lines of [11]. The treatment is separated into tensor, vector, and scalar sectors. For later use, we choose to represent the matter content of the Universe by two perfect fluids, thus labeled by ψ α , with α ∈ {1, 2}.

A. Tensor perturbations
The quadratic action for tensor perturbations (written in Fourier space) reduces to is the comoving momentum of a perturbation mode, and In obtaining this form, we have used both Friedmann equations. One obtains a simple no-ghost condition from the tensor sector, i.e.
The squared sound speeds of the tensor modes are c 2 T,1 = 1 and c 2 T,2 = c 2 for h ij andh ij , respectively. Due to the time dependence of the background geometry, the graviton mass cannot be defined without ambiguities of order O(H) in general. On the other hand, in de Sitter spacetime with ξ = constant and c = 1, it is the combinations h − ij and h + ij = h ij + κξ 2h ij that are the two eigenmodes of the mass matrix. In such a case, one can simply rewrite the Fourier space action in the form In this case m T is the mass of the massive mode, and both graviton sound speeds are equal to unity.

B. Vector perturbations
After integrating out two nondynamical vectorial degrees of freedom (e.g. B i andB i ), the quadratic action for vector perturbations reduces to (in Fourier space) where E − i = E i −Ẽ i is the only propagating (massive) vector mode, and The associated no-ghost condition in the UV regime is, using c > 0 and ξ > 0, The no-gradient-instability condition, c 2 V ≥ 0, implies C. Scalar perturbations The study of the quadratic action for the scalar perturbations requires more work than the vector and tensor sectors. Because of the size of the expressions, we do not give here the full Lagrangian. Instead, we give here the no-ghost conditions, which must be satisfied at all times during the numerical integration, and the squared sound speeds of the scalar sector, which must be positive at all times.
We start by integrating out four nondynamical degrees of freedom that enforce the Hamiltonian and (longitudinal part of) the momentum constraints (i.e., Φ,Φ, B,B). One can integrate out as well the would-be Boulware-Deser (BD) ghost, which is rendered nondynamical by the particular structure of the graviton potential term. One can further use the remaining gauge freedom to set, for instance, the spatially flat gauge, Ψ = E = 0. Eventually, one finds that in addition to the two matter perturbation modes, one has two scalar degrees of freedom, one from the chameleon scalar and the other from the massive graviton.
In order to find both no-ghost conditions and dispersion relations, we take the subhorizon limit k aH. Indeed, we are solely interested in checking the presence or absence of instabilities in the UV, any IR instability being less problematic [12].

No-ghost conditions
In the subhorizon limit, the action can be written schematically as where K = K, F = −F, M = M are 4 × 4 real matrices, and Y is a vector containing the four remaining dynamical scalar perturbations, each of which may or may not be rescaled by a positive constant coefficient. The kinetic matrix K can then be diagonalized, yielding the eigenvalues up to overall positive constant coefficients. Because of some field redefinition used to diagonalize the kinetic matrix, the indices in κ i are arbitrary, but roughly correspond to, respectively, the scalar graviton, the chameleon field, and both matter perturbations. While κ 2 ≥ 0 is trivial and κ 3 , κ 4 ≥ 0 translate into the null-energy conditions on matter fields, i.e., ρ α + P α ≥ 0 (where α is an index designing a specific matter field), κ 1 ≥ 0 yields a nontrivial no-ghost condition which will be checked at all times during the numerical integration. We also want to monitor the scalar sound speeds squared, which are read off from the dispersion relations in the subhorizon limit.

Scalar sound speeds
The scalar sound speed for high frequency modes can be found by studying the dispersion relation in the subhorizon limit. Two modes propagate with the usual squared sound speeds c 2 s,α of perfect fluids and can thus be identified with the matter modes. The product and the sum of the two remaining scalar sound speeds squared, c 2 s,1 , c 2 s,2 , are given by and where ρ = ρ 1 + ρ 2 and ρ = P 1 + P 2 . If one considers the vector sector no-ghost condition, J > 0, then Σ 2 < 0. The scalar sound speeds squared provide new stability conditions, as these need to be real and positive. We thus require that Although we do not give here the analytical expressions for the single squared sound speeds, which would be too large to write, we obtain their numerical value in the next section as part of our numerical example cosmology (see Fig. 4). The reader may find a discussion on the respective contributions of the chameleon and the scalar graviton to the scalar squared sound speeds in Appendix A.

A. Set of equations
Although in principle one can obtain several background equations -e.g. both Friedmann equations, both second Einstein equations, the scalar equation of motion, or the combination Eq. (14),-not all the equations will be directly integrated. For instance, this last equation can be used to fix the fiducial function c. Similarly, both Friedmann equations can be used to set two parameters or integration constants, as will be shown below. Of the equations cited above, only three will remain to be integrated: both second Einstein equations and the scalar equation of motion. In addition to finding the right set of equations, the choice of adequate initial conditions (ICs) is also essential. In what follows, a subscript i stands for the quantity evaluated at initial time.
Although in the previous section we were able to derive the results while keeping the functions β j (φ), j ∈ {0, · · · , 4}, and A(φ) completely general, these need to be specified for the sake of numerical integration. We will thus from now on use the example model defined in (6).
Several definitions help render the equations more practical for the purpose of numerical integration. First of all we consider the equations of motion in e-fold time with its initial value being N e,i = 0. We then define dimensionless variables. We start by using the dimensionless chameleon scalar field, ϕ, and Hubble parameter, h, as defined in Eq. (21). For the matter energy densities, we split the energy density of the matter fields (in the Jordan frame, for which a JF = A a) as where the subscripts r, d, and Λ, indicate the radiation, dust, and cosmological constant, respectively. We then define where r r , r m , and r Λ are dimensionless and constant throughout the evolution. Using these definitions, the Friedmann equation for the physical metric becomes 3h 2 = 1 2 h 2 ϕ 2 + e −λϕ c 0 + 3c 1 ξ + 3c 2 ξ 2 + c 3 ξ 3 + e β ϕ r d e −3Ne + r r e −4Ne + e 4β ϕ r l , while the Friedmann equation for the fiducial metric can be written where as noted previously a prime denotes differentiation with respect to N , and we have defined It is instructive to rewrite the physical Friedmann equation as For this, we have defined the Einstein frame density parameters (60) The new subscripts k and V indicate contributions from the chameleon kinetic energy and from the graviton potential term, respectively. We can also define the Jordan frame density parameters, using the fact that where η is the conformal time defined by η ≡ t 0 dt a(t ) . This allows us to write In a similar way, Therefore we can replace r r , r d , r Λ with either Jordan frame or Einstein frame density parameters, evaluated at initial time, i.e., In terms of the new variables we have that Eq. (14) can be rewritten as which defines c in terms of the other dynamical variables. When using this definition in the fiducial second Einstein equation, this reduces the degree of the equation to 1, with respect to the variable of interest ξ.
The set of dynamical equations to be integrated, the two second Einstein equations and the chameleon field equation, can be written as and, because of the choice of the variables/parameters, they do not explicitly depend on any scale, e.g., M g or m.

B. Requirements on initial data
Using a rescaling of the constants one can, without loss of generality, set the values of the ICs ϕ i , ξ i , and h i . In detail, this can be done, for example, by (i) redefining m 2 to set ϕ i = 0, (ii) redefining the c j and M f to set ξ i = 1, and (iii) an additional overall rescaling of the constants c j , which we use to set h i = 1. Once this is done, we only need to give one supplementary IC, i.e. ϕ i . Then the total set of yet needed ICs and parameters is c 0 , c 1 , c 2 , c 3 , c 4 , r r , r d , r Λ , λ, κ, β, ϕ i .
We can use the two Friedmann equations to set two of the parameters (or ICs, in principle). Without loss of generality, we solve them in terms of c 0 and κ (by linear equations). The initial conditions for the integration are set in a radiation-domination epoch, with the Universe obeying a scaling solution. These initial conditions allow us to recover a cosmology accommodating our Universe. In order to start with a radiation-domination phase, one simply needs to set 0 < 1 − Ω JF r,i 1. Since we also want to start from a scaling behavior during radiation domination, the remaining ICs and parameters are imposed so that the dynamics of the scale factor and the scalar field satisfy the scaling solution values found in Sec. III A, i.e., We choose to replace the parameters c 1 , c 2 , c 3 , and c 4 with new, more practical and transparent parameters. First, two of the constants can be chosen so that the condition (40) is always satisfied. This can be done, for example, by letting where both A and B are positive constants (and new parameters that replace two among c 1 , c 2 , and c 3 ), which is sufficient to guarantee that J > 0 for any ξ. Second, one may use Eq. (14), while approximating ξ i ≈ 0, to set c at the initial time to a specific value instead of one of the c i 's. Finally, we use the expression of the vector squared sound speed to replace the last parameter.

C. Results
Based on the previous section, we describe here a set of parameters which allows for an evolution similar to the usual Λ cold dark matter models. The values used in our example are , where the subscripts "in" or "i" mean the respective initial value. The initial density parameter for radiation, Ω EF ri , is directly determined by the Friedmann equation (59) at initial time, and all other parameters are fully determined by this set of choices. The only fine-tuned value is Ω EF Λi , which we have chosen in order to have Ω Λ of order unity today. In practice, this is the same as the cosmological constant problem today.
The simple choice of parameters in Eq. (71) is meant to show that it is possible to obtain a realistic cosmological evolution. It does not recover exactly today's observed values. However, it is possible, by an appropriate choice of constants -and without fine-tuning anything other than the cosmological constant -to obtain an evolution fitting more closely to data; e.g. one can reproduce today's abundances and other data. This, along with the constraints on the model from today's observational data, will be studied further in a future work.
For the sake of exposition, we present the evolution 1 of the density parameters in Fig. 1, while the evolution of other relevant variables is presented in Fig. 2. The evolution, starting from a radiation-dominated era, moves on to a matter-dominated era, finally attaining a final de Sitter phase. Given our set of initial density parameters, the system stays O(10) e-folds in each era before settling to a de Sitter epoch (roughly from 0 to 12 e-folds for radiation domination, from 12 to 19 e-folds for matter domination, from 19 to the end for the de Sitter era). However, by arranging these density parameters, one can achieve very different numbers of e-folds spent in each era. In order to have a handle on the precision of the numerical integration, we check all along the evolution to which extent the Friedmann equations are satisfied. For this purpose, one may define, for instance, inspired by both Friedmann equations (59) and (57), and where α stands for any of the species, i.e., indices {r, d, Λ, k, V }. The evolution of these two constraints is presented in Fig. 3. Both constraints are seen to be satisfied up to order O(10 −5 ). In our implementation, this has been achieved by using constraint damping, i.e., adding the constraint equations into the dynamical equations of motion (after normalizing the constraints by an appropriate factor), in order to damp any unwanted deviation. In addition to the Friedmann equations, we also present in Fig. 4 the evolution of the sound speeds and the fiducial lapse c. Together with the no-ghost conditions, which are found to be satisfied all along the evolution, the positivity of these shows that the background is stable under cosmological perturbations.
Finally, in order to demonstrate the purpose of the new scaling brought by the scalar field dependence in the graviton potential, we plot the Higuchi condition along the evolution in Fig. 5. The generalized Higuchi bound is seen to be well satisfied during the three eras, and in particular both at late and early times.

VI. DISCUSSION AND CONCLUSIONS
Following the recent proposal [9] of an extended massive bigravity theory supplemented by a chameleon scalar field as a means to cure or evade the fine-tunings of the original theory and improve its applicability, we have found it important to study further its validity and implications. For this reason, in this work we have explored the stability conditions of the model and confirmed its intended behavior by integrating numerically the equations of motion. In particular, we have numerically confirmed that at all times, the Higuchi ghost is never present: indeed, the presence of this ghost represented one of the most serious problems for a viable phenomenology of the original bigravity theory. In our model though, we have here shown that if no Higuchi ghost is present at one scale then the same ghost will not appear during the whole evolution of the Universe including the early epoch. This set of such allowed initial conditions is not of zero measure in general, so that we do not need to fine-tune the parameters of the theory.
The study of the action quadratic in perturbations with respect to a general flat FLRW background leads, in the UV, to no-ghost conditions for the tensor, vector, and scalar sectors. In addition to this, we have found the explicit action for the tensor and vector linear perturbations and for the scalar linear perturbations in the UV. From these the propagation speeds at short scales are easily extracted, thus leading to additional no-instability conditions. It is  found as expected that the theory propagates four tensor, two vector, and two scalar degrees of freedom (not including matter degrees of freedom), thus corresponding to the expected massive spin-2, massless spin-2, and chameleon scalar of the theory. In order to show the typical background time evolution, we have numerically integrated the background equations by using a choice of initial parameters consistent with an initial radiation-dominated era of the Universe. As supplementary input for the initial conditions, we have required that the stability conditions be satisfied and that the parameters of the theory are in the regime of interest for the expected scaling behaviors. The evolution displays an initial radiation-dominated era, followed by matter domination and a de Sitter era. The no-instability conditions are satisfied all along the evolution, and, in our implementation, the constraint equations show a numerical error of order O(10 −5 ) at most. This stable evolution comforts us into arguing that it may be possible to find a region of the parameter space allowing a close match with our cosmological observations. The recent binary neutron star merger observation, the first gravitational and electromagnetic wave multimessenger detection [13], has allowed us to set stringent bounds on the speed difference between gravitational and electromagnetic waves (see, e.g., [14][15][16][17]). Although in our model one of the gravitons propagates with a slightly modified sound speed c (see the lower panel of Fig. 4), the physical metric remains unaffected and the interactions between the two metrics are suppressed by the smallness of m 2 β i , i ∈ {1, 2, 3}. This implies that the propagation of gravitational waves in our model is essentially the same as that of photons as far as m 2 β i are small enough compared to the typical (squared) energy scales of the gravitational waves produced astrophysically. As a result, the constraint on our model from GW170817 is essentially the same as those from the previous GW observations (e.g., [18]) [16]. Concretely, the constraint is of the form of an upper bound on the mass of the graviton (which was not improved by GW170817) of m T < 1.2 × 10 −22 eV. While this bound has to and can be satisfied today, the scalar field dependence of the graviton mass in our model allows without problem for a larger mass at early times, rendering the cosmological evolution stable all the time. Therefore, our model can be considered as a unique testing ground of gravitational wave phenomenologies in bimetric theories of gravity. For example, it is intriguing to investigate the possible modification of the waveform of the gravitational wave signal due to the influence of the massive graviton.
As a clear avenue for future extension, the evolution of cosmological perturbations and an improved understanding of the viable parameter space will be considered in a future work. Furthermore, it may be interesting to study the detailed working of the screening mechanism for the chameleon scalar field and scalar graviton modes.