Gauge-invariant gravitational waves in matter beyond linearized gravity

Modeling the propagation of gravitational waves (GWs) in media other than vacuum is complicated by the gauge freedom of linearized gravity in that, once nonlinearities are taken into consideration, gauge artifacts can cause spurious acceleration of the matter. To eliminate these artifacts, we propose how to keep the theory of dispersive GWs gauge-invariant beyond the linear approximation and, in particular, obtain an unambiguous gauge-invariant expression for the energy–momentum of a GW in a dispersive medium. Using analytic tools from plasma physics, we propose an exactly gauge-invariant ‘quasilinear’ theory, in which GWs are governed by linear equations and also affect the background metric on scales large compared to their wavelength. As a corollary, the gauge-invariant geometrical optics of linear dispersive GWs in a general background is formulated. As an example, we show how the well-known properties of vacuum GWs are naturally and concisely yielded by our theory in a manifestly gauge-invariant form. We also show how the gauge invariance can be maintained within a given accuracy to an arbitrary order in the GW amplitude. These results are intended to form a physically meaningful framework for studying dispersive GWs in matter.


I. INTRODUCTION
Analytical studies of gravitational waves (GWs) usually focus on the nondispersive vacuum modes propagating at the speed of light [1], which are also the relevant modes for the observations so far [2][3][4].However, in the early Universe and possibly also near compact objects, the coupling of GWs with matter can be nonnegligible.Understanding this coupling is potentially important, for example, for understanding the electromagnetic signatures of the GW radiation [5,6].Also, since this coupling is small for the tensor modes and waves in cold matter [7], even small corrections due to thermal effects [8], viscosity [9], and matter-induced modification of the wave polarization [10,11] may be significant.
The tools for modeling dispersive GWs in matter can be imported from electrodynamics, plasma physics in particular, where the wave-induced oscillations of matter are commonly described in terms of electric susceptibility [12].Provided that matter is continuously distributed in spacetime, its effect on a GW can be similarly, and conveniently, described in terms of gravitational susceptibility [11,13].Then, by analogy with electrodynamics of continuous media [14], two questions arise: (i) How can one derive the equations of GW propagation at nonzero gravitational susceptibility of the ambient medium?(ii) How can one describe the exchange of the energy-momentum between dispersive GWs and matter?
For small wave amplitudes, the first question can be answered within linear theory [11], but the second question cannot.The energy-momentum of a linear wave produces nonlinear average forces on the underlying medium, and is itself quadratic in the wave amplitude.Defining it generally requires so-called quasilinear (QL) theory [15][16][17], specifically, QL theory of adiabatic (nonresonant) wave-medium interactions.Such interactions are commonly described using variational methods [18] that originate from Whitham's average-Lagrangian approach [19,20].Specifically, an adiabatic response of a medium to a wave is expressed through the wave field and the medium susceptibility; then the wave Lagrangian is identified as the part of the wave-matter Lagrangian that is quadratic in the wave amplitude, whence the wave energy-momentum and the nonlinear forces on the medium are readily inferred in the usual manner [21].This approach has been particularly fruitful in plasma theory in application to the electromagnetic interactions [22][23][24][25][26][27], where the wave energy-momentum and the nonlinear forces on the medium are naturally expressed through the electric field [21].Whitham's approach has also been used to describe GWs [28][29][30][31][32] and allows to circumvent the various problems [33][34][35][36][37][38] associated with averaging tensor fields in a self-consistent covariant fashion [39].Still, application of this approach to GWs remains deficient in that, unlike for electromagnetic waves, GW's average Lagrangians generally lack gauge invariance [28].
As a reminder, the gauge invariance for GWs derives from general covariance [40,Sec. 7.1].Consider a background metric g αβ = O(1) and a perturbation metric h αβ = O(a) on top of it, where a ≪ 1 is a small parameter.Then, a coordinate transformation x µ → x ′µ = x µ + ξ µ , with ξ µ = O(a), implies g αβ → g ′ αβ = g αβ and h αβ → h ′ αβ = h αβ −£ ξ g αβ , where £ ξ is the Lie derivative along the vector field within linearized theory (Sec.III B).If one considers h αβ as a tensor field on the background spacetime g αβ , this transformation can be viewed as a gauge transformation (with ξ µ being the gauge field) and, by general covariance, cannot have measurable effects.However, the GW Lagrangian is normally formulated in terms of h αβ [31] and not gauge-invariant in the presence of matter [28].That leads to a non-invariant energy-momentum of the wave and the possibility of spurious acceleration (or heating) of the matter with gauge artifacts, which is unacceptable.This brings the question: how can one construct a gauge-invariant QL theory of GWs and, in particular, define the energy-momentum of a dispersive GW such that it would be exactly free of gauge artifacts?
Here, we show how to do this using the tools that we have developed previously in Refs.[11,13,17,41,42].Let us briefly describe this series of papers to put our work into perspective.Our general formalism for studying dispersive GWs interacting with matter is outlined in Ref. [11] and, within a broader context, in Ref. [17].The work [13] also specifies some elementary blocks of our general theory, such as the so-called ponderomotive potential and the gravitational susceptibility.We discussed applications of our theory to GWs in neutral gases in Ref. [11] and non-magnetized plasmas in Ref. [42], focusing mostly on linear dispersion relations.A QL theory of GWs has been proposed in Ref. [17], but it is gaugedependent by construction and thus free of coordinate artifacts only asymptotically.We have also discussed the manifestly gauge-invariant approach to GWs in a general metric, but only in the context of a linear theory so far [41].Here, we merge those two approaches and, for the first time, formulate a manifestly gauge-invariant theory of dispersive GWs beyond the linear approximation.
Specifically, we show how to rewrite the linear-GW action, approximately yet without loss of accuracy, in terms of the gauge-invariant part (projection) of the metric perturbation that was derived in Ref. [41].We accomplish this for a general background metric, so our approach is applicable to waves propagating in matter arbitrarily (but smoothly) distributed in spacetime.For this, we assume the short-wavelength limit, which allows us to bypass the problems [33][34][35][36][37][38][39][43][44][45][46][47][48] associated with covariant self-consistent averaging on curved manifolds.1 Then, geometrical-optics (GO) equations are derived for linear GWs from the approximated action as usual [21,49,50].This approach is commonly known to provide a particularly simple and robust way of deriving reduced equations [22,26,27], especially compared to averaging differential equations directly [51,52].Next, we identify the action and the energy-momentum of a dispersive GW in a gauge-invariant form and derive the nonlinear effect of GW on the background medium. W also show how gauge invariance can be maintained within a given accuracy if nonlinearities are included up to an arbitrary order in the GW amplitude.An application to vacuum GWs is presented as a simple example that connects our general formulation with commonly known results.More importantly, though, our results show, for the first time, that a reduced QL theory of GWs in matter can be entirely freed of gauge artifacts and, as such, represents a physically meaningful framework to study specific GW modes in matter in the future.
Our paper is organized as follows.In Sec.II, we introduce our basic equations.In Sec.III, we propose a reduced variational formulation for weak short-wavelength GWs on a smooth background spacetime and put it in a gauge-invariant form.In Sec.IV, we introduce the equations for the perturbation metric and our gauge-invariant formulation of GO of GWs.We also discuss how vacuum GWs are described within this framework, and in particular, how their action and energy-momentum are naturally introduced in a gauge-invariant form.In Sec.V, we present the model for how weak dispersive GWs influence the background metric.We also discuss vacuum waves as an example.In Sec.VI, we present our second approach to gauge-invariant GW theory that extends beyond the QL approximation.In Sec.VII, we summarize our main results.

