Gravitational-wave cosmological distances in scalar-tensor theories of gravity

We analyze the propagation of high-frequency gravitational waves (GW) in scalar-tensor theories of gravity, with the aim of examining properties of cosmological distances as inferred from GW measurements. By using symmetry principles, we first determine the most general structure of the GW linearized equations and of the GW energy momentum tensor, assuming that GW move with the speed of light. Modified gravity effects are encoded in a small number of parameters, and we study the conditions for ensuring graviton number conservation in our covariant set-up. We then apply our general findings to the case of GW propagating through a perturbed cosmological space-time, deriving the expressions for the GW luminosity distance $d_L^{({\rm GW})}$ and the GW angular distance $d_A^{({\rm GW})}$. We prove for the first time the validity of Etherington reciprocity law $d_L^{({\rm GW})}\,=\,(1+z)^2\,d_A^{({\rm GW})}$ for a perturbed universe within a scalar-tensor framework. We find that besides the GW luminosity distance, also the GW angular distance can be modified with respect to General Relativity. We discuss implications of this result for gravitational lensing, focussing on time-delays of lensed GW and lensed photons emitted simultaneously during a multimessenger event. We explicitly show how modified gravity effects compensate between different coefficients in the GW time-delay formula: lensed GW arrive at the same time as their lensed electromagnetic counterparts, in agreement with causality constraints.

for a perturbed universe within a scalar-tensor framework. We find that besides the GW luminosity distance, also the GW angular distance can be modified with respect to General Relativity. We discuss implications of this result for gravitational lensing, focussing on time-delays of lensed GW and lensed photons emitted simultaneously during a multimessenger event. We explicitly show how modified gravity effects compensate between different coefficients in the GW time-delay formula: lensed GW arrive at the same time as their lensed electromagnetic counterparts, in agreement with causality constraints.

Introduction
The propagation of gravitational waves (GW) through cosmological distances offer promising new avenues for testing cosmology. For example, information on the GW luminosity distance as extracted standard siren events can be used to probe the distance-redshift relation [1][2][3][4][5][6], leading to measurements of the present-day Hubble parameter using GW observations [7]. Also, gravitational wave measurements allow us to test deviations from General Relativity (GR): in fact, GW results have been recently applied for excluding modified gravity models predicting a speed of gravitational waves different than light [8][9][10][11], as earlier suggested in [12,13]. Measurements of the GW luminosity distance can also be used to probe modified gravity friction effects for GW travelling through cosmological backgrounds, see e.g. 45]. In the future, precision GW measurements at high redshift, provided by LISA [37][38][39][40] and the Einstein Telescope [41] will offer new possibilities for testing our understanding of gravity and cosmology. With this aim in mind, it is imperative to further theoretically characterize the propagation of GW in alternative theories of gravity, also taking into account the implications of cosmological inhomogeneities [42][43][44] that might influence or be degenerate with modified gravity effects. This is the scope of this work, concentrating on high-frequency scalar-tensor theories of gravity in the limit of geometric optics. We focus our analysis on propagation effects only, assuming that at emission the properties of GW is identical to General Relativity.
In section 2 we use of symmetry principles based on coordinate invariance for characterizing our scalar-tensor system and the behaviour of propagating degrees of freedom. We make use of a fully covariant formulation, spelling out in detail symmetry properties under coordinate transformations for each of the sectors involved. This allows us to carry on a general, model independent analysis of scalartensor systems, also identifying physically reasonable conditions for decoupling the evolution equations of different sectors.
Basing our considerations on symmetry principles based on coordinate transformations, we then derive in section 3 the most general structure of the high-frequency GW evolution equations and energymomentum tensor, for a scalar-tensor set-up in the limit of geometric optics. We also define a covariant condition to express graviton number conservation in our framework (see also [45]). Modified gravity effects factorize in front of our expressions, and the overall factor has a simple physical explanation in terms of the modifications of the linearized evolution equations.
We obtain a set of covariant equations that can be used in a variety of situations. In section 4 we apply them to the study of GW propagating through a perturbed cosmological space-time. In fact, distinguishing implications of modified gravity from effects of cosmological perturbations will be a crucial step for extracting physical information from future GW detections. Building and extending the classic results by Sasaki [46] (developed for studying propagation of photons in a perturbed cosmological universe within General Relativity) we derive the expressions for the GW luminosity distance d for a perturbed universe within a scalar-tensor framework, for scenarios where graviton number is conserved. Since this relation is at the basis for relating angular and luminosity distances in GW measurements, it is of crucial importance to prove its validity in a general theory of gravity for GW propagation on a general space-time. (See [47] for a recent work discussing probes of Etherington reciprocity law using GW measurements.) Given that GW luminosity distances can be modified with respect to GR, also angular distances can receive corrections.
Values of angular distances d (GW) A are important in phenomena involving strong lensing of GW, for example for the time delay of lensed GW. Strong lensing of GW can be important in the future for providing alternative ways for determining cosmological parameters (see e.g. [48]). When focussing on the limit of geometric optics for studying the propagation of GW and electromagnetic waves, since they both follow null-like geodesics we expect that GW and light arrive at the same time at the detector, if they are emitted at the same time [49,50]. In section 4.3 we explicitly show how to express the GW time-delay formula in terms of combinations of d (GW) A , in such a way that all effects of modified gravity compensate and one finds identical time-delays for GW and electromagnetic signals, if they are emitted simultaneously during a multimessenger event.
Our conclusions can be found in section 5, and are followed by six technical appendixes.

Our set-up
We develop a covariant approach for investigating the dynamics of high-frequency modes in scalar-tensor theories of dark energy. Symmetry arguments based on coordinate invariance allow us to determine general formulas describing the evolution of high frequency gravitational waves. Our set-up is described by a covariant action coupling gravity with a scalar field φ -the dark energy (DE) field -and with additional matter fields, schematically indicated with matter in action (1). We make the hypothesis that this action is invariant under diffeomorphism transformations, i.e. coordinate reparameterization invariance: x µ → x µ + ξ µ (x) for arbitrary infinitesimal vector ξ µ . We do not need to further specify the structure of the Lagrangian L for our arguments, but in what follows we assume that matter fields are minimally coupled with the metric g µν (possibly after performing appropriate conformal transformations to select a Jordan frame). The dark energy field φ, on the other hand, can have non-minimal kinetic couplings with the metric, that generally influence the propagation of GW. See e.g. [51] for a comprehensive review on modified gravity models including scenarios with non-minimal couplings of scalars with the metric. One of the delicate issues in studying GW propagation in modified gravity is to distinguish tensor from scalar fluctuations, and correctly identify their roles in the evolution equations of high-frequency fields. This topic started with the classic papers [52,53], and has been recently reconsidered in [54][55][56][57][58] using a variety of methods. The issue can be subtle in theories where scalar and metric fluctuations propagate with different speed, a phenomenon associated with spontaneous breaking of global Lorentz invariance by means of a non-vanishing time-like gradient for the dark energy field. Here we develop a covariant approach to address the problem, more similar in spirit to the original works of Isaacson and to the effective field theory of inflation [59] and dark energy [60] (see e.g. [61] for a comprehensive review). Our framework is distinct from ones based on decomposing graviton helicities in terms of their rotational properties with respect to the GW axis of propagation.

