Gravitational wave oscillations in bimetric cosmology

Unlike general relativity, in bimetric gravity linear gravitational waves do not evolve as free fields. In this theory there are two types of tensor perturbations, whose interactions are inherited from non-trivial couplings between two dynamical metric tensor fields in the Hassan-Rosen action, and are responsible for the phenomenon of bigravity oscillations. In this work, we analyze the dynamics of cosmological tensor modes in bimetric gravity on sub-horizon scales and close to the general relativity limit. In this limit, the system has a characteristic length scale $L$ that is strictly contained within the comoving Hubble radius. Thus, depending on the magnitude of the comoving wavelength $\lambda$ relative to $L$, we identify two regimes of interest where the system can be studied analytically: (i) deep sub-horizon modes with $\lambda\ll L$, whose dynamics can be studied using multiple scale analysis and are characterized by small and slowly evolving super-imposed perturbations; (ii) sub-horizon modes with $\lambda\gg L$, where the dynamics is characterized by fast super-imposed oscillations that can be studied using asymptotic techniques for highly oscillatory problems. Furthermore, our analysis represents a substantial improvement compared to previous analyses based on a generalization of the WKB method, which, as we show, is ill-suited to study the system at hand.


Introduction
In the context of general relativity (GR) the gravitational interaction is mediated by the graviton, a massless spin-2 field with non-linear self-interactions [1].Despite the remarkable experimental and observational evidence in favour of GR [2][3][4] and the success of the ΛCDM cosmological model [5], the status of GR as a fundamental theory of gravity has been challenged by modified gravity theories [6,7].Such theories typically lead to modifications in the behaviour of gravity on large scales that could provide an alternative to dark energy and dark matter.Moreover, it has been argued that modified gravity may alleviate or even solve cosmological tensions recently emerged in the ΛCDM model [8].
Bimetric gravity is a highly non-trivial modification of GR that describes consistent nonlinear interactions between massive and massless spin-2 fields [9,10] (see also Ref. [11]).The theory features two dynamical metric fields with suitable interaction terms that ensure the absence of a Boulware-Deser ghost, which plagued previous attempts to formulate a consistent non-linear bimetric theory [12].In the limit where one of the two metrics becomes non-dynamical, one recovers de Rham, Gabadadze and Tolley (dRGT) massive gravity [13].The theory admits yet another independent limit, whereby GR is recovered [14].It has been shown that bimetric gravity is stable and well behaved in certain regions of parameter space [15], and that viable cosmological solutions that fit the expansion history of the accelerating universe also exist [14,[16][17][18].Furthermore, bimetric theory has been confronted with observations of Type Ia Supernovae, Baryon Acoustic Oscillations and the Cosmic Microwave Background [19][20][21], which constrain the parameters of the theory.Further constraints can be obtained from Big Bang nucleosynthesis [22] and from the requirement that there is a working screening mechanism [23].
The focus of this work is on the propagation of gravitational waves in an expanding universe in bimetric gravity, which leads to effects that are potentially testable with LIGO/Virgo gravitational wave observatories and with the forthcoming LISA mission, which will provide data in a yet unexplored frequency band [24].At a perturbative level, the bimetric structure of the theory leads to two distinct types of tensor perturbations (one for each metric) that couple to each other, with time-dependent couplings determined by the cosmological background.In a de Sitter universe, the couplings become constant and the evolution of tensor modes is characterized by an oscillatory behaviour analogous to neutrino oscillations, dubbed gravitational wave oscillations [25,26].Although cosmological perturbations and their propagation have been previously studied by different groups [27][28][29][30][31], general analytical methods are still needed to solve the dynamics for general background evolution and for generic values of the free parameters of the model.A proposal to generalize the standard WKB method to the case of interacting oscillators was made in Ref. [32], and applications to bimetric gravity have been studied in Refs.[32][33][34].However, unlike the standard WKB method, the regimes of applicability of such a generalization remain unclear, due to subtle issues that have been overlooked in previous analyses (these are discussed in detail in Appendix C).
In this work we propose a different method.We focus specifically on the evolution of subhorizon tensor modes on a background dominated by hydrodynamic matter and in the GR limit (as defined in Ref. [14]).Moreover, to keep our analysis as general as possible, we do not assume any specific values for the parameters of the theory in the derivation of approximate analytical solutions.Next, we identify two main regimes of interest for the dynamics of sub-horizon modes, that are separated by a critical scale L = α(aH) −1 determined by the product of the comoving Hubble radius and the parameter α that governs the GR limit.We show that deep sub-horizon modes with comoving wavelength λ ≪ α(aH) −1 can be studied using multiple-scale analysis [35], which is well suited for perturbative problems with dependence on very different evolution timescales and gives far more reliable solutions compared to the generalized WKB method used in Refs.[32][33][34].Moreover, we develop a novel method to study intermediate scales, i.e. for modes with α(aH) −1 ≪ λ ≪ (aH) −1 , which are characterized by fast superimposed oscillations.
The remainder of this paper is organized as follows.In Section 2 we review the formulation of bimetric gravity and its GR limit.The dynamical equations for the cosmological background and tensor perturbations are reviewed in Section 3. Significant simplifications in the structure of the equations are obtained in Section 4, where we exploit the properties of the background evolution on the so-called finite branch.Our original results are contained in Section 5 where we bring the system to a suitable form for the application of multiple-scale analysis.Here we identify the two regimes of interest and characterize them physically.After bringing the system of coupled oscillators to action-angle variables, we solve it using suitable analytical techniques in each regime.Our results are then compared with numerical simulations, showing that the agreement is excellent.Lastly, we review our results in Section 6.
Technical appendices are also included.In particular, Appendix B contains a derivation of the second order action for tensor perturbations in a curved background; although this result has been derived previously in Ref. [18], here we illustrate the calculations in full detail.Lastly, Appendix C includes a brief critical review of the generalized WKB method of Refs.[32][33][34] and a detailed comparison with the multiple scale solution for a toy model with the same structure of the coupled oscillators in bigravity.

Notations and conventions:
We assume the metric signature (− + ++) and units with c = 1.Greek indices are used for spacetime tensors, while Latin indices are used for spatial tensors.

Bimetric gravity: action and field equations
The Hassan-Rosen bimetric gravity action reads [36] S β n e n (S) , ( where R (g) and R (f ) are the Ricci scalars of the metrics g ab and f ab , respectively.Similarly, in the following, we will use the superscripts (g) and (f ) to distinguish elements corresponding to each metric.The constants M g , M f , and m have physical dimensions of mass, whereas the constants β n are dimensionless, and the matrix S is defined in terms of the two metrics g ab and f ab as Finally, the symmetric polynomials of S are defined as [37,38] e 0 (S) = 1 , (2.3a) ) ) ) with Tr[S] = S a a .
In this framework, matter fields are usually assumed to couple only to the metric g ab , which ensures the absence of the Boulware-Deser ghost1 [39,40].Thus, the total action reads where matter fields are collectively denoted as ψ, and are minimally coupled to gravity.The field equations obtained from the variation of Eq. (2.4) read where α ≡ M f /M g is the ratio between the 'Planck masses' of the two metrics, and T µν is the matter stress-energy tensor, defined as usual by The interaction potentials are given by where the matrices The matter action S m is assumed to be diffeomorphism invariant, which ensures the conservation of energy-momentum ∇ µ T µν = 0.As a consequence of the Bianchi identity g µρ ∇ ρ G (g) µν = 0, the divergence of (2.5a) gives the Bianchi constraint where ∇ is the covariant derivative compatible with g µν .The corresponding relations for f µν do not provide additional information, given that the divergences of the interaction contributions satisfy [11,41] Here ∇ is the covariant derivative compatible with f µν .Thus, Eq. (2.9) implies that also the right-hand side of Eq. (2.10) must vanish.

