Gravitational Wave Oscillations in Bigravity

TL;DR

This paper develops a theoretical framework for gravitational wave oscillations in bigravity, predicting observable effects similar to neutrino oscillations and providing new constraints on graviton properties.

Contribution

It derives consistent equations for gravitational wave oscillations in bigravity, highlighting the basis difference between wave propagation and detection, and discusses phenomenological implications.

Findings

Predicted oscillatory behavior of gravitational waves in bigravity.

Derived new limits on graviton mass and coupling parameters.

Highlighted potential observational signatures in gravitational wave data.

Abstract

We derive consistent equations for gravitational wave oscillations in bigravity. In this framework a second dynamical tensor field is introduced in addition to General Relativity and coupled such that one massless and one massive linear combination arise. Only one of the two tensors is the physical metric coupling to matter, and thus the basis in which gravitational waves propagate is different from the basis where the wave is produced and detected. Therefore, one should expect -- in analogy to neutrino oscillations -- to observe an oscillatory behavior. We show how this behavior arises explicitly, discuss phenomenological implications and present new limits on the graviton parameter space in bigravity.

Figures

Equations

$$R_{\mu \nu }-\frac{1}{2}{g}_{\mu \nu }R+B_{\mu \nu }(g)=\frac{1}{M_g^2}{T}_{\mu \nu },$$

$${\tilde{R}}_{\mu \nu }-\frac{1}{2}{\tilde{g}}_{\mu \nu }\tilde{R}+{\tilde{B}}_{\mu \nu }(\tilde{g})=0,$$

$${\nabla }_{\mu }{B}_{\nu }^{\mu }=0,{\tilde{\nabla }}_{\mu }{\tilde{B}}_{\nu }^{\mu }=0,{\nabla }_{\mu }{T}_{\nu }^{\mu }=0,$$

$$\frac{3}{a^2}\left(H^2+k\right)=\Lambda (y)+\frac{\rho (\eta )}{M_g^2},$$

$$\frac{3}{b^2}\left(J^2/{\tilde{c}}^2+k\right)=\frac{\tilde{\rho }(y)}{M_{\tilde{g}}^2},$$

$$\begin{array}{c}{\beta }_1{\cos }^2\theta y^{-1}+(3{\beta }_2{\cos }^2\theta -{\beta }_0{\sin }^2\theta )+\\ +(3{\beta }_3{\cos }^2\theta -3{\beta }_1{\sin }^2\theta )y+\\ +({\beta }_4{\cos }^2\theta -3{\beta }_2{\sin }^2\theta )y^2-{\beta }_3{\sin }^2\theta y^3=\frac{\rho (\eta )}{M_g^2{m}^2}.\end{array}$$