The perturbative expansion in high-frequency fields
We base our considerations on a double perturbative expansion for the metric and the scalar field around quantities solving the background equations, as [62,63]. Schematically, we expand metric and scalar fields as and we are interested to study the dynamics of the metric and scalar perturbations h µν and ϕ. In the previous expression fluctuations are distinguished from the background both for their absolute sizewe call it expansion in the amplitude, controlled by a parameter α -and for the size of their gradients -we call it expansion in gradients, controlled by a parameter . More specifically: -The α-expansion in the amplitude is used to define the so-called linear (first order) and quadratic (second order) approximations, and is common in cosmology. The parameter α is a book-keeping device to denote the order of amplitude expansion.
-The −expansion in gradients is controlled by the physical quantity controlling the ratio among the typical (small) wavelength λ of the high-frequency fields versus the (large) scale L B of spatial variation of slowly-varying background quantities. Among the latter, we include a dark energy scalarφ(x) whose time-like profile varies on scales of order L B .
The fluctuations h µν and ϕ are thought as high-frequency fluctuations whose gradients are enhanced by a factor of 1/ with respect to the background; moreover, they are small perturbations whose amplitude is suppressed by a factor of order O(α) with respect to the background. The possibility to use as small parameter to organize a perturbative expansion is one of the key observations of Isaacson: his approach is reviewed and expanded in the textbooks [64,65]. We adopt it here, extending the discussion of [55]. This framework allow us to implement a geometric optics limit where a generic small fluctuation σ(x) (scalar or metric) is decomposed into a slowly-varying amplitude, and a rapidly-varying phase (we understand the 'real part' symbol in what follows) is the small parameter of eq (4) controlling the rapid phase variations. The evolution equations contain up to second order derivatives in the fields: hence, substituting Ansatz (5) in such equations, we expect contributions scaling as 1/ 2 , 1/ , plus positive (or null) powers of . The geometric optics framework focusses on the leading (1/ 2 ) and next-to-leading (1/ ) orders in expansion in the small parameter -controlling respectively the evolution of phase and amplitude -and neglects the higher-order terms.

The symmetry transformations
We consider coordinate transformations acting on the quantities h µν and ϕ. We denote the gradient of the low-frequency scalar mode profile as spontaneously breaking global coordinate reparametrization along the direction of the time-like vector v µ . The vector v µ can be thought as being associated with cosmological acceleration, analogously to the approaches of the effective field theory for inflation and dark energy [59,60], and plays a special role in our discussion. From now on, all covariant derivatives are taken with respect to the lowfrequency metric fieldḡ µν of eq (2), used also to raise and lower indexes. The non-vanishing gradient (6) has important implications for diffeomorphism transformations. Under a change of coordinates the linearized fluctuations transform as for infinitesimal vector ξ µ . The scalar symmetry transformation (8) corresponds to a non-linearly realized diffeomorphism transformation, after the spontaneous space-time symmetry breaking associated with the scalar gradient v µ . We assume from now on that the amplitude of h µν is of order O( 0 ) in a gradient expansion, and we neglect in what follows possible contributions of order O( 1 ) and higher in the -parameter. (We checked that, even including those contributions, the arguments we develop are all still valid.) In order to actively apply the transformation on the fast moving modes h µν we assume that the size of ξ µ is reduced by a factor of with respect to h µν . i.e.
The gradients of ξ µ in eq (7) enhance its contributions by a factor O(1/ ), so that the result is of order of the same order of h µν in an -expansion. Again, for simplicity we assume that ∇ µ ξ ν does not receive contaminations at order O( 1 ), since as mentioned above we neglect contributions of order O( 1 ) and higher to the metric fluctuations h µν .
What can we say about the size of ϕ? We start noticing that the symmetry transformation (8) 'turns on' high-frequency scalar excitations even if they are initially absent. The dynamics of the two sectors, metric and DE perturbations, is inevitably coupled in their path from emission to detection. Even if DE fluctuations are not produced at the source (for example thanks to some screening mechanism), they can be generated by metric fluctuations that are travelling from source to detection. We then expect that propagation effects are able to excite scalar modes with an amplitude suppressed by a factor of O( ) with respect to metric fluctuations: Then, scalar modes transform non-trivially under the non-linearly realized diffeomorphism transformations controlled by the quantity v µ ξ µ (which is of order O( 1 )).
Our expectation encoded in the hierarchy (10) is also supported by interpreting scalar fluctuations ϕ as 'Goldstone bosons' of global space-time symmetries broken by the scalar profile v µ [59,60]. The size of the background gradient is of order ∇ µφ ∼ L −1 B . Keeping a fixed high-frequency wavelength λ for the metric fluctuations h µν , in the limit ∇ µφ → 0 (or equivalently L B → ∞) we expect the scalar Goldstone modes ϕ to be absent, since the symmetry is restored, and Goldstone bosons do not propagate. The size of the scalar excitation ϕ can be then expected to be suppressed by a factor λ/L B ∼ with respect to metric fluctuations, in agreement with eq (10). Motivated by these arguments, we impose the hierarchy (10) for linearized fluctuations.
Under these hypothesis, in the technical appendix A we build combinations of scalar and metric fluctuations that transform conveniently under coordinate transformations. After appropriate gauge fixings, we single out a transverse-traceless dynamical fluctuation h (T T ) µν of order O( 0 ), orthogonal to the vector v µ , that we identify with the high-frequency GW. Its dynamics is invariant under the residual transformation where ξ (T ) µ is a vector orthogonal to v µ , which satisfies additional gauge conditions we spell out in appendix A. The remaining high-frequency degrees of freedom are scalar modes. We develop arguments to show that, under the condition that scalar and GW propagate with different velocities along different geodesics, such scalar perturbations decouple from GW modes at linear order in perturbations. From now on, for definiteness, we then concentrate on studying the dynamics of the transverse-traceless GW modes h (T T ) µν , leaving the study of the independent scalar sector, when propagating 1 , to a separate work.

GW evolution equations and energy momentum tensor
Symmetry considerations provide a powerful tool for constraining the dynamics of our system. In fact, we can use symmetry arguments to determine the structure of the linearized GW equations of motion and energy-momentum tensor, with no need to rely on specific models. This is the aim of this section.