The general-relativity limit
In this section we review the argument given in Ref. [14,42] that shows how general relativity can be recovered from bimetric gravity in the limit where α tends to zero.First, note that the bimetric interaction potentials (2.7) satisfy the following identity: where V is the interaction term in the bimetric action (2.1), that is, β n e n (S) . (2.12) Plugging the field equations (2.5) into the relation (2.11), we obtain which, in the limit α → 0, reduces to Next, taking the covariant divergence of (2.14), as the Einstein tensor and the source T µν are both covariantly conserved, we deduce that ∇ µ V α→0 = 0. Therefore, in this limit V is constant on-shell, and Eq.(2.14) can thus be interpreted as Einstein equations for the metric g µν with a non-zero cosmological constant m 2 V .From here, we conclude that the general-relativity limit of bimetric gravity is defined by α → 0, and that bigravity's self-acceleration survives in this limit-even for β 0 = 0 [14].As can be seen in (2.5b), this is a strong-coupling limit for the f µν sector, which is completely determined in terms of g µν [42].Let us end this section by commenting on another limit of the theory: m → 0. Looking at the form of the Hassan-Rosen action (2.1), it is clear that in the weak-coupling limit m → 0, each metric g µν and f µν evolves independently following its corresponding Einstein equations.However, this limit suffers from the so-called van Dam-Veltman-Zakharov discontinuity [43,44].On the contrary, the limit α → 0 is not affected by such discontinuity.