II. PRELIMINARIES A. Einstein equations
Let us consider a metric g αβ with signature (− + ++) on a four-dimensional spacetime with coordinates denoted as x α .The dynamics of this metric is governed by the least-action principle [53] δS = 0, S = S m + S EH . ( Here, S m is the action of the matter, S EH is the vacuum action of the gravitational field, called the Einstein-Hilbert action, R is the Ricci scalar, g .= det g αβ , and .= denotes definitions.By default, we assume units such that the Einstein constant κ ≡ 8πG/c 4 and the speed of light c are equal to unity, The equations for g αβ , called the Einstein equations, are obtained from where g αβ is the inverse metric (g αβ g βγ = δ α γ ) and [g] denotes that the result is evaluated on g αβ .Using where G αβ is the Einstein tensor and T αβ is the energymomentum tensor of the matter, Eq. ( 4) can be represented as We will assume, for clarity, that the matter is not ultrarelativistic; then each element of T αβ is O(ρ), where ρ is the mass density.

B. GW and average metric
Here, we consider a class of problems within which the GW propagation is typically studied.Specifically, let us assume the existence of, and restrict our consideration to, the coordinates such that g αβ can be decomposed (in a way to be specified shortly) as where each element of g αβ is O(a 0 ) and each element of h αβ has a magnitude that does not exceed a small constant a ≪ 1.Furthermore, we assume that h αβ has two separate scales: ℓ h , which is the characteristic wavelength of GWs of interest, and ℓ g ≫ ℓ h , which is the scale of the background or its radius of curvature. 2Then, our definition of the scale separation is the same as that of the high-frequency limit in Ref. [28], where a detailed discussion about its physical relevance can also be found.
Because of the assumed scale separation, one can find a scale ℓ a that satisfies and one can introduce a small "GO parameter" More formally, we assume that any h αβ of interest is a superposition of quasiperiodic functions, i.e., functions of (ǫx, θ(x)), where ǫ ≡ (ℓ h /ℓ a ) 2 is a small dimensionless parameter and the dependence on θ is 2π-periodic.Then, any function A (generally, a tensor) that contains g αβ is also quasiperiodic with scales that satisfy (8) and can be assigned a local average . . .over a spacetime volume of size ℓ a as Here, Ψ is a fixed function that is localized on the scale ℓ a near x ′ = x and is normalized such that d 4 x ′ Ψ(x, x ′ ) = 1.The shape of Ψ can be arbitrary, because of the following.As a quasiperiodic function, A can be represented as a Fourier series in θ: Then, A n (ǫx)e inθ(x) is exponentially small with respect to √ ǫ except for n = 0, and A 0 (ǫx) ≈ A 0 (ǫx) with error of order √ ǫ, because on the scale ℓ g , the function Ψ is close to delta function.This means that A ≈ A 0 , which becomes an exact equality at ǫ → 0 independently of the shape of Ψ and ℓ a as long as Eq. ( 8) is satisfied.(Corrections caused by nonzero ǫ will not be important for our purposes.)This subject is discussed in detail in Ref. [17].Furthermore, notice the following.The aforementioned restriction on the choice of background coordinates implies that any coordinate transformation that is O(a 0 ) is assumed to have the scale ℓ g or larger.[See also Eqs. (44) and the ensuing discussion.]This means that, at ǫ → 0, such transformations preserve the expansion (10) in that each A n remains the nth Fourier harmonic of A with respect to the phase θ, which is a true scalar.Then, each coefficient A n transforms as a tensor independently from the other coefficients in the expansion (10), and thus, in particular, our averaging procedure is frame-invariant at ǫ → 0. (It is also possible to use other averaging schemes [39,47,55], and those are known to produce equivalent results [43][44][45][46]48] under the condition (8a).) We now specify the splitting ( 7) by requiring that h αβ satisfies the following condition that is invariant with respect to allowed background transformations at ǫ → 0: We call such a perturbation a GW.Then, g αβ can be understood as the background metric for the GW or as the average part of the total metric: For any pair of (complex) fields u 1 and u 2 on the background space, we introduce the following inner product: where g .= det g αβ .We also introduce the inverse background metric g αβ via g αβ g βγ = δ α γ .Then the inverse total metric g αβ can be expressed as follows: (Here and further, the indices of the perturbation metric h αβ are manipulated using the background metric and its inverse, unless specified otherwise.)We will assume that the matter responds adiabatically to GWs, so g αβ and h αβ are the only degrees of freedom, meaning that resonant interactions of waves with matter are ignored (but see Ref. [17] for these effects).This is a standard approach that was used in the past for electromagnetic and other interactions, for example, in Refs.[13,18,24,25,56,57] and is reviewed, for example, in Ref. [26].
The sign conventions will be assumed as in Refs.[40,58].Then, in particular, the commutator of two covariant derivatives acting on any vector field u α can be expressed as where R αβ ≡ g αγ R γ β is the the Ricci tensor of the background metric.In addition to these, we also introduce the operator and we define Ξ α β as the Green's operator of the following equation for u α : where q α is a given vector field.In other words, Ξ α β is such that Notably, Q α β is a hyperbolic operator similar to the one from the driven Maxwell's equation for the Lorenzgauge in vacuum [59] except for the opposite sign in front of the Ricci tensor.We will assume the adiabatic limit, when Q α β is approximately invertible due to the scale separation (8), in which case one can write Further details can be found in Ref. [41].