The linearized evolution equations
Isaacson, working in the context of the geometric optics limit of General Relativity (GR), shown that the original diffeomorphism invariance is preserved order-by-order in the gradient expansion, and at each order in the system is invariant under coordinate transformations [62,63]. This property further demonstrates the utility of the perturbative scheme based on gradients, which can be made compatible with the symmetries of the original theory, at least within the limits of geometric optics. We now make use of this fact in the scalar-tensor framework we are interested in. We change perspective and impose the symmetry invariance of the evolution equations at each order in the -expansion. As we will see, this viewpoint allows us to write the most general structure for the equations governing the GW dynamics, and to encode the effects of modified gravity in few physically transparent parameters.
Our starting point are the Einstein equations for high-frequency GW fluctuations, expanded at first order O(α 1 ) in the amplitude. Calling G µν the Einstein tensor, and with the suffix (n) the order of the amplitude α-expansion, the system of equations can be written as where h are the transverse-traceless fluctuations, orthogonal to the vector v µ , that we identify with GW -see the discussion around eq (11). We focus our attention to the leading and next-toleading orders −2 and −1 in the gradient expansion, which define the geometric optics framework as discussed after eq (5). Such contributions are obtained by singling out terms containing respectively second and first derivatives on the fields involved. As stated at the end of the previous section, we focus on the evolution of the GW tensor h (T T ) µν only, and postpone an analysis of dynamically independent high-frequency scalar modes to a separate work.
In the usual geometric optics Ansatz of GW propagation in GR, it is costumary to assume that matter fields are slowly-varying, and one considers the evolution equations at order 1/ 2 and 1/ as free equations (T µν = 0). Here we go beyond this hypothesis. In fact, a dark energy scalar field plays an important role in determining the behaviour of gravity at large cosmological distances. Kinetic couplings between scalar and metric lead to derivatives acting on the high-frequency modes, contributing to the effective linearized EMT T (1) µν . Non-minimal couplings between the dark-energy scalar and the metric are in fact common and well-motivated in theories of dark energy and modified gravity.
Nevertheless, symmetry considerations allow us to determine the general structure of T (1) µν , without relying on specific models. We demand that GW propagate with the speed of light. Since h (T T ) µν satisfies a transverse-traceless gauge, as well as the orthogonality requirement v ρ h (T T ) ρσ = 0, one finds that the left-hand-side of (12) reads We can now discuss the allowed structure for the linearized energy-momentum tensor T (T ) µν contributing to the GW evolution equation at orders −2 and −1 . It should be transverse-traceless, and orthogonal to v µ at orders −2 and −1 ; moreover, it should be conserved at order −2 : ∇ µ T (T ) µν −2 = 0, and it should be invariant under the transformation of eq (11). Finally, we demand that it ensures that GW propagate with light speed, to be consistent with GW170817 constraints [67], and have then a standard dispersion relation.
The only allowed structure of the linearized T

(T )
µν (h ρσ ) that satisfies these requirements at orders −2 , −1 is where τ A,B depend only on slowly varying fields. In fact, at order −2 , T (T ) µν contains second derivatives, but the unit-speed condition only allows for the combination proportional to τ A in the formula above. At order −1 it contains first derivatives, but the gauge conditions we impose allow only for the contribution proportional to τ B in eq (14). Calling the combination which depends on slowly-varying fields only, we rewrite the linearized evolution equation for GW fluctuations in terms of a single parameter characterizing deviations from GR: The deviations from GR on the propagation of high-frequency GW only appear as a first-order gradient of the GW high-frequency fluctuation, proportional to the parameter T depending on slowly-varying fields. Such contribution can be thought as a 'friction term' for the GW, and is common to find it in scalar-tensor systems with non-minimal kinetic couplings between scalar and metric degrees of freedom. In the context of gravitational wave cosmology in modified gravity several groups explored the consequences of such friction term in specific cosmological models, see e.g. [14, 16-24, 26-30, 45]. (see also [68] for a review), finding it is related with the parameter called α M in the effective field theory approach to dark energy. Also, this friction term arises in cosmological models with time-varying Planck mass: see Appendix B for the analysis of a representative model 2 . It is interesting to find that the fully covariant 'beyond-GR' friction term in eq (16) is the only one allowed by our symmetry principles and our physical considerations. In what follows, we consider the quantity T as an effective parameter controlling deviations from General Relativity.

Evolution equations in the limit of geometric optics
We now discuss how our covariant equations at leading and next-to-leading orders in an -gradient expansion allow us to derive the evolution equations for the physical degrees of freedom in the limit of geometric optics. As stated above, we focus only on the GW sector controlled by the transverse-traceless tensor h µν . The eikonal Ansatz for the GW reads with A T the amplitude, ψ (T ) the phase, the small parameter of eq (4), and e µν a polarization tensor normalized such that e µν e µν = 1. The gradient of the phase defines the GW 4-momentum 3 : We now apply Ansatz (17) to the covariant evolution equation (16), and separate the geometric optics analysis of the orders 1/ 2 and 1/ in our gradient expansion. The transverse-traceless condition, and the condition of orthogonality with respect to v µ impose the following requirements on the polarization tensor: where notice that there is a degeneracy among the last two conditions, so the previous equations provide 8 instead of 9 independent conditions. The order 1/ 2 of the equations, obtained from singling out second derivatives on the fields, control the evolution of the GW phase and the GW dispersion relations, leading to I.e. the GW 4-momentum is a null vector, propagating along a null geodesics. It is convenient to define the affine parameter λ controlling the evolution along the GW geodesics: for any function f , the derivative along the affine parameter is defined by The integral curves of the vectors k µ define the GW-rays: an important quantity for what follows.
While so far nothing changes with respect to General Relativity, at order 1/ -obtained from first derivative contributions -the effects of modified gravity become manifest. The evolution equation for the amplitude is where the quantity T , given in eq (15), depends on slowly-varying fields only. Recalling that v µ = ∇ µφ , the previous equation can be 'integrated' to The schematic expression T denotes the following integral with λ s corresponding to the value of the affine parameter at the source position. The quantity (25) represents a cumulative integration of modified gravity effects (the friction term in eq (16)) over the the GW geodesic's affine parameter. In integrating eq (23) we have chosen boundary conditions so that modified gravity contributions vanish at the location λ = λ s of the source, as expected since near emission modified propagation effects do not have time to develop. Modified gravity effects get exponentiated and appear as an overall factor inside the parenthesis in equation (24): importantly, we do not need to demand that T is 'small' for writing the equation. The exponential structure above will have several interesting consequences for our discussion.