Linear tensor perturbations of cosmological solutions in bimetric gravity
This section is divided in two subsections.In Sec.3.1 we consider a homogeneous and isotropic solution of the equations of motion of bimetric gravity (2.5), with the matter content given by a perfect fluid.Subsequently, in Sec.3.2, we analyze the dynamics of tensor perturbations evolving on the commented cosmological background.

Dynamics of the cosmological background
Let us assume the closed Friedman-Lemaître-Robertson-Walker (FLRW) geometry as a background for the two metrics f ab and g ab .As shown in Ref. [45], the Bianchi constraints imply that the curvature radii for the spatial geometries must be the same for the two background metrics.Thus, we can write Here a and b are the scale factors, c is the lapse of the conformal time of f ab relative to g ab , and γ ij is the metric on the three-sphere γ ij dx i dx j = r 2 0 dχ 2 + sin 2 χ dθ 2 + sin 2 θ dϕ 2 , with a fixed radius r 0 .The conformal Hubble rates of g ab and f ab will be denoted as H ≡ a ′ /a and H f ≡ b ′ /(bc), respectively, where the prime stands for a derivative with respect to the conformal time η.
The equations of motion for the two metrics, obtained from (2.5), read where ρ m and p m are the energy density and pressure of the perfect fluid, respectively, and satisfy the continuity equation Deviations from general relativity due to bigravity interactions are encoded in the effective energy densities and pressures, which are defined as ) ) ) with y ≡ b/a representing the ratio between the scale factors.We assume that matter obeys the null energy condition as a strict inequality, ρ m + p m > 0. The Bianchi constraint (2.9) is equivalent to the continuity equation for the effective fluid, where Γ(y) ≡ β 1 + 2β 2 y + β 3 y 2 .As it is clear written in the form (3.5), there are two possible solutions to the Bianchi constraint.On the one hand, it is satisfied if the conformal Hubble rates for the two metrics coincide, H f = H: this condition defines the so-called dynamical branch.On the other hand, (3.5) is also satisfied if Γ(y) = 0, which defines the so-called algebraic branch [11].However, it is well known that solutions on this latter branch are unstable [18,46].Hence, we will not take them into account and, in the remaining of the paper, we will focus on the dynamical branch with H f = H.Therefore, using the the condition H f = H, in combination with the two effective Friedmann equations in (3.2), it is possible to obtain an algebraic relation between the different energy densities, Replacing in this relation the definitions of the effective energy densities (3.4) leads to the following polynomial equation for the ratio between the scale factors y An evolution equation for the ratio y can also be obtained from the definitions of the conformal Hubble rates for the two metrics Using the Bianchi constraint on the dynamical branch, H f = H, this equation simplifies to which can equivalently be written as

Dynamics of tensor perturbations
As in general relativity, also in the theory of bimetric gravity, scalar, vector and tensor perturbations evolve independently at the level of the linearized field equations.In this work we focus exclusively on tensor perturbations evolving on the background described in the previous section.Thus, our perturbative ansatz can be written as where h ij and H ij are tensor perturbations, i.e., they are transverse and traceless ) Here D i denotes the Levi-Civita connection associated with γ ij , D i γ jk = 0.In the following, spatial indices will be raised and lowered with the metric γ ij .
The equations of motion can be obtained by expanding the field equations (2.5) to linear order in h ij and H ij , using the perturbative ansatz (3.11), where we have defined the third-order polynomial in the ratio y, and π ij is the anisotropic stress, i.e., the tensorial component of the linear perturbations of T µν .