III. VARIATIONAL APPROACH
A general approach to QL theory involves the Weyl calculus, as detailed in Ref. [17].However, it requires cumbersome machinery that is redundant for describing adiabatic waves that we are interested in here.For the purposes of this paper, it is sufficient to use a simpler and more intuitive Whitham's method, as motivated in Sec.I.This is the approach we adopt below.That said, readers are also encouraged to consider Ref. [17] as an affirmation of Whitham's method for adiabatic GWs (in particular, see Sec. 9.4 there), even though Ref. [17] is not concerned with exact gauge invariance that we pursue here.

A. Basic equations
Using Eq. ( 14) for the inverse total metric, the total action S can be expanded in h αβ as follows: Here, the square brackets [g] emphasize that the corresponding expansion coefficients, such as the matrix function σ αβ and the operator D αβγδ , are evaluated on g αβ .Also, without loss of generality, one can assume that σ αβ = σ βα and where the dagger denotes Hermitian adjoint with respect to the inner product (13).
Because h αβ = 0, the linear term in Eq. ( 20) vanishes.[To reiterate, this is possible only within the shortwavelength approximation (8).]Let us also assume for clarity that the three-wave interactions are negligible (a sufficient condition for this is that the GW spectrum must be not too broad); then the average of terms cubic in h αβ vanishes too [60,Sec. 1.1.2].This leads to where ∆S = O(a 4 ).The term S (2) can be represented as a sum of the vacuum action S vac and the term S m that describes the GW-matter coupling: Assuming a 2 = o(ρ), the term ∆S in Eq. ( 22) is small not only with respect to S vac , but also with respect to S (2) m .(The formulation presented below is valid also at arbitrarily small ρ to the extent that the GW interactions with matter can be neglected.)Then, the system can be described with the following truncated action that retains not only the vacuum linearized gravity but the GW-matter interactions too: An example of S (2) [g, h] that emerges from the quasimonochromatic-GW interaction with neutral gas was derived in Ref. [13] as a part of studying the ponderomotive effect of GWs on matter.