The energy momentum of GW at second order in perturbations
Isaacson [63] proved that GW can be associated with their own energy-momentum-tensor (EMT), defined at second order in the α-expansion, which can influence the background dynamics. Schematically, we can write where T (2) µν denotes the EMT associated with GW. As stated above, we focus on the contributions associated with the transverse-traceless tensor fluctuations h (T T ) µν only, and do not discuss scalar contributions in this work since, under our hypothesis, the two sectors evolve independently. The quadratic terms in the GW energy-momentum-tensor have equal-size momenta in opposite directions which compensate each other, hence contributing at zeroth order in the -expansion. Using only symmetry arguments we are able to determine the structure of the GW contribution to the tensor T (2) µν in a general class of scalar-tensor systems.
When focussing on transverse-traceless excitations, Isaacson's result for the EMT is The symbol . . . denotes the so-called Brill-Hartle spatial average, see [63,64]. Among other things, this average procedure ensures that the EMT is diffeomorphism-invariant, and conserved. Interestingly, the condition of coordinate invariance fixes the structure of T µν associated with GW. In fact, recall the EMT is quadratic in h (T T ) µν , and contains two derivatives in total acting on the transverse-traceless GW excitations (by 'integration by parts', we can place one derivative per field). The structure in the combination (27) within the average is the only one with these properties, and that is compatible with the condition of invariance under symmetry (11). See the discussion in Appendix C.
The only freedom we are left with is in the overall factor in front of the Brill-Hartle average appearing in eq (27). In fact we can change perspective, and use the condition of invariance under symmetry for determining the structure of T (2) µν in the scalar-tensor set-up we are interested in. In other words, we do not compute the EMT using a 'top-down' approach starting from a given theory, but instead we deduce its structure from the symmetry conditions imposed in the theory. In the scalar-tensor framework we are focussing on, the previous considerations allow for the following structure for the EMT in modified gravity where # is a function (to be determined) of the slowly-varying fields, metric and scalar. We now proceed to determine this quantity, making use of the condition that the energy-momentum tensor T (2), MG µν should be conserved by virtue of the Bianchi identity, and of the geometric optics evolution equations of section 3.2. We substitute the geometric Ansatz of section 3.2 to the previous formula, and get Using the evolution equation (24), as well as the condition (20) that GW follow null-like geodesics, the condition of conservation of the EMT fixes # to the value e − T as defined 4 in eqs (24), (25): Hence we find that the second order GW energy-momentum tensor in our scalar-tensor framework, in the geometric optics limit, reads These general considerations then allow us to single out transparently the effects of modified gravity in the overall factor depending on the quantity T of eq (25), a cumulative integral of modified gravity contributions along the GW geodesics from source to detection. As we will learn in what follows, phenomenological implications of our results, as well as the explicit example discussed in section B, further support the structure (31) for the EMT in the scalar-tensor systems under consideration.

Conservation of graviton number
We can do some further steps following [72], and relate the properties of the quantities above with graviton number conservation. we expect graviton number to be conserved within a GW ray bundle, and we are going to prove this fact in our setting within geometric optics. We express T (2), ST µν in terms of quantities k µ -interpreted as graviton 4-momentum -and N µ , as: where N µ is interpreted as the graviton number density, and is defined as Graviton-number conservation ∇ µ N µ = 0 is ensured by relation (24) 5 . See also [14] for a perspective on graviton number conservation in a cosmological setting in a modified gravity framework.
Moreover, by making use of well-known geometric relations (see Figure 3), we can express this condition in a geometrically more direct way, which further supports our identification of N µ with graviton number density. We call S(λ) the cross-sectional area of a GW bundle, and λ the affine parameter along each GW ray (the graviton trajectory) with four-momentum k µ . A geometric optics theorem (see [64], exercise 22.13) states that with λ the affine parameter associated with the GW 4-momentum k µ . Together with (33), relation (34) implies the important identity that makes more manifest the required flux conservation for a stream of gravitons crossing the S-areas along the GW evolution parameterized with λ. Notice the presence of the overall exponential factor e − T due to modified gravity -see the discussion after eq (25) -that changes the flux of gravitons through a given surface. Such coefficient plays the role of 'damping term' in the GW amplitude during propagation in a modified gravity set-up, as expected given that we interpret T as friction in the evolution equations. The result (35) will be important for cosmological applications in what follows.

Cosmological distances and GWs
We now apply the general findings of the previous sections to GW propagating through a perturbed Friedmann-Robertson-Walker (FRW) space-time. We prove the validity of Etherington reciprocity law between GW luminosity and angular distances in the scalar-tensor framework developed in the previous sections, and we discuss the implications of our findings for GW lensing.
Cosmologists use various different definitions of distance depending on the context, and the observables they are interested in (see e.g. [73,74] for enlightening reviews). While usually definitions make use of light detected from distant sources, GW offer new tools for measuring cosmological distances. We consider here two distinct GW distance probes: is defined in terms of the ratio of GW power emitted at source position (intrinsic GW luminosity), versus the GW flux at detector location -see section 4.1. The luminosity distance depends on the universe expansion rate, and enters into the GW waveforms and can be directly measurable by detecting GW from distant sources. Following early important works [1][2][3][4][5][6], d (GW) L is being recognized as a key observable to independently measure cosmological parameters by means of GW, as well as testing theories of modified gravity (see e.g. the review [68]).
The angular diameter distance d A is usually understood as being related with the luminosity distance d L through the so-called duality-distance relation, or Etherington reciprocity law d L = (1 + z) 2 d A . On the other hand, the theoretical validity of this relation should be explicitly proved, and this is of the aims of this section, together with applications to GW lensing. In fact, since we learned in the previous section that GW evolution is affected by the friction term proportional to T in eq (23), we expect that both luminosity and angular distances are influenced by modified gravity. We show that, thanks to graviton number conservation (see section 3), these quantities are related by Etherington reciprocity law (see e.g. [72,106] for the case of photon propagation) for GW propagating through perturbed FRW space-times in the scalar-tensor scenarios we are focussing on.
Our treatment in this section follows very closely the discussion of the classic paper by Sasaki [46] that for the first time discussed the concept of luminosity and angular distances for photons propagating in a perturbed FRW universe. Sasaki's early work was followed by many articles that further generalized it extending the analysis of luminosity distance for photons in a perturbed background -see e.g. [42][43][44][107][108][109][110]. Nevertheless, as we are going to discuss, the formalism developed in [46] is sufficiently flexible to be applied to GW propagation on our scalar-tensor systems, with little adaptations needed along the way.
We start by introducing some geometric tools we need for our arguments. The space-time metric we are interested in is written as a conformally flat FRW universe, and reads where from now on in this section we denote with a hatĝ µν the physical space-time metric, while g µν is the 'comoving' part of the metric tensor. Analogously to eq (2), we can write and h µν the high-frequency field. In this cosmological context, the comoving metric componentsḡ µν corresponds to Minkowski space-time, plus the long wavelength perturbations. Gravitational waves correspond to transverse-traceless metric fluctuations, whose null-like 4-momentum propagates along null geodesics. The physical GW energy momentum tensor is given by eq (31): where the overall factor depending on T controls the modified gravity contribution, and the GW physical 4-momentum isk withλ the physical affine parameter along the GW ray, as associated with the vectork µ . Following the procedure developed in section 3.2, we conveniently express the evolution eq for the GW amplitude whereθ is the expansion parameter along the GW geodesics.
It is convenient to pass from physical (hat) to comoving (no hat) quantities. A conformal transformation maps a null GW geodesics inĝ µν into null geodesics in g µν [111]. The GW affine parameter scales as dλ → dλ = a −2 dλ , with a the scale factor. The evolution equation in the comoving frame then results The evolution equation for the comoving expansion parameter θ can be easily determined [46], finding where σ is the shear along the GW geodesics: and R µν the space-time Ricci tensor. The graviton number conservation (35) remains unchanged and reads d dλ in this comoving frame (λ referring to the comoving affine parameter). These are the geometric ingredients we need for the physical considerations we develop next.