Harmonic decomposition on S 3
In order to simplify the above equations, and taking into account the symmetries of the background, we decompose the perturbations in tensor spherical harmonics [47], The index s represents the polarity (s = 0 for even/polar harmonics and s = 1 for odd/axial harmonics), and the Q nlms ij are eigenfunctions of the Laplacian on the three-sphere with radius r 0 , By definition, these are transverse and traceless: Note that there are two main differences with tensor modes in a flat universe: (i) the spectrum of the spatial Laplacian is discrete rather than continuous; (ii) the minimum eigenvalue is strictly positive [48].
A similar decomposition is introduced for the anisotropic stress-energy tensor The equations for the harmonic components of the tensor perturbations h ij and H ij are then obtained after replacing the decompositions (3.16) and (3.18) into the evolution equations (3.14), In order to make the notation lighter, from this point on, we will drop the labels n, l, m, s.Thus, the dynamics of tensor perturbations is governed by the following two second-order non-autonomous linear ordinary differential equations,

Mukhanov-Sasaki-like variables
In analogy with the standard Mukhanov-Sasaki variables, which are commonly used in general relativity, it turns out to be very convenient to rescale h and H by their respective scale factors, thus defining the new dynamical variables µ ≡ ah and ν ≡ bH.In this way, the system (3.20) can be recast as where we have defined k 2 ≡ (n 2 − 1)/r 2 0 .In general, the effects due to a spatially closed universe are most relevant for the lowest lying modes in the spectrum since, as n increases, the difference between consecutive eigenvalues tends to zero.However, in the next sections we will focus on the analysis of the dynamics of sub-horizon modes (with large n), for which the effects of spatial curvature are unappreciable.
4 Sub-horizon dynamics of tensor modes in the small y-limit In this section we further specify our assumptions on the background evolution, which will be needed to obtain approximate analytical solutions of Eq. (3.21) describing the propagation of gravitational waves in the next sections.In particular, our aim is to study the dynamics of sub-horizon tensor modes in the GR limit of the theory, specifically in cosmic epochs where the background energy density is dominated by hydrodynamic matter and the ratio y is small.Our motivation to consider the GR limit is that one expects bimetric effects to be small, since GR agrees with observations on all scales.However, even in the regime where departures from GR are strongly suppressed at the background level, the bimetric field equations still give rise to interesting phenomenology for the propagation of tensor perturbations, which is the subject of the following sections.
More precisely, we will consider an expanding universe (with an initial singularity such that the scale factor a vanishes at η = 0), and assume a barotropic equation of state p m = w ρ m , with a constant parameter w > −1, for matter.Then, the continuity equation (3.3) implies that the matter density behaves as ρ m ≈ a −3(w+1) , and thus it diverges at early times a → 0. Concerning the ratio between scale factors y, the quartic equation (3.7) provides four different solutions, which define different branches with distinctive behaviors at early times.Since three of those branches have been shown to be unphysical [49], we will restrict ourselves to the so-called finite branch [50], where y is a monotonically increasing function of η and tends to zero as a → 0 [21].Note that, strictly speaking, for solutions on the finite branch the limit y → 0 is only approached in the very early universe, where the redshift z diverges.Nonetheless, it is possible to satisfy the condition y ≪ 1 until relatively small redshift values, depending on the initial conditions.For instance, the analysis in Ref. [21] shows that there are viable solutions where y ≪ 1 as late as z ∼ 100.
Let us thus solve the background equations under the commented assumptions.For small y, Eq. (3.6) can be approximated as with ρ m > 0 provided that β 1 > 02 .Differentiating, this gives Replacing relation (4.2) in the matter continuity equation (3.3), and making use of Eq. (3.9), we obtain c in terms of the equation-of-state parameter, Using the fact that ρ m ∼ a −3(w+1) , for y ≪ 1 the effective energy density ρ g can be neglected in Eq. (3.2a), whose solution provides the evolution of the background scale factor and the Hubble rate, Substituting this in (3.9), one obtains In these expressions η ∈ [0, ∞) and η * is an arbitrary reference value of the conformal time.Note that all the integration constants have been fixed at the point η * , and that the particular case w = −1/3 is not included in this solution.The results we will present will thus be completely generic for any constant w > −1, excluding w = −1/3.Finally, concerning the tensor perturbations, as already commented, we are going to focus our analysis on sub-horizon modes, defined by the condition k η ≫ 1.For simplicity, in the following we will also assume matter with zero anisotropic stress, π = 0. Using (4.4), it is easy to see that for these modes a ′′ /a ≈ 1/η 2 ≪ k 2 , and, since c is constant given by (4.3), the system (3.21) can be approximated as ) Our last assumption for the physical scenario under consideration is that the bimetric effects are small, and thus we are near the GR limit, so that α ≪ 1.However, this system of equations is still very involved and it is not possible to obtain exact analytical solutions.Therefore, in the next section we will consider well-suited methods to construct approximate analytical solutions.