$$\rho (\eta )={\rho }_0\{\begin{array}{ll}1& \text{if }\omega =-1,\\ {\left(\frac{a(\eta )}{a({\eta }_0)}\right)}^{-3(1+\omega )}& \text{if }\omega >-1,\end{array}$$

$$\delta y(\eta )=-\frac{\rho (\eta )}{3m^2{M}_g^2}\frac{y_{\ast }^3}{{\Gamma }_{\ast }({\cos }^2\theta +y_{\ast }^2{\sin }^2\theta )-2\frac{{\tilde{\rho }}_{\ast }{y}_{\ast }^4}{3m^2{M}_{\tilde{g}}^2}},$$

$$\tilde{c}(\eta )=$$

$$\simeq$$

$$h^{\prime \prime }+2H{h}^{\prime }+k^{2}h+{\sin }^2\theta m^2{\Gamma }_{\ast }{a}^2(h-\tilde{h})$$

$${\tilde{h}}^{\prime \prime }+2H{\tilde{h}}^{\prime }+k^2\tilde{h}+{\cos }^2\theta \frac{m^2{\Gamma }_{\ast }}{y_{\ast }^2}{a}^2(\tilde{h}-h)$$

$$h_1^{\prime \prime }+2H{h}_1^{\prime }+k^2{h}_1$$

$$h_2^{\prime \prime }+2H{h}_2^{\prime }+k^2{h}_2+a^2{m}_g^2{h}_2$$

$$h(t,k)=$$

$$\tilde{h}(t,k)=$$

$${⟨h^2(z,k)⟩}_{T_0\ll T\ll T_{\ast }}=\frac{{\cos }^4\theta }{{\left({\cos }^2\theta +y_{\ast }^2{\sin }^2\theta \right)}^2}\left[1+y_{\ast }^4{\tan }^4\theta +2y_{\ast }^2{\tan }^2\theta \cos \left(\frac{\delta \omega }{H}\log (1+z)\right)\right].$$

$$L_{\text{coh}}\sim 0.1\mathrm{s}\frac{2E^2}{m_g^2}={\left(\frac{10^{-22}\mathrm{eV}}{m_g}\right)}^2\mathrm{Gpc}.$$

Full text

We derive consistent equations for gravitational wave oscillations in bigravity. In this framework a second dynamical tensor field is introduced in addition to General Relativity (GR) and coupled such that one massless and one massive linear combination arise. Only one of the two tensors is the physical metric coupling to matter, and thus the basis in which gravitational waves propagate is different from the basis where the wave is produced and detected. Therefore, one should expect – in analogy to neutrino oscillations – to observe an oscillatory behavior. We show for the first time how this behavior arises explicitly, discuss phenomenological implications and present new limits on the graviton parameter space in bigravity.

I Introduction

The question whether a theory of a massless spin two particle can have a consistent massive extension was a longstanding open problem. The quest that lead to formulate this theory took place in the second half of the last century Fierz and Pauli (1939); Fierz (1939); van Dam and Veltman (1970); Zakharov (1970); Vainshtein (1972); Boulware and Deser (1972). Only recently it has been proven that a consistent framework of massive gravity exists and relies on the existence of multiple spin two fields with non-linear interactions de Rham and Gabadadze (2010); de Rham et al. (2011, 2012); Hassan and Rosen (2012a, 2011); Hassan et al. (2012); Comelli et al. (2012a); Deffayet et al. (2013a, b); de Rham et al. (2016).

In this letter, we study a setup with two dynamical spin two fields corresponding to two metrics Hassan and Rosen (2012b, c) known as bigravity. The coupling of the metrics to matter is a delicate problem and has been discussed in de Rham et al. (2014) as an arbitrary choice of coupling reintroduces inconsistencies. Demanding the absence of a ghost in the theory translates into an asymmetric coupling of the metrics to matter, and this asymmetry is at the core of the physical phenomenon we will discuss in this letter. The simplest choice of matter coupling which permits a ghost-free theory is minimal coupling of one metric tensor to matter, which we will call the physical metric and no coupling of the second metric tensor to matter. This second metric tensor is a reference or sterile metric which only interacts with the physical metric via the non-linear terms in the Lagrangian. This situation is analogous to the introduction of a sterile neutrino which carries no electroweak charges.

In this theory the gravitational interactions are mediated by two gravitons, one massless and one massive. Since the two are superpositions of the physical and the sterile metric, their effective coupling to matter is different and depends on the mixing angle between the metrics. This leads to an oscillation phenomenon, first mentioned in Berezhiani et al. (2007) in a theory of massive gravity and Hassan et al. (2013) in bigravity. In this work we will study the propagation of gravitational waves (GW) in this bi-metric theory which are produced in the ‘flavor basis’ at the source, namely only as perturbations of the physical metric. Describing the wave propagation we find a close analogy to neutrino oscillation described in the wave-packet formalism.

This phenomenon is presented for the first time in a consistent approach. Attempts have been previously made in De Felice et al. (2014a); Narikawa et al. (2015), however only in a specific setting, and leading to an unphysical result; in particular, during propagation the mode coupling to matter exhibits an enhancement of the strain in violation of (local) energy conservation. We show that in the parameter space we consider physical no such behavior is found, as one should expect in a healthy theory. The novelty of our work in the bigravity setup is that we consider graviton masses corresponding to length scales which can be probed by astrophysical tests, while the majority of prior works have focused on much smaller graviton masses, i.e. of the order of the Hubble scale today. This approach makes it possible to confront direct detection data of GW signals as seen by the LIGO experiment Abbott et al. (2016a) with the oscillation hypothesis. The corresponding parameter space of $m_{g}=10^{-22}-10^{-20}$ eV and comparably large mixing angle $\theta$ is studied, in close resemblance to the effects of pure massive gravity. Note that previous studies have found instabilities which plague the parameter regime in which the graviton mass is of the Hubble scale today, and specific parameter choices are needed to obtain viable solutions De Felice et al. (2014b); Akrami et al. (2015). However, for the larger graviton masses probed here, this problem is considered to be less restricitve Babichev et al. (2016).

II Gravitational wave oscillations

In this model, the oscillation of metric perturbations is driven by classical dynamics of the Friedmann equations De Felice et al. (2014a). They are extracted from the Einstein field equations of bigavity von Strauss et al. (2012),

[TABLE]

with $B_{\mu\nu}(g)=m^{2}\cos^{2}(\theta)\sum_{n=0}^{3}\beta_{n}V^{(n)}_{\mu\nu}$ and $\tilde{B}_{\mu\nu}(\tilde{g})=m^{2}\sin^{2}(\theta)\sum_{n=1}^{4}\sqrt{g^{-1}\tilde{g}}\,\beta_{n}\tilde{V}^{(n)}_{\mu\nu}$, $\cos^{2}(\theta)=M_{\text{eff}}^{2}/M_{g}^{2}$, and $\sin^{2}(\theta)=M_{\text{eff}}^{2}/M_{\tilde{g}}^{2}$. The $V^{(n)}_{\mu\nu},\ \tilde{V}^{(n)}_{\mu\nu}$ encode the variation of the interaction terms in the action w.r.t. $g,\tilde{g}$. Furthermore, by applying the covariant derivatives to Eqs. (1), we obtain the conservation laws,

[TABLE]

the first two of which are known as Bianchi constraints.

II.1 Background cosmology

We now calculate the cosmological implications on a static background. For both $g$ and $\tilde{g}$, we assume an FRW background metric, $\mathrm{d}s^{2}=\,a(\eta)^{2}(-\mathrm{d}\eta^{2}+\mathrm{d}\vec{x}^{2})$ and $\mathrm{d}\tilde{s}^{2}=\,b(\eta)^{2}(-\tilde{c}(\eta)^{2}\,\mathrm{d}\eta^{2}+\nolinebreak\mathrm{d}\vec{x}^{2})\,$. The lapse function $\tilde{c}(\eta)$ determines the light cone for the second metric and plays a role for the propagation speed of the massive gravitational wave excitations. This is the most general ansatz compatible with a homogeneous and isotropic Universe Comelli et al. (2012b).

Plugging this ansatz into Eqs. (1) and omitting explicit dependencies yields the cosmic evolution equations,

[TABLE]

where $\Lambda(y)\equiv m^{2}\sin^{2}\theta\left[\beta_{0}+3\beta_{1}y+3\beta_{2}y^{2}+\beta_{3}y^{3}\right]$ and $\tilde{\rho}(y)\equiv M_{\tilde{g}}^{2}m^{2}\cos^{2}\theta\left[\beta_{1}y^{-3}+3\beta_{2}y^{-2}+3\beta_{3}y^{-1}+\beta_{4}\right]$. Here, a prime denotes a derivative w.r.t. $\eta$, $y=b/a$, and $H=a^{\prime}/a$ as well as $J=b^{\prime}/b$ are the Hubble parameters for both metrics in conformal time.

Moreover, Eqs. (2) imply that $\rho^{\prime}(\eta)=-3H(1+\omega)\rho(\eta)$ and $(\tilde{c}H-J)\left[\beta_{1}y+2\beta_{2}y^{2}+\beta_{3}y^{3}\right]\equiv(\tilde{c}H-J)\Gamma(y)=0,$ for a perfect fluid with equation of state $P=\omega\rho$. It was shown that only the vanishing of the round brackets yields a physically meaningful solution von Strauss et al. (2012). Thus, we find $J(\eta)=\tilde{c}(\eta)H(\eta)$.

Using this result, we can derive an algebraic equation for $y$, by subtracting Eq. (3a) from Eq. (3b),

[TABLE]

By assumption, $\rho$ is the density of a perfect fluid with $\omega\geq-1$, which behaves as von Strauss et al. (2012)

[TABLE]

such that any fluid of type $\omega>-1$ is diluted, i.e. $\rho\nobreak\to\nobreak 0$ for $\eta\to\nobreak\infty$. It is in fact sufficient to consider such densities, since any cosmological constant (CC) type of energy density may be included in the interaction terms of the bigravity theory.

In this limit, we denote the solution of Eq. (4) as $y_{*}$. An exact expression is in principle feasible, however not very enlightening. Therefore, and since we are interested in late times, we leave $y_{*}$ undetermined and linearize Eq. (4) as $y=y_{*}+\delta y$ to obtain,

[TABLE]

with the short-hand notation $\Gamma_{*}\,{=}\,\Gamma(y_{*})$ and $\tilde{\rho}_{*}\,{=}\,\tilde{\rho}(y_{*})$.

This manipulation allows us to rewrite Eq. (3a) as $a(\eta)^{-2}(H(\eta)^{2}+k)=\frac{1}{3}\Lambda_{*}+\frac{\rho(\eta)}{3M_{\text{Pl}}^{2}}$ with the physical CC $\Lambda_{*}=\Lambda(y_{*})$ and Planck mass, $M_{\text{Pl}}^{2}=\left(M_{g}^{2}\cos^{2}\theta+y_{*}^{2}\sin^{2}\theta-\frac{2\tilde{\rho}_{*}y_{*}^{4}}{3m^{2}M_{\tilde{g}}^{2}\Gamma_{*}}\right)\left(\cos^{2}\theta-\frac{2\tilde{\rho}_{*}y_{*}^{4}}{3m^{2}M_{\tilde{g}}^{2}\Gamma_{*}}\right)^{-1},$ which approaches $M_{g}^{2}\left(1+y_{*}^{2}\tan^{2}\theta\right)$, as $\tilde{\rho}_{*}\to 0$, in agreement with De Felice et al. (2014a); Narikawa et al. (2015).

Finally, we may use that $y^{\prime}=\left(\frac{b}{a}\right)^{\prime}=y(J-H)$ and $J=\tilde{c}\,H$ to find that,

[TABLE]

Note that $\tilde{c}$ can be both larger or less than 1, depending on the choice of the $\beta$-parameters. However, $\tilde{c}>1$ would introduce GWs propagating with a speed larger than the speed of light. In certain frameworks this might be acceptable, e.g. the present case is similar to the framework studied in Babichev et al. (2008), where all matter propagates on the $g$ background and no causal paradoxes arise.

From Eq. (II.1) we obtain $|\tilde{c}-1|\approx 10^{-20}$ for typical values in the parameter region of interest. This motivates the limit where $\tilde{c}=1$ and $y$ takes the constant value $y_{*}$, which we apply in the following.

II.2 Gravitational wave oscillations

We now address the propagation of tensor perturbations around the background metric.111The scalar mode couples to the trace of a conserved source and will thus in principle be excited, too. However, it is suppressed due to the Vainshtein effect de Rham (2014). Defining the transverse traceless components, the equations of motion are Comelli et al. (2012c),

[TABLE]

where $k=|\vec{k}|$ denotes the three-momentum and the polarization indices $+/\times$ are implicit. For the linear combinations $h_{1}\equiv\cos^{2}\theta\,h+\sin^{2}\theta\,y_{*}^{2}\,\tilde{h}\,\ \text{ and }\ h_{2}\equiv h-y_{*}^{2}\,\tilde{h}\,,$ one of the equations decouples and we obtain

[TABLE]

where we have defined the physical graviton mass $m_{g}^{2}\nolinebreak\equiv\nolinebreak m^{2}\Gamma_{*}(\sin^{2}\theta+\frac{\cos^{2}\theta}{y_{*}^{2}})$, and the source term is proportional to $\kappa(\theta,y_{*})\equiv(1-y_{*}^{2})/(\cos^{2}\theta+y_{*}^{2}\sin^{2}\theta)$. We observe that Eqs. (9) comprise one massless and one massive propagating tensor perturbation, where the latter is sourced by the former. Ignoring the Hubble rate, which is typically much smaller than the wave numbers $k$ under consideration, and setting the expansion rate to a constant, $a=1$, we can solve these equations and rotate back to the physical basis,

[TABLE]

where $\eta$ has been replaced by cosmic time $t$ as per $a=1$.

Since the graviton mass is restricted to be much smaller than the typical wave number $k$, we may expand $\sqrt{k^{2}+m_{g}^{2}}\simeq k\left[1+{m_{g}^{2}}/{(2k^{2})}\right]\equiv\omega_{0}+\delta\omega$. We see that the numerator in Eq. (10a) is minimized when the second cosine acquires a total phase shift of $\delta\omega\,T_{*}\,\pi$, and thus, $T_{*}(\omega_{0})=\frac{2\pi\,\omega_{0}}{m_{g}^{2}}\,,$ which coincides with the expression for the oscillation length for neutrinos, confirming our naïve expectation.

In order to make a quantitative statement about the modulation of the strain observed in GW observations, we average this expression over a timescale $T$, which is bigger than the period of one massless mode’s inverse frequency, $T_{0}=\frac{2\pi}{\omega_{0}}$, but much smaller than the period of the modulation induced by the mass term, $T_{*}=\frac{\pi}{\delta\omega}$. Squaring the strain, we find its envelope function where the normalization is determined by the condition $\left.\left\langle h^{2}(t,k)\right\rangle\right|_{T=0}=1$, i.e. initially a pure perturbation of the physical metric has been excited.

Finally, we aim to express the strain in terms of the cosmic redshift $z$, which is defined as $1+z=a(t_{0})/a(t)$. For a universe dominated by a CC, we find that $H=\mathrm{const.}$ and $a(t)=e^{Ht}$. We therefore express the time as $t=-\frac{1}{H}\log(1+z)$.222Note that we have reinstated $a(t)\neq$ const. in conflict with the condition $a=1$ used in the analytic derivation of Eq. (11). Thus Eq. (11) is only a valid approximation for small $z$. In summary, the squared amplitude of the GW signal in bigravity is modulated as

[TABLE]

At this point, we would like to point out that the phenomenon has previously been studied in De Felice et al. (2014a); Narikawa et al. (2015), where the authors find a modulation that is proportional to $\tilde{c}-1$. As we will outline in the following, this is not the leading effect in our analysis, where oscillations occur also in flat space. Furthermore, we find that the phenomenon leads to a reduction rather than amplification of the amplitude compared to GR, as expected from neutrino oscillations. Both are physically sensible outcomes.

II.3 Phenomenology

Given that we have reached a quantitative understanding of GW oscillations in terms of the modulation (11), we now ask whether this effect is visible in realistic scenarios. To this end, we have made use of the available data for the events GW150914 Abbott et al. (2016b) and GW151226 Abbott et al. (2016c) obtained by means of numerical simulations Wardell et al. (2016); Loffler et al. (2012); Pollney et al. (2011); Schnetter et al. (2004); Thornburg (2004); Ansorg et al. (2004); Dreyer et al. (2003); Goodale et al. (2003); Brown et al. (2009); Husa et al. (2006); Thomas and Schnetter (2010); SXS (2016). This yields the strain as it would be observed in a detector on Earth. We then modulate the strain according to Eq. (11). Two such examples for GW150914 are shown in Fig. 1, where the parameters are chosen such that one obtains a maximally visible effect, i.e. $\theta=\pi/4$ and $y_{*}=1$. One observes that a graviton mass of $m_{g}=10^{-22}\,\mathrm{eV}$ strongly changes the shape of the signal, where the modulation is at first strongly suppressing the amplitude and then gradually approaching the GR amplitude towards the typical merger peak, commonly referred to as chirp. On the other hand, a larger graviton mass $m_{g}=10^{-19}\,\mathrm{eV}$ leads to a global suppression of the amplitude by a constant factor. This effect is similar to the decoherence of oscillating neutrino wave packets and we will now briefly discuss this effect.

The massive and the massless modes propagate in wave packets with different group velocities $v_{g}=\frac{\partial\omega}{\partial k}$. As for very light, relativistic neutrinos, the difference of group velocities is approximately given by $\Delta v_{g}\simeq\frac{m_{g}^{2}}{2E^{2}}$. The wave packets will de-cohere, i.e. interference will be absent and the frequency dependence of the suppression is lost to a constant reduction, once the time of propagation exceeds $T_{\text{coh}}\sim L_{\text{coh}}/c\sim\sigma_{x}/\Delta v_{g}$, where $\sigma_{x}$ is the spatial/temporal width of the wave packet Beuthe (2003). Since its determination would involve an exact solution of the full set of Einstein equations for the system, it will be practically impossible to obtain $\sigma_{x}$. However, from the shape of the signal, we estimate $\sigma_{x}\sim 0.1\,\mathrm{s}$ for GW150914. Therefore, we find that for $E/\hbar\sim 100\,\mathrm{Hz}$,

[TABLE]

This rather heuristic argument is nevertheless in good agreement with Fig. 2, where for $m_{g}=10^{-22}\,\mathrm{eV}$ no averaging is observable at distances of the order 100 Mpc, while for $m_{g}=5\cdot 10^{-22}\,\mathrm{eV}$, or even $m_{g}=10^{-19}\,\mathrm{eV}$, the amplitude levels out for distances below the Gpc scale. Note that the longer time scale of GW151226 has little effect on the mass scale relevant for decoherence by virtue of Eq. (12). The resulting $\mathcal{O}(1)$ correction is not relevant for the estimate presented here.

Once the distance increases beyond the scale set by $L_{\text{coh}}$, the strain suppression relative to the prediction of GR caused by oscillations levels out. E.g. for $y_{*}=1$, $\theta=\pi/4$ we find $\left\langle h(t,k)\right\rangle_{T\gg T_{\text{coh}}}=\frac{2}{\pi}$, which predicts a suppression factor constant in frequency and distance of about 64 % at large redshifts, which is clearly confirmed in Fig. 2.

Note that higher graviton masses lead to shorter length scales before the amplitude averages out, in complete analogy to neutrino oscillations. In practice, such a frequency-independent suppression is indistinguishable from ordinary GWs of GR and would be interpreted as a larger redshift, i.e. one would generally overestimate the redshift on such BBH merger events. However if the source of the GW can be localized e.g. by electromagnetic observations, a discrepancy between the inferred redshift and the one obtained from the GW amplitude within GR could hint at graviton oscillations in the decoherence regime. Additionally, if a larger set of events becomes available, this can be used to constrain the larger-graviton mass regime by comparing expected distribution of BBH systems with the observed event rates.

For the low-mass regime, we can constrain the parameters of the model by demanding that the waveform be in agreement with the error bars of the observed events. We have used a simple $\chi^{2}$-analysis to obtain Fig. 3, where we set $y_{*}=1$ exploiting the parameter redundancy of $m$ and the $\beta_{i}$. For very small $m_{g}$, or $\theta\approx 0,\pi/2$, the suppression vanishes. Similarly, all events that lie beyond $L_{\text{coh}}$ are indistinguishable from an equivalent event in GR at larger $z$. From the remaining events the waveform in bigravity is clearly distinguishable from the GR strain, and we draw our conclusions on the excluded parameter space. We note that GW150914 gives stronger constraints than the second event GW151226. But even with only one observation, we find that for large enough mixing angles we may exclude values of $m_{g}\gtrsim 10^{-22}\,\mathrm{eV}$, comparable to the bounds set by GW150914 via a modified dispersion relation Abbott et al. (2016a). We have adopted the model-independent mass bound from solar system tests, $m_{g}<7.2\cdot 10^{-23}\,\mathrm{eV}$ de Rham et al. (2017), to the present case by multiplying the mass with a factor $\sin\theta$ to account for the bigravity modification of the classical Newtonian potential, see e.g. Platscher and Smirnov (2017). We find that GW oscillations give stronger constraints for smaller mixing angles, where the bound from local gravity tests quickly becomes weaker. In conclusion, GW oscillations offer excellent prospects to probe the bigravity parameter space once more events at higher precision become available.

III Conclusions

We have studied the oscillatory behavior of gravitational waves in the framework of bigravity In full analogy to neutrino oscillations, we have seen that a non-diagonal coupling of the two modes to matter gives rise to potentially significant modulations of the strain that would be observable e.g. in the LIGO or LISA detectors. Using the first ever detected gravitational wave signals GW150914 and GW151226, we illustrated that the bigravity modification of GR can lead to drastic modulations of the strain compared to the predictions of GR. Using this, we have constrained the parameter space of the model in the low-mass regime, and pointed out that, once more events are available, the high-mass regime can be constrained, too.

In this letter, we have made several approximations and assumptions in order to be able to give compact analytic expressions that allow the reader to understand the mechanisms behind gravitational wave oscillations. Nevertheless, the fully general results are obtained easily by following our approach such that future analyses may directly use the results of this work.

acknowledgments

We would like to thank Evgeny Akhmedov for very useful discussions on the fundamentals of neutrino oscillations. We are also grateful to Angnis Schmidt-May and Mikael von Strauss for very useful comments on the manuscript. MP is supported by IMPRS-PTFS.

Note added

Simultaneously to this manuscript, Brax et al. (2017) appeared, where GW oscillations in doubly coupled bigravity are studied. Note that there, the leading effect is proportional to $\tilde{c}-1$ because of the democratic coupling of the tensors to matter.

.