The GW luminosity distance
To determine the GW luminosity distance we proceed step-by-step as [46]. We introduce an observer whose physical four-velocity we denote withû µ . The measured GW energy flux by such observer readŝ whereT µ ν is the GW energy momentum tensor (39), and while the GW flux amplitude and GW frequency measured by the observer are The notion of GW frequency allows us to define the GW redshift z at the value λ of the comoving GW geodesics affine parameter λ For computing the luminosity distance, we assume that GW are emitted by an approximately spherically symmetric system, with characteristic radius R s (this assumption is nevertheless not important since we send R s → 0 at the end of the calculation). The flux amplitude F measured at the source position is related with the intrinsic source luminosity by the relation with λ s the conformal affine parameter at the source. See Fig 2, left panel. The luminosity distance to the source as measured by an observer located at λ = 0 is defined as Substituting relation (50), we find the following expression Notice the role of the modified gravity friction term in the overall exponential factor, containing the cumulative integral of the friction parameter T along the GW geodesics path (see the discussion after eq (25)). Before proceeding, to make contact with the literature, it is interesting to consider the ratio between the GW luminosity distance (55) versus the electromagnetic luminosity distance. This ratio provides an interesting observable in case of multimessenger events. We find singling out the modified gravity contribution as an integral from λ = 0 (the position of the observer) to the source at λ = λ S . Once substituting in the integrand of (56) the explicit form of friction terms used for parameterizing deviations from GR in cosmological models of dark energy, one finds exactly the same formulas used in the literature -see Appendix D for a discussion of such comparison.
We now proceed expressing the GW luminosity distance (55) in an alternative way, that is more useful for explicitly including effects of cosmological perturbations, and for then comparing with the angular distance in section 4.2. From now on, we denote a a perturbed quantity with a tilde, and unperturbed without tilde. For example, we write for the comoving metric meaning that δg µν are long-wavelength perturbations. We now introduce a null vectorK µ proportional tok µ , and use it to define a corresponding affine parameterization: where λ from now on is the affine parameter associated withK µ . This vector is normalized in such a way that, once evaluated at the source λ = λ s , we find whereũ µ is the perturbed comoving observer 4-velocity.
The introduction of the vectorK µ . is technically convenient to easily relate the physical size of the source with the affine parameter along the GW geodesics. In fact, as shown in [46], the characteristic size R s of the source can be expressed as with ∆λ s the infinitesimal affine parameter associated with the source size. See Fig 2, where the suffix L is included to associate the expansion parameter θ with the luminosity distance. The deviation for the expression of θ L at first order in cosmological inhomogeneities can be expressed as The evolution equation for the first order perturbation δθ L (λ) can be obtained from from (43). It reads where R µν is the perturbed space-time Ricci tensor at the position λ along the GW geodesics. Integrating eq (63) along the affine parameter, imposing the boundary condition δθ L (λ s ) = 0, we get Collecting all the results so far, we can integrate eq (42), and get the relation This result can be inserted into eq (55): taking ∆λ s → 0 (i.e. considering a source of negligible size) we end with the compact expressioñ Notice that all the effects of modified gravity friction term are implicitly included in the expression (65), which relates the affine parameter λ s with the remaining quantities. The compact expression (66) (accompanied by relation (65)) is exact and include the effects of cosmological fluctuations -on the other hand is implicitly expressed in terms of λ s , and is not easy from it to extract in a physically transparent way the implications of cosmological fluctuations and of modified gravity. Such implications are more easily studied by using the cosmic ruler formalismsee [55]. This approach explicitly identifies contributions from peculiar velocities, weak lensing, Sachs-Wolfe effects, volume effects, and Shapiro time delay, and allows to appreciate the contributions due to modified gravity. We refer the reader to [55] for more details: for our purposes to prove the validity of Etherington reciprocity law -our aim for the next section -formula (66) will be sufficient.

The GW angular distance, and Etherington reciprocity law
We now prove the validity of Etherington reciprocity law connecting luminosity and angular GW distances. This relation is expected to hold in scenarios where graviton number is conserved, as our scalar-tensor set-up (see section 3.4).
The GW angular distance d GW A is formally defined in terms of the ratio between the angular diameter d s of the source located at conformal affine parameter λ s , and the source apparent angular size ∆φ as measured by an observer at λ = 0. In formulas: Following [46], it is convenient to re-express d with S(λ) the cross-section area of GW rays at λ, and its diameter by d(λ). ∆λ is the affine parameter in proximity of the observer. [46] proved the relation connecting the ratio d(∆λ)/∆φ with ∆λ. See Fig 2, right panel. Fig 2 shows that in evaluating the angular distance one considers GW bundles expanding from the observer position (while, on the contrary, the luminosity distance considers bundles expanding from the source). Hence, the expansion parameter θ A associated with angular distance reads, when neglecting effects of cosmological inhomogenities, can be obtained integrating eq (43): When including the contributions of perturbations, we find the formal solution for the first order perturbation δθ A to the angular expansion parameter.
Integrating eq (45) along the GW geodesics, and comparing with the definition (68), we get the relation Moreover, integrating eq (42) using the result (70), we now obtain as the relation between affine parameter and angular expansion parameter. Substituting the results of eqs (73) and (69) into eq (72), we obtain the expressioñ Comparing with the expression for the luminosity distance, we get The second line, eq (76), is the desired Etherington relation, valid including first order perturbations. (The step between eq (75) and (76) requires technical calculations that we defer to Appendix E.) Hence we proved that in the scalar-tensor framework discussed in this work, with conservation of graviton number, luminosity and angular distances for GW are connected by the classic Etherington law (76). We have seen that GW and electromagnetic luminosity distances can differ -see eq (56)and this fact is important in case of multimessenger events. Then eq (76) tells us that the same is true for angular distances, and we can schematically write a relation analog to eq (56): In what comes next we briefly discuss some applications of these results to GW lensing.