Approximate analytical solutions: multiple-scale analysis
In order to apply approximate methods, let us first rewrite the system (4.7) in a dimensionless form.In particular, we define a dimensionless time variable t ≡ k η.Furthermore, we introduce a dimensionless parameter ϵ ≡ 1/(kη * ), where η * will be interpreted as the value of the conformal time η at the end of the considered epoch, where the small-y assumption still holds.Since we are assuming sub-horizon modes, kη * ≫ 1 and thus ϵ ≪ 1.However, for the mode to be sub-horizon all along the considered evolution, one would need to define a minimum value of η = η 0 , so that kη 0 ≫ 1.In this way, the following analysis will in principle be valid for all η in the interval η 0 ≲ η ≲ η * .Then, after some straightforward algebra, the system (4.7)takes the form where t ≡ ϵ t is a slow time variable.The time-dependent couplings in Eq. (5.1) are obtained from Eqs. (3.15) and (4.7) after substituting the background evolution.Their expressions as functions of t are where we have introduced the dimensionless parameters B n ≡ m 2 a 2 * β n /H 2 * [51].As discussed in Sec.2.1, the general-relativity limit of bimetric gravity is recovered as α → 0 [14,42].Thus, the system (5.1)naturally involves two small parameters, α and ϵ: the former is fixed at the outset by the values of M f and M g , whereas the latter is by definition mode-dependent (i.e., for a given mode, it depends on the corresponding eigenvalue of the spatial Laplacian).On closer inspection of the system (5.1),we realize that the behavior of solutions is quite sensitive to the hierarchy satisfied by α and ϵ.In fact, the limit of Eq. (5.1b) where both α and ϵ tend to zero clearly depends on the order in which these limits are taken.Thus, we can identify two different types of modes: For modes of Class I the couplings are small and evolve over long-time scales determined by t.Therefore, the system can be studied using techniques from multiple-scale analysis presented in Ref. [35].However, for modes of Class II by no means can the couplings (5.1b) be considered small, and therefore the methods of Ref. [35] are not applicable.Such large couplings are responsible for fast superimposed oscillations that require different methods to be appropriately described.These methods will be developed in Sec.5.2.Finally, modes with ϵ ≈ α, that is, with k ≈ M g /(M f η * ), do not lay in any of the mentioned classes.For such modes, the explicit perturbative parameters in Eq. (5.1b) are absent, which makes difficult to design an appropriate perturbative analysis.
Ultimately, note that the value of ϵ is epoch-dependent and gets smaller as the universe moves to later epochs (which have larger values of η * ).Therefore, Class II modes eventually become Class I in later epochs (as long as the background dynamics is dominated by hydrodynamic matter with w > −1).
Before we proceed further with our analysis, it is convenient to recast the system in a different form.We note that Eq. (5.1) has the same structure as the non-autonomous linear system of coupled oscillators considered in Ref. [35] (see Chapter 4 therein), with ω 2 µ ( t, ϵ) = 1 + ϵ 2 p( t) and ω 2 ν ( t, ϵ) = c 2 + ϵ 2 r( t)/α 2 .This system can then be written in the so-called standard form, transforming to the action and angle variables defined as In terms of these variables, Eqs.(5.1) read, ) This is indeed the form of the system we will consider to obtain approximate solutions.