B. Coordinate transformations
Let us explore how the action ( 22) is transformed under near-identity coordinate transformations where ξ µ is a small vector field.A transformation (25) induces the following transformation of the total metric: where £ ξ is the Lie derivative along ξ µ [40], Let us specifically focus on the case ξ µ = 0.Then, using Eq. ( 14) together with Eqs. ( 11) and ( 12), one obtains the following transformation of the perturbation metric, or a gauge transformation with ξ µ serving as a gauge field: Also, the background metric is transformed accordingly as From Eq. ( 24), S[g ′ ] can be written as ), (30) where the difference between D αβγδ [g ′ ] and D αβγδ [g] has been absorbed into O(a 4 ).Because S[g ′ ] = S[g] by general covariance, this leads to where we used Eqs.(21a) and (28).Notice that where we used Eq. ( 29) and the fact that σ αβ [g] = O(a 2 ) by the Einstein equations.After substituting this into Eq.( 31), one obtains Because ξ µ is arbitrary, one can also flip the sign of ξ µ to obtain Together, Eqs. ( 33) and (34) lead to As a side remark, note that within the linear approximations, the O(a 4 ) terms in Eqs.(30) [and the O(a 2 ) term in Eq. ( 29)] are negligible.Then the right-hand side of Eqs. ( 35) is zero, which shows that linear theory is covariant under transformations h µν → h ′µν = h µν + £ ξ g µν , or gauge-invariant.
In contrast, the reduced field theory induced by the approximate action (24) is not strictly gauge-invariant; i.e., S generally does change under coordinate transformations (25).These changes are O(a 4 ), so they are beyond the accuracy of our approximation and therefore do not invalidate the reduced theory per se.Still, the lack of covariance makes the reduced theory not entirely satisfactory, as discussed in Sec.I. Below, we show how to derive an alternative approximation of S that, on one hand, is equivalent to Eq. ( 24) within the accuracy of our theory but, on the other hand, is exactly gauge-invariant.