Implications for GW lensing
Strong GW lensing from large-scale structures between GW source and detector is an important phenomenon that -although not yet observed -is likely to offer new ways to probe cosmological parameters with future gravitational wave detections. For example LISA, by observing sources from high-redshift sources, will likely detect lensed events [92]. See e.g. [49, 50, 58, 76-92, 92-105, 119] for works discussing this topic. We consider strong GW lensing from point-like lenses in the geometric optics limit, valid when the GW wavelength is well shorter than the Schwarzschild radius of the lens. In this limit, we do not need to discuss interference effects that, although very interesting, go beyond the scope of this work. We focus on the specific observable associated with the time-delay that the presence of the lens induces on the propagation time of the GW from source to detector. We compare the GW time-delay induced by the presence of the lens with the electromagnetic (EM) time delay of lensed light received in a multimessenger detection.
The works [49,50] shown conclusively that GW and EM lensed signals arrive at the same time at the detector, if both waves propagate at the same speed and are emitted at the same time. In the geometric optics limit this is expected when photons and GW travel through null geodesics, since by definition both sectors cover the minimal possible distance from source to detector. Causality arguments based on Fermat principle allow one to prove this statement in full generality. [49] also argues that the same result should be valid in any theory of gravity, to respect causality.
Said this, it is interesting to analyze the topic in an explicit modified gravity set-up, for understanding how effects of modified gravity balance so to ensure the same time delay for GW and light. This is the scope of this section. We find this topic interesting since the expression for the time-delay commonly used in the literature (see e.g. [72], as well as the recent [103] in the context of gravitational waves) explicitly contains factors depending on the angular distance d (GW) A , which can be modified with respect to the standard case (see eq (77)). In fact, the GW time delay ∆t (GW) has a geometrical contribution, and a Shapiro contribution t (GW) Φ due to the presence of inhomogeneities in the background cosmological space-time crossed by GW in their path from source to detection. We express it in the following form In the previous expression, z is the redshift, d (GW) OL the GW angular distance as measured from the observer to the lens, d is the aforementioned Shapiro contribution due to inhomogeneities. θ is the observed angular position of the source, θ S the would-be angular position of the source in absence of the lens. In the electromagnetic case, the corresponding time-delay ∆t (EM) has exactly the same structure, changing the suffixes from GW to EM. See Appendix F for a derivation of the geometric part (the first term) of the previous formula.
In [55] we shown explicitly that the Shapiro time-delay t Φ is exactly the same for GW and EM observations in a scalar-tensor framework, t : we refer the reader to this work for full details. However, eq (78) also contains explicitly angular distances in the first geometric term, which we dub ∆t geo . Using eq (77) we can in fact understand how the effects of modified gravity in the geometrical time-delay ∆t geo compensate. We can write: To pass from the second to third line we used relation (77), while the fourth line is a simple consequence of summing the integrals in the exponent. So we find that modified gravity contributions carefully compensate, leading to the equality ∆t geo . Together with the fact that the Shapiro time delay coincides in the GW and EM sector, we find that by using the expression eq (78) we are ensured that GW and EM experience the same time-delay from strong lensing in a scalar-tensor theory of gravity, as we wish to demonstrate. In particular, even if formula (78) is expressed in terms of angular distances that, when taken individually, can be different with respect to General Relativity, one finds cancelations between modified gravity coefficients that lead to the result (79).

Conclusions
In this work we studied the propagation of high-frequency gravitational waves (GW) in scalar-tensor theories of gravity, with the aim of examining properties of cosmological distances as inferred from GW measurements. We first developed a covariant set-up for our scalar-tensor systems, which by hypothesis are characterized by symmetry properties based on coordinate invariance. Symmetry considerations allowed us to extract transverse-traceless components of the high-frequency scalar-tensor fluctuations, identified with GW. In scenarios where scalar and tensor components propagate with different speeds the two sectors decouple at the linearized level around an arbitrary background, and the evolution of high-frequency GW and scalar modes can be studied independently.
We then determined the most general structure of the GW linearized equations and of the GW energy momentum tensor, assuming that GW move with the speed of light. Modified gravity effects are encoded in a small number of parameters, and we studied the conditions for ensuring graviton number conservation in our set-up. We then applied our general findings to the case of GW propagating through a perturbed cosmological space-time, deriving the expressions for the GW luminosity distance d . Both luminosity and angular distances can be modified with respect to General Relativity. We proved for the first time the validity of Etherington reciprocity law d in a perturbed universe within a scalar-tensor framework. We discussed implications of this result for gravitational lensing, focussing on time-delays of lensed GW and lensed photons emitted simultaneously in a multimessenger event. We explicitly wrote an expression for the time-delay formula, showing how modified gravity effects carefully compensate between different contributions. As a result, lensed GW arrive at the same time as their lensed electromagnetic counterpart, in agreement with causality constraints.
It would be interesting to consider scenarios where graviton number is not conserved, as done in alternative theories of cosmology where the photon number is not conserved (see e.g. [112]). Another avenue for future research is to extend our arguments to scenarios with non-standard dispersion relations where the speed of tensor modes is not equal to one, at least in some frequency ranges. More in general, it would be interesting to explore in more details how observables depending on the GW angular distance and GW lensing can be used to probe alternative theories of gravity. We plan to study these topics in future works.

Acknowledgments
It is a pleasure to thank Carmelita Carbone for input and discussions.

A A covariant approach to high-frequency fluctuations
In the main text we derived the covariant evolution equations and the energy-momentum-tensor for transverse-traceless GW excitations. In this appendix we spell out the technical arguments used for identifying the GW sector, and for distinguishing it from the scalar sector. We work within the covariant approach of section 2.

A.1 Decomposing the gauge transformations
Assuming a time-like 6 direction v µ for the scalar gradient of eq (6), we introduce the vector where We decompose the gauge vector ξ µ in eqs (7), (8) in a part orthogonal, and a part parallel to v µ : This decomposition defines what we call a T -gauge transformation, proportional to the vector ξ (T ) µ orthogonal to v µ , and an a S-gauge transformation, depending on the scalar ξ (S) . The metric and scalar field perturbations transform under a gauge transformation as Hence, ϕ transform only under an S-gauge transformations, while h µν both S and T -gauge transformations. We now show how to use the decomposition (82) to consistently distinguish GW from scalar excitations in model-independent framework. We introduce the quantitỹ Each of the contributions toh µν are of the same order in the gradient expansion , since we are assuming O(ϕ) ∼ O(h). The field combination in eq.(85) is S-gauge invariant: We define the orthogonal projection operator relative to the vector X µ such that X µ Λ µν = 0, and we apply it toh µν . We find The quantities Therefore, using eq. (86), we find that, under a T -gauge transformation and up to O( 0 ), the quantities h (S, V, T ) ...