Class I
We will begin analyzing modes of Class I, which have ϵ ≪ α ≪ 1.The asymptotics of the solutions of system (5.5) in this regime is controlled solely by ϵ.We observe that the system (5.5)belongs to the more general class (5.6b) In our particular case, we have N ∈ {µ, ν} for the labels, while the unperturbed frequencies are constant and given by ω ν = c, and the remaining functions where we have stressed the dependence on t.Moreover, we observe that the condition ϵ ≪ α ≪ 1 implies that (5.5) describes, in the language of Ref. [35], a standard-form system.In order to apply the techniques developed in Ref. [35], based on the method of multiple scales [52], to obtain asymptotic solutions for this kind of systems, the function F N and G N in (5.6) must satisfy the following two conditions: (i) all the functions F N and G N are O(1) in the ϵ → 0 limit; (ii) F N and G N are periodic functions of the φ i with period 2π.Note that, for condition (i) to follow by our functions (5.7), we must require the leading order in r( t), proportional to t−2 , to remain sufficiently small during all over the time interval we are considering.This imposes B 1 ∼ t0 2 α 2 y * (1 + 3w) 2 /4c.The method of multiple scales is able to account for the cumulative effect of small perturbations over long timescales, of order t ∼ ϵ −1 .Thus, the effect of perturbations can be understood as a further dependence of the dynamics on a slow time variable t, in addition to the fast time t of the unperturbed system.The method of multiple scales treats t and t as independent, exploiting the additional freedom to remove secular terms (which would grow monotonically in time) and thus achieve a uniform perturbative expansion over long timescales.
The dependence of the dynamics on the timescales t and t, along with the smallness of ϵ, motivates the following perturbative ansatz for the solution of (5.5), N (τ i , t) + O(ϵ 4 ) . (5.8b) The τ i are the fast times associated with ωN frequencies.As such, they are in principle slowly varying, and thus can be formally defined as where the quantities Ω N are functions of t to be determined later on.Thus, treating the time variables τ i and t as independent, we can represent the total time derivative as (5.10)This is the central equation underlying the method of multiple scales.Note that to zero-th order, φ N depends on τ N -although not on the remaining τ i with i ̸ = N 3 .
3 In fact, it can be shown that a more general ansatz where φ N a priori depends on all τi gives unwanted secular terms.The condition for the elimination of such secular terms then leads to φ N (τN , t) [35].
We substitute the ansatz (5.8) into Eq.(5.6), expand in powers of ϵ and collect the contributions to equal order in ϵ.The next step is to decompose the functions F N and G N into their average and oscillatory parts (the latter are defined as having zero average with respect to all the φ i ).When solving the averaged equations at each order, we must ensure that our solutions are free of secular contributions, following the steps outlined below.Explicitly, the decomposition of F N (J i , φ i , t; ϵ), that is periodic in each φ µ and φ ν , is decomposed into an average term F N (J i , t; ϵ) and an oscillatory term F N (J i , φ i , t; ϵ): where (5.12) and, thus, the oscillatory part can be computed as The decomposition of G N is entirely analogous.Next, each function also needs to be expanded in powers of ϵ.Applying the steps we just outlined to the particular case represented by Eq. (5.7), we obtain i , J i , φ i , J i , φ i , t) + O(ϵ 3 ) . (5.13d) We observe that at order ϵ 2 all the contributions come from the oscillatory terms.This way, we obtain the following differential equations for equal powers of ϵ: O(ϵ) : dJ N , (5.16a) N , (5.16b) N , (5.17a) N . (5.17b) Note that, as ϕ ν is real, we must require r( t) > 0 to ensure consistency of the δ → 0 limit.From the definitions in (5.2) we conclude that, since B 1 > 0 and under our assumptions y * is small and c > 0 (see (4.3)), this condition is satisfied for t up to order 1 as long as We remark that, upon integration of the highly oscillatory functions in Eq. ( 5.28), one gets a O(δ 2 ) contribution since dt F (t) sin δ −2 ϕ(t) ∼ O(δ 2 ) (see Appendix A for details on the assumptions over the functions ϕ and F ).This is a key property of the system that (upon integration) allows us to keep track of the order of the various terms that arise in the expansion for small δ.In particular, we observe that even though the right-hand sides of Eqs.(5.28a) and (5.28c) involves negative powers of δ, these are cancelled out upon integration, so that there are no divergences in the δ → 0 limit.We integrate Eqs.(5.29a) and (5.29b) where u, v are dummy variables (with ũ = ϵ u, ṽ = ϵ v) and t 0 is some initial time where we set initial conditions A 0 = A(t 0 ), B 0 = B(t 0 ).The integral equations (5.31) can be solved iteratively, by repeated substitutions.To include all contributions of order ϵ 2 , δ 2 , we need two iterations4 , which gives (5.32b) Substituting these expansions for A and B into (5.29c), and retaining terms up to second order in ϵ and δ, we obtain (5.33) We proceed similarly for ϕ ν , but in this case we need to include terms up to order δ 4 ; this is necessary for consistency since ϕ ν always appears in the combination ϕ ν /δ 2 in the argument of trigonometric functions.Thus, we obtain We can write the solutions to Eqs. (5.33) and (5.34) as having defined The integrals in (5.35) include higher-order terms that can be removed by iteratively substituting the right-hand sides of (5.35a) and (5.35b) in the argument of the trigonometric functions.At this point, all the arguments in the trigonometric functions are φ µ or ϕ ν /δ 2 -both of them being non-oscillatory functions, which allows us to make use of the asymptotic formulae in Appendix A. Thus, evaluating the remaining integrals we obtain Having finally obtained approximate solutions for φ µ and ϕ ν , A and B can be obtained substituting the expressions (5.37) in Eq. (5.32), and expanding to order ϵ 2 and δ 2 .After integration, the solutions read (5.38b) Figures 3 and 4 show the evolution of the approximate solutions obtained in (5.37) and (5.38), along with the corresponding numerical solutions of (5.5).As in the previous section, the absolute error of the perturbative solutions is consistent with the order of the asymptotic expansion considered.Nevertheless, we note that the approximation here considered breaks down at the nodes of the envelopes of the solution, as it is clear from Fig. 5.This is because the leading-order asymptotics of the highly oscillatory integrals in (5.38) vanish at the nodes; therefore, in order to accurately describe this region one should include higher-order terms.More specifically, in our case this means including O(δ 4 ) terms in the asymptotic expansions (see Appendix A).Thus, since the distance between two nodes is ∆t = π, our current solutions describe correctly the superimposed oscillations, although only during a half-period of the µ modes.