C. Projecting on the invariant subspace
As shown in our recent Ref. [41], the perturbation metric can be uniquely decomposed into the gauge-invariant part ψ αβ and the remaining gauge part that is the Lie derivative of some vector field ζ α : (The notation here is slightly different than in Ref. [41].)Specifically, where the linear operator Ξ α β is defined in Eq. (19).Also, the linear operator Π αβ γδ is given by and has the following properties: where u µ is any vector field.Substituting Eqs.(36) into Eq.(20c) for S (2) yields where we used Eqs.(21).By Eqs.(35), the second and the third term here are O(a 4 ) and therefore negligible within the accuracy of our theory.Hence, only the first term should be retained.Using Eq. (36b), one thereby arrives at the following approximation: where we have repeated Eq. ( 24) for completeness.Equation (41) defines a field theory where g αβ and h αβ are independent fields.The advantage of this approximate field theory is that it is exactly invariant under gauge transformations for any ξ µ .Specifically, by Eq. (39b), one has Let us also consider a class of transformations where x = x(x ′ ) is a prescribed four-dimensional function.This can be understood as the zeroth-order coordinate transformation on the original spacetime (or a regular transformation of the background coordinates).Since Π αβ γδ [g] is defined covariantly with respect to such transformations, the object Π γδ ελ h ελ is a true tensor on it.Similarly, covariance of D αβγδ [g] and S[g] ensures that S (2) and the whole action (41a) are true scalars on the background spacetime.Thus, in addition to the gauge symmetry, the field theory ( 41) is also invariant with respect to transformations (44).

D. Metric perturbation as a vector field
Following Ref. [41], let us consider h αβ as a 16dimensional field h a , or h in the index-free notation, of the form h = (h 00 , h 01 , h 02 , h 03 , h 10 , . . ., h 32 , h 33 ) ⊺ , where ⊺ denotes transpose.(Here and further, Latin indices from the beginning of the alphabet range from 1 to 16.)In other words, where the index function ι is defined via assuming a, b = 1, 2, . . ., 16. Accordingly, assuming the same notation as in Eq. ( 47).
Let us define H 1 as a Hilbert space of one-component functions on the background spacetime with the usual inner product • , • .Then, the 16-dimensional fields (45) can be considered as vectors in the Hilbert space H 16 that is the tensor product of 16 copies of H 1 , with the inner product where • , • is given by Eq. ( 13).(Unlike in the rest of the paper, summation is shown explicitly here in order to emphasize the difference between • | • and • , • .)Then, Π αβ γδ induces an operator Π a b on H 16 defined via By Eq. (39a), one has where φ αβ .= £ ξ g αβ and ξ µ is any vector field.This shows that Π is a projector of the metric perturbation on the gauge-invariant subspace.However, note that Π † = Π, so Π is not an orthogonal but oblique projector.[Here, the dagger denotes Hermitian adjoint with respect to the inner product (50).]Using this machinery, one can rewrite Eqs.(41) in the following compact form: Here, we have introduced the operator which is Hermitian by Eqs.(21).The indices here can be raised and lowered using γ ab [Eq.( 49)] as a metric.

E. Metric perturbation as a complex variable
For dynamics governed by Eq. ( 24) to be consistent with the assumption of slow g αβ , the integrand in S (2)  must be properly averaged.This is done as follows.Because h αβ is assumed to be rapidly oscillating, one can unambiguously split it into the part hαβ that corresponds to positive frequencies and the complex-conjugate part hαβ * that corresponds to negative frequencies: which is always doable provided the scale separation (8) [61].Also, one can write hαβ as a sum (possibly, integral) over all quasimonochromatic waves present in the system: where the envelopes a αβ s and the local wavevectors kµ .= ∂ µ θ s are slow compared with θ s and the sum is taken over all the waves that are present in the system.In Eq. (53d), contributions of hαβ hγδ and h * αβ h * γδ average to zero, and so are terms like h * αβ s ′ hγδ s with s ′ = s.Then, Using Hermiticity of D, one can drop "Re " and rewrite S QL .