transform as
thus, h (S) is also T -gauge invariant at order O( 0 ).

A.2 Gauge fixing
Since we demand that the linearized equations are invariant under coordinate transformations -i.e. separately S and T -gauge invariant -we can assume they can be organized in terms of the S-gauge invariant combinations h (S, V, T ) ...
. After identifying the gauge-invariant quantities, we can now make use of the T -gauge freedom of eqs (92) and (93) for fixing convenient gauge conditions to study the physics of the system.
The first gauge fixing condition we impose is by choosing ξ in eq (92). This condition is compatible, at order O( 0 ), with the orthogonality requirement X µ h Eq. (94) leaves the residual T -gauge freedom We separate h µν into a traceless plus trace components. Since X µ h (T ) µν = 0 we select the trace in the subspace orthogonal to X µ , with The quantity γ µν satisfies Λ µν γ µν =ḡ µν γ µν = 0 (99) since X µ γ µν = X µ h (T ) µν − 1 3 Λ µν h (tr) = 0. Starting from eq. (93), we can find how γ µν and h (tr) transform under a gauge transformation. For those vector fields ξ (T ) µ that satisfy eq (96) we have We now use eq. (101) to impose the transversality condition∇ µ γ µν = 0: This gauge transformation is valid up to order O( −1 ). We retain only contributions up to this order in : this is consistent with keeping only terms up to O( 0 ) in the gauge transformation of γ µν , since we are gauge fixing its gradient. After such gauge choices, the quantity γ µν is transverse and traceless. We dub it and we identify it as the high-frequency GW discussed in the main text. At this stage, we point out that is not possible to choose h (tr) = 0, within the residual gauge freedom given by (96), if h (tr) depends on the coordinate in the direction of X µ . For simplicity, we can exhaust the gauge freedom imposing ∇ µ ξ (T ) µ = 0, such that the trace h (tr) is gauge-invariant, while the transverse-traceless GW excitations h (T T ) µν are invariant under the residual transformation that can be read from eq (101): To sum up, after imposing the gauge conditions above, the S-gauge invariant metric perturbations (85) we started with readh We can also use the S-gauge for setting the unitary gauge ϕ = 0: then, after these gauge fixings,h µν coincides with the original metric fluctuations h µν .
Some words on the number of propagating degrees of freedom (dof). The quantityh µν , before we make any gauge choice, has 10 non-vanishing components, each of them a potential dof. Making gauge fixings as explained above we imposed 6 conditions, since both h

A.3 Separating the evolution equations
We now discuss the evolution equations for the metric perturbation of eq. (105), as obtained from the linearized Einstein equations. We build arguments to show that at the linearized level, under physically reasonable conditions, different sectors evolve independently one from the other.
Since by hypothesis our system is invariant under gauge diffeomorphisms, i.e. under coordinate transformations, Einstein equations can be expressed in terms of the fields discussed above, obtained after gauge-fixing appropriate gauge-invariant quantities. Therefore, in an unitary gauge with ϕ = 0, we write the linearized Einstein equations as withh ρσ the S-gauge invariant combination given in eq (105). Taking the trace of the previous linear equation we eliminate from the left-hand-side the dependence on the transverse-traceless fluctuation h (T T ) µν . In fact, the linearized Ricci scalar reads where Λ µν is the projector introduced in eq (87). We notice that while the trace scalar h (tr) receives a kinetic contribution controlled by the d'Alembertian operator , second derivatives acting on the scalar h (S) are always weighted by the projector operator Λ µν , hence they are directed on the space orthogonal to the vector v µ . Let us now use the time-like vector v µ for slicing the space-time into a family of space-like surfaces, as in the ADM approach to General Relativity (see e.g. [64]). As a result, one finds that Ricci-scalar derivative contributions to the evolution equations for h (S) are not sufficient for propagating this field: its dynamics is in fact constrained to live on the hypersurfaces orthogonal to v µ , with no components on the direction of the system evolution. Indeed, h (S) plays a role analogous to the lapse function N in the ADM formalism. This can also be deduced by the definition of h (S) in eq (89): this quantity collects the contribution to h µν from the components along the vector v µ , exactly as the lapse constraint in ADM. (We also discuss in Appendix B an explicit, simple example where h (S) is manifestly non-dynamical.) Can the energy-momentum tensor in the right-hand-side of eq (106), or the Ricci tensor in its left-hand-side qualitatively change these considerations? Not if it is derived from a covariant action as (1), where non-minimal couplings of dark energy scalar to the metric are expressed in a covariant form in terms of the metric, Riemann, and Ricci tensors. This is always possible in manifestly covariant and diffeomorphism invariant formulations of scalar-tensor systems. They contribute to the scalar kinetic terms with a structure as the one we obtained from the Ricci scalar (107): the second derivatives acting on h (S) appear in combinations weighted by the projector Λ µν as the parenthesis in the last term of eq (107). Hence the same considerations as above apply. Since for such space-time slicing it is a constrained field, h (S) does not propagate, and its role is to impose conditions on the slow-frequency part of the system, or on the remaining high-frequency modes.
Let us then assume to solve the equations of motion for the non-dynamical field h (S) , and substitute its solution on the original action. We are left with an action containing h (T T ) µν and h (tr) as potentially propagating high-frequency degrees of freedom. Thanks to linearity, the linearized Einstein equations for the propagating modes can be decomposed as We expect that second derivatives contributions on the scalar sector have a rich structure, with different coefficients in front of contributions orthogonal or parallel to the vector v µ , associated with spontaneous breaking of Lorentz invariance associated with the vector v µ . As a consequence, tensor and scalar fluctuations normally propagate with different velocities 7 . Given the strong experimental bounds on the GW velocity associated with the GW170817 event, we set the speed of GW to the one of light. Within this hypothesis, we implement a geometric optic Ansatz to both the h (T T ) µν and h (tr) sectors, using the geometric optic approach explained in the main text. We write: The amplitudes of both modes are slowly varying, while the phases are rapidly varying thanks to the factors of 1/ in the exponent. When plugging Ansatz (109) and (110) into eq (108), one gets a linear combination of terms with rapidly oscillating phases and slowly varying overall coefficients. Schematically, we expect that the geometric optics limit of Einstein equations has a structure as where within the parenthesis we collect slowly varying contributions at order −2 and −1 in a gradient expansion. The −2 contributions depend on derivative of the phases ψ (T T ) and ψ (tr) : they control the dispersion relations for the two species of excitations, scalar and GW (see section 3.2 for a geometric optics analysis of the GW sector). Since in general h and h (tr) propagate with different speed, they are characterized by distinct dispersion relations, hence the phases ψ (T T ) and ψ (tr) are different. Equation (111) is a linear combination of two contributions weighted by two distinct phases which rapidly oscillate over space and time: in order to satisfy it, we need to impose that the coefficients of each of these two terms separately vanish. Within the geometric optics limit, this procedure effectively separates the evolution of scalar modes (characterized by the phase ψ (tr) ) and GW modes (characterized by the phase ψ (T T ) ).
Given these considerations, in dark energy scenarios where the scalar and GW modes have different phases due to different dispersion relations, we can effectively separate the high-frequency GW and scalar sectors in eq (108), and write Within our hypothesis, GW sector is decoupled from the scalar sector at the linearized level, and we can study its dynamics as done in the main text. The arguments discussed above might be made more rigorous with a more systematic and detailed analysis of perturbations evolution equations, for example using the approach of the recent work [44]. We leave this analysis to separate investigations.
B A simple example: F (φ) R Let us make a specific, simple example of the friction-term contributions found in our general formula of eq (16), which arises in models characterized by a time-varying Planck mass controlled by the dark energy scalar field φ. We consider the following non-minimal kinetic coupling between scalar φ and metric which can be considered a part of the classic Brans-Dicke action [113]. We are not interested to study in detail the system, but only apply to its corresponding GW evolution equations the approach explained in Appendix A.
We decomponse the linearized Einstein equations in terms of the high-energy fluctuations, and focus on orders 1/ 2 and 1/ in a gradient expansion, as described in the previous section A. We find that GW modes obey the equation An evolution equation governing scalar modes can be determined by taking the trace of the Einstein equations where the vector X µ is defined in eq (80), and the projector Λ µν in eq (87). These equations have the structure expected from our considerations in the main text and in Appendix A. In fact, comparing the GW evolution equation (115) with the general expression in eq (16), we notice that the former has a friction term T controlled by the derivative of F along the dark energy field: T = 2F ,φ /F . Using the results of section 3.3, we find that the energy-momentum-tensor at second order in the transverse-traceless fluctuations reads (we choose the extreme of integration φ in such that F (φ in ) = 1) which is the expected structure associated with the Lagrangian of eq (114). The scalar modes have a kinetic structure depending on the Lorentz violating vector v µ . The kinetic term for h (S) in eq (116) is projected by the tensor Λ µν in the direction orthogonal to v µ , as discussed in Appendix A, hence it does not contribute to the dynamical degrees of freedom.
We can apply these findings to cosmology, and consider the case of GW propagating through a conformally flat Friedmann-Robertson-Walker (FRW) universe, with metric ds 2 = a 2 (η) η µν dx µ dx ν , and for a homogeneous scalar fieldφ =φ(η). Eq (115) results then (H = a /a, with prime denoting derivative along time) The effect of the friction term due to the non-minimal scalar-tensor couplings has the expected structure and is manifest within the parenthesis of the previous expression. It would be interesting to further apply our covariant approach to fluctuations to more complex scenarios, from Horndeski [114] to beyond Horndeski [115] and DHOST [116][117][118].