Conclusion
In this paper we have analyzed the propagation of cosmological tensor perturbations in bimetric gravity in the GR limit M f /M g ≪ 1, focusing on sub-horizon modes kη ≪ 1 in a universe dominated by hydrodynamic matter.Mathematically, the system can be described as a pair of linear oscillators with time-dependent frequencies and couplings.The dynamics is governed by two dimensionless parameters, α ≡ M f /M g and ϵ ≡ 1/(kη * ) (where * denotes the end of the cosmic epoch at hand).We have identified two main regimes of interest (corresponding to two distinct classes of modes) determined by the relative magnitude of ϵ and α .There are significant physical differences between such regimes and each requires specific analytical techniques for its study, particularly for the derivation of approximate analytical solutions.For Class I modes, with k ≫ (M g /M f )H * , the couplings of the dynamical system are slowly evolving, and solutions can thus be obtained with multiple-scale analysis, using the methods of Ref. [35].In this regime, α can be treated as a finite quantity whereas ϵ is infinitesimal.After rewriting the system in action-angle variables, a perturbative scheme is adopted where each term of the perturbative series in ϵ depends on two time variables t and t ≡ ϵt (respectively, fast and slow times)-treated as independent.The method enables us to find uniform approximations to the solutions over a time interval of the order 1/ϵ without introducing secular terms.Even though the couplings are slowly evolving, the action variables for the two oscillators are not conserved over time, thus signalling a breakdown of adiabaticity, as shown in Figs 1 and 2. In terms of the original oscillator variables, we observe that to zeroth and first order in ϵ the main effect consists of a slow modulation of the amplitude and frequency of the oscillators, while the couplings between oscillators only give rise to higher-order corrections.
On the other hand, for Class II modes (H * ≪ k ≪ (M g /M f )H * ) the system is strongly coupled and different analytical techniques are thus required for its study.In this regime both parameters α and ϵ must be treated as infinitesimal and the evolution of the system is characterized by fast super-imposed oscillations, whose frequency is divergent in the limit where δ 2 ≡ α/ϵ tends to zero.In order to capture this crucial feature of the system, we have developed a suitable perturbative   scheme based on the asymptotics of highly oscillatory integrals, as a double expansion in δ and ϵ.Also in this case the action variables are not conserved, and there are significant deviations from adiabaticity due to rapid oscillations on short time scales, as illustrated in Figs. 3 and 4. Also in this case, our approximation scheme is accurate over time scales of the order t ∼ 1/ϵ.Our approximation could be further improved by systematically including higher-order corrections, which are necessary to properly account for the behavior of the action variables near the nodes of the envelopes, as shown in Fig. 5. Lastly, we remark that the dynamics of tensor modes in an expanding universe in bimetric gravity has been studied earlier in Ref. [32] (see also [33,34]), using a scheme that attempts to generalize the standard WKB method to the case of interacting oscillators with time-dependent couplings and friction.The salient features of this approach are reviewed in Appendix C.Here we show that the method has not been developed in a fully consistent fashion, and provide a concrete example where it fails to provide a reliable approximation to the numerical solution, and it is in fact much less accurate than the multiple-scale method.Thus, future work should be devoted to a careful re-evaluation of the projected constraints on gravitational wave oscillations derived in [33] using the methods here developed.