F. Continuum wave spectrum
As usual [16,17], in the case of a continuum wave spectrum, S (2) can be expressed through the spectrum of the correlation matrix of h: where s is a 4-tuple (not a vector), so x ± s/2 is to be understood as the index-free notation for x α ± s α /2, and h a (x±s/2) is the metric perturbation as a 16-dimensional vector field h a evaluated at the coordinates given by x ± s/2.Note that W ab is not a truly tensorial object in the general case.[See the discussion following Eq.( 61).]It is also known as the (average) Wigner matrix and can be understood as the (average) Weyl symbol of the "density operator" |h h|, at least up to a constant coefficient [17,49].Here, J is a metric factor given by [62,63] J(x, s) .
As a reminder, g(x) is the determinant of the background metric evaluated at x.Much like in Refs.[16,17], one obtains where D ab is the Weyl symbol of D ab .In the index-free form, this can be written as As discussed in Ref. [50], the matrices D ab and W ab are not necessarily true tensors.However, notice the following.In the GO limit, for any quasimonochromatic wave with a given local wavevector k α , one has Dh → D(x, k)h.Since D is a covariant object, D is thereby a true tensor in this limit.Likewise, W ab becomes a tensor density in the GO limit.This is, in particular, because typical s α that contribute to the integral in Eq. ( 58) are of the order of the correlation length (which in the GO limit is much less than ℓ g ), so J(x, s) can be replaced −g(x).Then, s can be approximately interpreted as a vector and x ± s/2 can be understood as a small displacement of x along this vector.More details about Wigner matrices and their role in QL theory can be found in Ref. [17] and references cited therein.We are not focusing on this topic here because it is not needed for the rest of our paper, but see also Sec.V C.

IV. EQUATIONS FOR THE PERTURBATION METRIC A. Basic equations
First, let us consider linear theory, where g αβ is treated as a prescribed function.In this case, the field equations yielded by the action (57) are as follows: These two equations are the complex conjugates of each other, so only one of them will be considered.Specifically, one obtains so D is understood as the dispersion operator for GWs.Also note that because this equation is linear, one can equally write it as and also as By the definition of D [Eq.(53e)] and by the properties of Π [Eq.( 52)], this equation is invariant with respect to the gauge transformations (42).