C Gauge invariance of the GW energy-momentum tensor
We show that the structure of eq (28) is fixed by the gauge invariance. The most general structure for the energy momentum tensor quadratic in the high-frequency transverse-traceless (TT) modes h where C αβγδ depends on slowly-varying fields. This is the most general structure for T (2), MG µν compatible with the TT gauge imposed on h (T T ) µν . After fixing the TT gauge as discussed in section A, we are left with invariance under the transformation in eq (11), as derived around eq (104). We need to ensure that T (2), MG µν is invariant as well. Notice that the tensor C αβγδ has to be invariant under an interchange of α and β. Under such gauge transformation, at the linearized level we get a contribution which must vanish for any ξ (T ) β , and for any choice of µ, ν. The only way to ensure this is to use the transverse gauge condition, and require that the contraction with C αβγδ forces the condition α = γ (or alternatively α = δ). This implies C αβγδ = C δ αγ C βδ (122) for some constant C. Plugging this result in eq (119) we get Since the EMT is symmetric in the indexes, the quantity C βδ is symmetric. Applying the transformation (11), we find the only non-vanishing contribution and the only way to make it always vanishing is to have C βδ ∝ δ βδ . Symmetry arguments force the EMT to be proportional to Isaacson's form (28), up to the overall constant, as we wish to prove.

D Comparison with the literature
We now that our expression (56) for the ratio among GW and electromagnetic luminosity distances: coincides with analog expressions found in the literature, once we specialize to a cosmological setting. Let us then analyze GW travelling through an unperturbed, conformally flat Friedmann-Robertson-Walker (FRW) universe with metric ds 2 = a 2 (η) −dη 2 + d x 2 . The dark-energy background scalar is time-dependent only,φ = φ(η), and we have v µ = (φ , 0, 0, 0). Let us use, for definiteness, the notation of [39]. We understand tensorial indexes and call h (T T ) µν = h(η, x). The evolution equation for the GW modes is expressed as (H = a /a) where modified gravity friction contributions are contained in the time-dependent parameter δ(η). The expression for the ratio between GW and electromagnetic luminosity distances is found to be (see []) We now apply our formula (56) to this cosmological example. Equation (16) reads Comparing with eq (127), we can then identify δ = Tφ /(2H). Using also the fact that dη/dt = 1/a, H = H/a, 1 + z = a(0)/a(t), we can write eq (128) as (the suffix 0 indicate present-day observers) which coincides with eq (126).
E Proof of the step between eq (75) and eq (76) We prove the validity of the step between eq (75) and eq (76), using results from [46]. For shortening the notation, we call q(λ) = δ(R µν K µ K ν ) λ the geometrical combination depending on the background geometry, which implicitly appears in both formulas (75) and eq (76). We call We wish to prove I L + I A = 0 .
First we re-express I L more conveniently. We change variable λ = λ s − σ, λ = λ s − σ . Then Then by an integration by parts we can re-write I A as a single integral as Moreover, simple steps lead to as desired.

F The geometric time-delay
We derive the expression for the geometric time delay of waves whose trajectory is bended by a pointlike lens, in the limit of geometric optics and of Euclidean geometry. The GW moves with the speed of light, and we consider (we should take care at the position of the indexes in the angular distances). We work in the limit of infinitesimal angles, so we can expand trigonometric functions. Hence Similarly, one has OB = D OL . We compute step by step the GW time delay, corresponding to the quantity (recall that GW travel at the speed of light, set to one) Expanding the cosine for small angles, we can reassemble the previous formula as Then the quantity we are after is which is the formula used in eq (78) of the main text.