B Second-order action for tensor perturbations
We expand the Hassan-Rosen action (2.1) to second order in the metric perturbations h ij , H ij , around a solution of the background field equations HR + . . ., (B.1) where S (r) HR denotes the r-th order variation.Moreover, by definition, the first variation S HR vanishes if the background obeys the field equations.In this appendix, we will assume the vacuum field equations for reasons of simplicity; however, the generalization to the case where matter is included is straightforward.
In order to evaluate the right-hand side of Eq. (B.1), we compute the second order expansion of the matrix S, which appears in the definition of the bigravity interactions, Then, recalling the definitions of the symmetric polynomials (2.3), we obtain the following expansions For the volume elements we have The first two terms in Eq. (2.1) have the same structure as the Einstein-Hilbert action; therefore, their expansion proceeds as in general relativity.We have, up to a total divergence Recalling that the conformal time for f ab satisfies dη = c(η)dη, from Eq. (B.7) we obtain, after some straightforward substitutions Then, combining Eqs.(B.5), (B.7) and (B.8), and using the background field equations (3.2) in vacuo (i.e., with ρ m = p m = 0) we compute the zeroth-order action evaluated on a solution (Hamilton's principal function) Similarly, from Eqs. (B.6), (B.7) and (B.8), and using the (vacuum) background field equations (3.2) to simplify the coefficients of the terms proportional to h ij h ij and H ij H ij , we can finally compute the second-order action for tensor perturbations on a curved FLRW background EH [g ab ] + S (2) (B.10) The above derivation of the second-order action (B.10) also holds with minimal modifications in the presence of matter fields that contribute to the background field equations.However, in this case there are additional pressure contributions further to the effective pressure terms in (B.6), leading to analogous cancellations with the background quantities in Eqs.(B.7) and (B.8).Furthermore, even though we focused specifically on a closed FLRW model, the derivation of the second-order action is entirely analogous in open and flat cosmologies.Specifically, the corresponding action for an open universe can be formally obtained from (B.10) with the replacement r 2 0 → −r 2 0 , while the case of a flat universe is recovered by letting r 2 0 → +∞.
Lastly, we observe that in the presence of matter perturbations with a non-vanishing anisotropic stress (i.e., a non-zero transverse traceless component of the perturbed matter stress-energy tensor), here denoted as π ij , the second-order action includes the extra term Here we have assumed that all matter fields only couple to the metric g ab , in order to avoid the Boulware-Deser ghost.Therefore, we do not need to introduce an analogous source term for the perturbations of f ab .

C Shortcomings of the generalized WKB method
Here we review previous attempts at a generalization of the WKB method, proposed in Refs.[32][33][34] to solve similar differential equations as the ones considered in this paper.A general model is assumed for the cosmological dynamics of coupled tensor modes, which reads (in matrix form) Here Î is the identity matrix, k is the spatial momentum, and ν, Ĉ, Π, M are the friction, velocity, chirality and mass matrices, respectively.The bigravity case, Eq. (3.20), can be obtained as a special case of Eq. (C.1).In analogy with the standard WKB method [54], a small dimensionless parameter δ (in this case unrelated to the physical parameters of the model) is then introduced by hand to suppress time derivatives In this way, the determination of Φ is then framed as finding an asymptotic solution to the singular perturbation problem (C.3) in the δ → 0 limit.The following ansatz is then made, as a generalization of the usual WKB ansatz that also accounts for mixing Φ = Ê e where Ê is the matrix of eigenvectors and θ is the diagonal matrix of eigenfrequencies. 5he ansatz (C.4) is then substituted in the equations of motion (C.3), which gives with a( t) = p t 3 + q t 4 , b( t) = t 2 a( t) , w( t) = r t−1 .The plot clearly shows that the generalized WKB method fails to give an accurate approximation of the numerical solution, which is instead well approximated by the multiple-scale solution.

Figure 5 :
Figure 5: Nodes of the envelope corresponfing to the solution of I ν in Fig. 4.
B.4) where √ γ is the spatial volume element associated with γ ij .Finally, combining Eqs.(B.3) and (B.4), and recalling the definitions of the effective pressures (3.4), we obtain the following expansion for the interaction terms in the bimetric action (2.1), S