B. Geometrical optics of linear quasimonochromatic waves
Let us consider a single quasimonochromatic wave where the envelope a αβ and the local wavevector kµ are slow compared with θ by a factor of ǫ ≪ 1.Using the Weyl expansion of the dispersion operator [50,64], Eq. (63b) can be reduced to the following first-order equation: where The field equation ( 65) is gauge-invariant only in the limit ǫ → 0, but a strictly gauge-invariant theory can also be constructed at nonzero ǫ, namely, as follows.
As seen from Eq. ( 65), D h ≈ D(x, k) h, so Eq.(57b) can be reduced to At small ǫ, one also has where Π is the Weyl symbol of Π as in the Minkowski space [41]: [Here, Greek indices are converted to Latin indices as in Eq. ( 47), k 2 .= k µ k µ , and the brackets in indices denote symmetrization, as usual.]Using this, let us approximate Eq. (67) as where L is the Lagrangian density given by and a † is a row vector with elements a b * .By the properties of Π, this L is exactly invariant with respect to the gauge transformation (where ξ α is an arbitrary vector field), which represents the GO limit of transformations (42).Hence, so is the whole field theory generated by Eq. (70).Because S (2) has the standard form of the GO action that governs a linear dispersive wave, this theory is readily constructed as usual [21,49]. 5First, let us express the vector h through the scalar amplitude a .= |a| and the local polarization e .= a/|a|.The polarization and the dispersion relation are found from Eq. (63b), which in the GO limit yields Then, L can be expressed as Hence the dispersion relation can be expressed as The corresponding flux density of the wave action is or equivalently, J α = −a 2 e † V α 0 e/4, where The corresponding action-conservation theorem is [21] the energy-momentum tensor of the wave is [21] and satisfies Note that T α β is the canonical energy-momentum of a wave and, as such, does not have to be symmetric.On the definition of the wave energy-momentum, see, e.g., Refs.[21,24] and, particularly, Ref. [17] for the relevance of the canonical energy-momentum in so-called oscillation-center QL theory.Together with the GO ray equations [49] where τ is a parameter along the ray, the above equations provide a complete gauge-invariant GO theory of adiabatic GWs in a general medium.
C. Example: quasimonochromatic vacuum GWs in the normal coordinates In the limit of vanishingly small density (and smooth enough background metric), one has . Then, using Eq. ( 52), one obtains Also, the GW dispersion relation in this case is k 2 → 0, so Π becomes singular.However, T α β must remain finite, which means that e † V α 0 e must remain finite.Let us assume the normal coordinates, so that the metric is locally close to the Minkowski metric.Then, by performing a direct calculation, one finds that polarization in this case must satisfy e 00 − 2e 03 + e 33 = 0, (84a) e 01 − e 13 = 0, (84b) e 02 − e 23 = 0, ( at kα → (ω, 0, 0, ω).One can recognize these as the conditions for certain metric invariants [41] to be zero.One can also express L in terms of metric invariants directly.Any combination of those are also invariants, and we choose to work with the following set (cf. Ref. [41]): and also where we assumed kα = (ω, 0, 0, k) with nonzero ω.Then, a direct calculation shows that Eq. ( 71) can be rewritten as follows: From δS (2)  δΨ a = 0, δS (2)  δΨ * a = 0, (88) one finds that Ψ 3 = Ψ 4 = Ψ 5 = Ψ 6 = 0; in particular, Ψ 5 = Ψ 6 = o(k 2 ) at k 2 → 0. Assuming the notation Ψ .= (Ψ + , Ψ × ) ⊺ Eq. (87) can be simplified as which describes the two well-known vacuum GW modes with dispersion relation k 2 = 0. Also, by Eqs. ( 76) and (79), the density of the wave action flux and the canonical energy-momentum tensor are given by D. Quasimonochromatic vacuum GWs in general background coordinates The above results for the normal coordinates can be readily generalized to an arbitrary background coordinates.As usual [20,26], in order to obtain leading-order GO equations, one can approximate the wave Lagrangian density to the zeroth order in ǫ. 6 Then, the dispersion relation is still given by k 2 = 0 for both modes.By standard variational theory of linear waves [20], this immediately yields Here, A + and A × are real functions quadratic in the wave amplitude (we define them as positive semi-definite for clarity), A .
, and the minus sign ensures that the wave energy be non-negative.Then, the corresponding action integral (70) becomes where √ −g no longer can be replaced with unity.Because S (2) and k 2 are scalars, and so is d 4 x √ −g, the functions A + , A × , and A are scalars.One can readily link them to the quantities that we have introduced above for the normal coordinates by comparing the Lagrangian densities ( 89) and (91): and the elements of a αβ can be inferred from A + and A × using that we already know the wave polarization (to the zeroth order in ǫ, which is of interest here) in the normal coordinates.Finally, we find that and can describe the evolution of the wave amplitude using Eqs.( 78) and (79).[The precise formulas for J α and T α β inhomogeneous medium contain small local corrections compared to Eq. ( 94), but those are bounded to remain O(ǫ) and vanish in the GO limit.]The former readily shows that A respond nontrivially to the background inhomogeneity: while Eq. ( 79) becomes ∇ β T α β = 0.
V. QUASILINEAR THEORY

A. Basic equations
Now let us consider the self-consistent evolution of g αβ as an independent variable.Equations (63) do not change in this case, but an additional equation emerges: which can also be written as where G αβ and T αβ are the background Einstein tensor and the background energy-momentum tensor respectively, given by substituting the full metric with the background metric in Eqs. ( 5).The equation above ( 97) is invariant with respect to the gauge transformations (42) by the invariance of S (2) .It is nonlinear in that N αβ is quadratic in h αβ and describes the slow evolution of the background metric in response to GWs.In contrast, direct interactions between GWs are neglected in this model, as seen from the linearity of Eq. ( 63).We call this model QL by analogy with the QL models that are used in plasma physics to describe the interactions between plasmas and electromagnetic waves [15,17,25,66].More specifically, it is an adiabatic QL model in the sense that g αβ and h αβ are assumed to be the only independent variables in the system (see also Sec.II B).

B. Alternative independent variables
The expression for S (2) in Eq. ( 98) can be taken from Eq. (57b).In this case, it may be more convenient to switch from independent variables [g, h, h * ] to the independent variables [g, h, h * ], where where the second equality is introduced by analogy with Eq. ( 54).Then, Eq. ( 58) can be expressed as so W depends only on h but not on g αβ .This shows that the functional derivative in Eq. ( 98) can be taken at fixed W : where [g] and [h, h * ] have been included in arguments to emphasize the dependence on the corresponding fields.Also, Eq. (63c) can be written as where we have introduced an alternative Hermitian dispersion operator C. Example: vacuum GWs As an example, let us consider quasimonochromatic vacuum GWs in the general smooth background, which we previously discussed in Sec.IV D. Let us rewrite Eq. (92) as where kα .= ∂ α θ and we have also introduced Ā .= (−g) 1/4 A as a new measure of the wave amplitude.By taking variation in Eq. ( 98) with respect to g αβ at fixed θ and A, one obtains where we used Eq. ( 90).Because T αβ = 0 in vacuum, Eq. (97) becomes, expectedly, Similarly, for a continuous spectrum, we have Here, W Ψ .= W + + W × , W s are the Wigner functions of the two vacuum modes: and s = +, ×.Hence, Eq. (97) can be expressed as Also note that the Wigner function is delta-shaped in this case, W Ψ ∝ δ(k 2 ), because the waves are constrained by the dispersion relation k 2 = 0.This means, in particular, that by taking the trace of Eq. ( 109), one finds that the corresponding Ricci scalar is zero.

VI. GAUGE INVARIANCE BEYOND THE QUASILINEAR APPROXIMATION
There is also an alternative path to a gauge-invariant theory of GWs, which extends beyond the QL approximation and is formulated as follows.Let us represent the perturbation metric as h αβ = aϑ αβ , where ϑ αβ = O(1).Then, one can search for g αβ in the form and a is now a constant that represents the characteristic GW amplitude.This implies because, by the previously stated assumption, h αβ = O(a).The variation σ αβ that enters the Einstein equations (4) can be represented as Then, using the Einstein equations σ αβ = 0, one obtains an infinite set of equations for ϑ These equations are invariant with respect to coordinate transformations (25), which is seen as follows.The functions σ Similar arguments apply when ξ µ scales with a nonlinearly, say, as ξ µ = n a n ξ µ (n) .In this case, Eq. ( 114) can be written (exactly) in the form where L(ξ)σ αβ denotes a term that is linear in σ αβ and, generally, nonlinear in ξ µ .The relations between σ ′(n) αβ and σ (n) αβ depend on specific ξ µ (n) , but one can still apply the same arguments as in the previous case.Then, one again finds that Eqs.(116) are equivalent to Eqs. (113).
Finally, notice the following.The functions σ αβ , is not the background metric g αβ .In order to find g αβ to the second order in a, one must solve the equation σ αβ and take the average as described in Eq. (111a).This makes our formulation different from that in Ref. [28], where gauge invariance was lost because the second-order corrections were introduced directly into the zeroth-order equation.

VII. SUMMARY
In summary, we show how to keep the theory of dispersive GWs gauge-invariant beyond the linear approximation and, in particular, obtain an unambiguous gauge-invariant expression for the energy-momentum of a GW in dispersive medium [Eq.( 79)].Using analytic tools from plasma physics, we propose an exactly gaugeinvariant QL theory, in which GWs are governed by linear equations (Sec.IV) but also affect the background metric on scales large compared to their wavelength (Sec.V).As a corollary, the gauge-invariant GO of linear dispersive GWs in a general background is formulated (Sec.IV B).As an example, we show how the well-known properties of vacuum GWs are naturally and concisely yielded by our theory in a manifestly gauge-invariant form (Secs. IV C and V C).We also show how the gauge invariance can be maintained within a given accuracy to an arbitrary order in the GW amplitude (Sec.VI).These results are intended to form a physically meaningful framework for studying dispersive GWs in matter.
This material is based upon the work supported by National Science Foundation under the grant No. PHY 1903130.
αβ only with i ≤ n.Hence, the set of equations (113) can be truncated by considering only 0 ≤ n ≤ m with some finite m.In particular, note that the equation with n = 1, which governs ϑ(1)αβ , is solved on the solution of σ