Gravitationalwave tests of Gravity and Dark Energy
Dark Energy in light of MultiMessenger GravitationalWave astronomy
Abstract
Gravitational waves (GWs) provide a new tool to probe the nature of dark energy (DE) and the fundamental properties of gravity. We review the different ways in which GWs can be used to test gravity and models for latetime cosmic acceleration. Lagrangianbased gravitational theories beyond general relativity (GR) are classified into those breaking fundamental assumptions, containing additional fields and massive graviton(s). In addition to Lagrangian based theories we present the effective theory of DE and the  parametrization as general descriptions of cosmological gravity. Multimessenger GW detections can be used to measure the cosmological expansion (standard sirens), providing an independent test of the DE equation of state and measuring the Hubble parameter. Several key tests of gravity involve the cosmological propagation of GWs, including anomalous GW speed, massive graviton excitations, Lorentz violating dispersion relation, modified GW luminosity distance and additional polarizations, which may also induce GW oscillations. We summarize present constraints and their impact on DE models, including those arising from the binary neutron star merger GW170817. Upgrades of LIGOVirgo detectors to design sensitivity and the next generation facilities such as LISA or Einstein Telescope will significantly improve these constraints in the next two decades.
Contents
I Introduction
The Standard Model of Cosmology (or CDM) stands as a robust description of our universe. It is based on the theory of General Relativity (GR), which dictates the longrange gravitational interactions, together with the Cosmological Principle, which describes the geometry as homogeneous and isotropic on large scales. Standard matter (baryons, photons, neutrinos…) represents only a small fraction of the energy budget of the universe. The main ingredient is dark energy (DE), an unknown substance causing the late time acceleration. The other major component is dark matter (DM), an undetected constituent that seeds cosmic structures. The last piece of the SM of Cosmology are the initial conditions, which are thought to be set by an early period of quasiexponential expansion known as inflation. Despite the observational success of this model Aghanim et al. (2018), it remains as a puzzle the fundamental origin of each piece, which could be associated to new physics (see Fig. 1 for a summary of the different ingredients).
In the SM of Cosmology, the current accelerated expansion is explained by a constant energy density acting as a perfect fluid with negative pressure. Such a cosmological constant (CC) term is perfectly consistent with present observations but notoriously disagrees with theoretical expectations for the vacuum energy Weinberg (1989); Martin (2012). If this energy density is let to evolve in time, one naturally arrives to a dynamical description of DE sourced by a cosmological scalar field Copeland et al. (2006). If this field is now allowed to interact (nonminimally) with gravity, the possibilities to describe the cosmic expansion escalate Clifton et al. (2012). Alternatives to CDM offer the possibility to alleviate some of its tensions. For instance, DE models with an effective equation of state more negative than the cosmological constant could ease the tension between the local measurement of the Hubble constant and the inferred value from the cosmic microwave background (CMB). Exploring the largest scales with galaxy surveys like Euclid or LSST will help us understanding the expansion history of the universe and will provide new insights about gravity.
Gravity can be tested at different scales and regimes. Classical tests of gravity range from laboratory experiments to Solar System distances, and cover gravity in its weak field regime Will (2014). Astrophysical observations provide new avenues to improve these tests Berti et al. (2015). Pulsars in particular can be specially constraining, for instance with the recent observations of a triple stellar system Archibald et al. (2018). Tests in a much stronger regime has been performed tracking stellar orbits around the galactic center Hees et al. (2017). Altogether, these observations severely constrain modifications of GR. Theories beyond Einstein’s theory should thus resemble GR at small scales, e.g. hiding fifth forces with screening mechanisms Brax (2013); Joyce et al. (2015). At large scales, however, present constraints are considerably weaker. Combining different probes could be crucial to set an observational program to test gravity from cosmology Weinberg et al. (2013).
Gravitational wave (GW) astronomy offers the possibility to test gravity both in the strong regime and at large scales. So far there has been six individual detections, five binary blackholes (BBH) Abbott et al. (2016a, b, 2017a, 2017b, 2017c) and one binary neutron star (BNS) Abbott et al. (2017d). No GW background Abbott et al. (2018a), periodic source Abbott et al. (2018b) or longduration transient Abbott et al. (2018c) have been detected.
GWs could be critical in resolving the open problems of the SM of Cosmology. For instance, (non) observations of cosmological backgrounds of primordial GWs test inflation. BBHs events teach us about the population of BHs, which constrains their possible contribution to DM and their possible primordial origin Sasaki et al. (2018). Moreover, if DM is described by ultralight bosons or axions, it could resonate with pulsars Blas et al. (2017) or form clouds around BHs observable with GWs Arvanitaki et al. (2017). Finally, BNS with an associated counterpart such as GW170817 Abbott et al. (2017e, f) become standard sirens Abbott et al. (2017g) and allow to probe DE. In this review we will focus on this last case, exploring the possibilities of multimessenger GW astronomy to probe the nature of DE and the fundamental properties of gravity (see a schematic timeline of present and future facilities in Fig. 2).
i.1 Summary for the busy reader
Dark energy is the major component of the universe and yet its nature escapes our present understanding. Beyond the cosmological constant paradigm, a plethora of alternative theories of gravity has been proposed to explain the current cosmic acceleration (see Fig. 3 for a roadmap of possible modifications of gravity). We present an overview of the landscape of theories in Sec. II, as well as a summary of the different approaches to cosmological gravity (see Fig. 4 for a schematic diagram).
Gravitational wave astronomy opens new possibilities to probe gravity and DE. For readers unfamiliar with the basics of GWs, we provide a short introduction in Sec. III. For the purpose of cosmology, the most promising GW events are those that can be observed by other messengers (either EM waves or neutrinos). There are four main tests one can do with multimessenger GW events:

Standard sirens (Sec. IV): the amplitude of GWs is inversely proportional to its luminosity distance. If a counterpart of the GW is observed, a redshift measurement of the source is possible and the cosmic expansion history can be constrained. For close by sources, only the Hubble constant is measured. Future standard sirens measurements could help resolving the present tension in (see Fig. 8).

GW speed (Sec. V): the propagation speed of GWs follows from the dispersion relation. Once the location of a GW event is known, it is possible to compare the speed of GWs with respect to the speed of light. Many alternative gravity theories predict that GWs propagate at a different speed either by modifying the effective metric in which GWs propagate, by inducing a mass for the graviton or by introducing higher order terms in the dispersion relation.

GW damping (Sec. VI): modified gravity interactions can also alter the amplitude of GWs. In addition to the cosmic expansion, effective friction terms can damp GWs. This introduces an inequality between the GW and the EM luminosity distance that can be tested.

Additional polarizations (Sec. VII): in alternative theories of gravity, there could be additional modes propagating. These extra polarizations could be directly tested if the source is localized and there is a network of detectors online. Moreover, these modes could mix with the tensor perturbations leading, for instance, to GW oscillations.
In this review we aim at summarizing current bounds on gravity theories and dark energy models from the first multimessenger GW detection, GW170817. Up to date, the most constraining test is the GW speed. We also survey the prospects of different multimessenger tests with future detectors. Significant improvements can be achieved in probing the GW Hubble diagram with an increasing number of events. A schematic timeline of multimessenger GW astronomy is presented in Fig. 2 (the reader should be aware that expectations far in the future are very preliminary). The theoretical implications of present and future observations are discussed in Sec. VIII. We close the work in Sec. IX with an outlook of prospects and challenges of multimessenger GW tests of gravity and DE.
Ii Theories of gravity and Dark Energy
The quest to test gravity and find alternatives to the cosmological constant has produced many theories beyond Einstein’s General Relativity (GR) and other descriptions of gravity on cosmological scales. We will classify the different means to modify Einstein’s theory and review their status as descriptions of cosmic acceleration. Then we will review other general approaches to describe gravity on cosmological scales, namely through the effective theory of dark energy and phenomenological parameterizations of the gravitational potentials. The landscape of alternative theories is summarized in Fig. 3 and the approaches to cosmological gravity are schematically described in Fig. 4.
ii.1 Theories of Gravity
The actionbased approach to modify gravity is based on generalizing the EinsteinHilbert action
(1) 
where is Newton’s constant and denotes the action of matter, universally and minimally coupled to the metric . Variation of the action (1) with respect to the metric leads to Einstein’s field equations
(2) 
where is the Ricci tensor, the Ricci scalar and is the matter energymomentum tensor. Einstein’s equations can be used to obtain solutions for the spacetime () given the matter content () in any physical situation, including cosmological solutions relevant to study dark energy.
The structure of gravitational theories is severely restricted and several results can be used to prove the uniqueness of General Relativity under quite broad assumptions. Weinberg’s theorems restrict possible infrared (low energy) interactions of massless, Lorentz invariant particles, which for spin2 lead unavoidably to the equivalence principle Weinberg (1964) and the derivation of Einstein’s equations Weinberg (1965).^{1}^{1}1In addition to GR, there is another theory for massless, spin2 fields in 4D, Unimodular Gravity, which is invariant under diffeomorphisms preserving the 4D volume element van der Bij et al. (1982). At the classical level, the results of Lovelock imply that the EinsteinHilbert action is unique in 4D Lovelock (1971, 1972).
According to the above results, alternative theories of gravity can be classified into those that

Break the fundamental assumptions.

Include additional fields.

Make the graviton massive.
Note that those descriptions are not exclusive, and many theories fall within several categories. For instance: bimetric gravity has an additional field (tensor) and contains a massive graviton, EinsteinAether is both Lorentzviolating and includes a vector field, TeVeS has a scalar in addition to a vector, and many extradimensional models can be described in terms of additional fields in certain limits. Also, when referring to massive gravitons, we will be considering only classical spin2 fields.
ii.1.1 Breaking fundamental assumptions
The theorems that fix the structure of General Relativity assume a four dimensional pseudoRiemannian manifold and local interactions satisfying Lorentz invariance. Any departure from these principles offers a way to construct modified theories of gravity.^{2}^{2}2A class of GR extensions include additional geometric elements like torsion or nonmetricity. These elements can be viewed as either breaking the fundamental assumptions or including additional fields.
Extra dimensions:
Additional spatial dimensions allow the inclusion of new operators constructed only from the metric tensor. The canonical example are Lovelock invariants Lovelock (1971), such as the GaussBonnet term (a topological term in 4 dimensions which does not contribute to the equations of motion). The lack of observation of extra dimensions requires some mechanism to hide them. One example is compactification, when extra dimensions are sufficiently small that they are not accessible to experimental tests Bailin and Love (1987); Overduin and Wesson (1997). A radically opposite view consist on Braneworld constructions, in which the standard model fields live in a 3+1 dimensional brane, embedded in the higher dimensional space ArkaniHamed et al. (1998); Antoniadis et al. (1998); Randall and Sundrum (1999). The DvaliGabadadzePorrati (DGP) model Dvali et al. (2000); Nicolis and Rattazzi (2004) is one such construction in which selfaccelerating solutions^{3}^{3}3Selfaccelerating solutions are those in which there is a late time acceleration without a cosmological constant (). can be obtained. However, this branch of solutions is plagued by a ghost instability. The 4D effective theory can avoid this problem and was the origin of Galileon gravity Nicolis et al. (2009).
Lorentz Invariance Violation
Gravity can be extended by breaking Lorentz invariance. In many of these alternatives a preferred time direction emerges spontaneously breaking Lorentz symmetry (see Blas and Lim (2015) for a review). Hořava gravity Horava (2009) implements Lorentz violation through a preferred foliation of spacetime, with the attractive property that Lorentz symmetry can be recovered at low energies (see Blas et al. (2010, 2009); Sotiriou et al. (2009a, b); Sotiriou (2011) for extensions/variants) and leading to a powercounting renormalizable theory of gravity. Another class of Lorentzviolating theories is EinsteinAether, in which a vector field with constant norm introduces a preferred direction Jacobson (2007). The special case of EinsteinAether theories in which the vector field is the gradient of a scalar is known as Khronometric Blas et al. (2011a). Khronometric theories describe the lowenergy limit of some extension of Hořavagravity, linking the two frameworks Jacobson (2010). These ideas have been studied as cosmological scenarios Audren et al. (2013, 2015).
Nonlocal theories
Nonlocal theories include inverse powers of the Laplacian operator. These models can involve general functions (e.g. ) Deser and Woodard (2007); Koivisto (2008) or be linear (e.g. ) Jaccard et al. (2013). The latter class of models lead to phantom dark energy Maggiore (2014); Maggiore and Mancarella (2014) and are compatible with cosmological observations Dirian et al. (2015) (see Maggiore and Mancarella (2014) for a review). However, their viability on the solar system is disputed due to the time evolution of the effective degrees of freedom and the lack of a screening mechanism Barreira et al. (2014a). Nonlocal interactions have been also proposed as a means to improve the ultraviolet behavior of gravity Biswas et al. (2012); Modesto (2012). Nonlocal models are constructed using the Ricci scalar, since nonlocal terms involving contractions of the Ricci tensor give rise to cosmological instabilities Ferreira and Maroto (2013); Nersisyan et al. (2017).
ii.1.2 Additional fields
Gravity can be extended by the inclusion of additional fields that interact directly with the metric. These theories will vary by the type of field (scalar, vector, tensor) and the interaction with gravity it has. Theories with additional tensors (bigravity and multigravity) are extensions of massive gravity and will be described in Sec. II.1.3. We will assume a minimal universal coupling of matter to the metric. For a very complete review of gravity theories containing additional fields, see Ref. Heisenberg (2018a).
Scalar field
A scalar is the simplest field by which gravity can be extended. Scalars do not have a preferred orientation and thus a macroscopic, classical state can exist in the universe without affecting the isotropy of the spacetime if it depends only on time. Moreover, a potential term can mimic a cosmological constant very closely in the limit in which the field is varying very slowly (e.g. if the potential is very flat), which is the foundation of the simplest singlefield inflation and dark energy models (quintessence). Scalar fields may also arise in effective descriptions of fundamental theories belonging to other categories, such as braneworld constructions de Rham and Tolley (2010); Goon et al. (2011); Koivisto et al. (2014). These properties had lead to a proliferation of scalarbased models to describe accelerating cosmologies, both in the context of inflation and dark energy.
Recent efforts to study scalartensor theories have lead to a classification based on the highestorder derivatives of the additional field present in the action and the equations of motion, with three generations of theories

Oldschool scalar tensor theories: 1 order derivatives in the action, order in equations.

Horndeski theories Horndeski (1974): 2 order derivatives in the action and order in equations.

Beyond Horndeski: 2 order derivatives in the action and higher order in equations.
The classification is motivated by Ostrogradski’s theorem, which states that theories with second and higher (time) derivatives in the action generically introduce unstable degrees of freedom Ostrogradski (1850); Woodard (2015). While most physical theories belong to the first class, known loopholes to Ostrogradski’s theorem exits, for instance in effective or nonlocal theories (in which the ghost degrees of freedom are removed) Simon (1990) or when the theory is degenerate (that is, the inversion to canonical variables is not possible). The degeneracy condition is automatically satisfied if the equations of motion are second order, but that is not strictly necessary (different conditions appear when there are additional degrees of freedom Motohashi et al. (2016)). ^{4}^{4}4Scalartensor theories can be reformulated in terms of differential forms in which the second order equations follow naturally from the antisymmetry of this language Ezquiaga et al. (2016). This approach can be generalized to gravity theories with additional vector and tensor fields as well Ezquiaga et al. (2017). Known viable beyond Horndeski theories are known as Degenerate Higher Order Scalar Tensor (DHOST) Langlois and Noui (2016), which have second derivatives in the action (higher derivatives in the equations), but recently toy models with higher derivatives in the action have been proposed Motohashi et al. (2018).
Oldschool scalartensor theories contain at most first derivatives of the scalar in the action. They can be seen as a generalization of the JordanBransDicke theory of gravity Brans and Dicke (1961)
(3) 
where is the canonical kinetic term of the scalar field. This theory includes GR (), quintessence () Wetterich (1988); Ratra and Peebles (1988), BransDicke models Brans and Dicke (1961) (, ), kessence ArmendarizPicon et al. (1999, 2001) (, ). Archetypal modifiedgravity models such as Carroll et al. (2004); Hu and Sawicki (2007); Sotiriou and Faraoni (2010) are equivalent to instances of these theories De Felice and Tsujikawa (2010a). Chameleons Khoury and Weltman (2004) and symmetrons Hinterbichler and Khoury (2010) also belong to this class of theories (see Burrage and Sakstein (2016) for a review). Certain freedom exists in writing the theory due to the possibility of rescaling the metric and redefining the scalar field, i.e. the Jordan frame in which the metric is minimally coupled (3) and the Einstein frame in which is constant but matter is explicitly coupled to the scalar Flanagan (2004). Current cosmological observations constrain the BransDicke parameter (99%) Avilez and Skordis (2014).
Horndeski’s theory contains the best understood examples of scalartensor theories. The Horndeski action encompasses all local, 4D Lorentz invariant actions whose metric and field variation leads to second order equations of motion Horndeski (1974) (Horndeski’s theory is also known in the literature as Generalized Galileons Deffayet et al. (2011); Kobayashi et al. (2011)). Horndeski’s action reads
(4) 
where we have assumed minimal and universal coupling to matter in . The sum is over the four Lagrangians
(5)  
(6)  
(7)  
(8)  
where and are functions of and , and the subscripts and denote partial derivatives. Horndeski theories include all the generalized JordanBransDicke type, plus new additions that involve second derivatives of the scalar at the level of the action. These include kinetic gravity braiding (KGB) () Deffayet et al. (2010); Kobayashi et al. (2010); Pujolas et al. (2011), covariant galileons (, ) Nicolis et al. (2009); Deffayet et al. (2009), disformal Koivisto et al. (2012) and DiracBornInfeld gravity () de Rham and Tolley (2010); Zumalacarregui et al. (2013), GaussBonnet couplings Ezquiaga et al. (2016) and models selftuning the cosmological constant Charmousis et al. (2012a); MartinMoruno et al. (2015). Just as BransDicke is invariant under rescalings of the metric, Horndeski theories are invariant under fielddependent disformal transformations , which amount to a redefinition of the Horndeski functions (and the introduction of an explicit coupling to matter) Bettoni and Liberati (2013).
Theories beyond Horndeski have higher order equations of motion without including additional degrees of freedom. The first examples of these theories Zumalacárregui and GarcíaBellido (2014) were related to GR by a metric redefinition involving derivatives of the scalar field Bekenstein (1993),
(9) 
applied to the gravity sector. The simplest such beyond Horndeski theory emerged from the metric rescaling with derivative dependence , and was dubbed kinetic conformal gravity Zumalacárregui and GarcíaBellido (2014)
(10) 
where is an additional scalar field Lagrangian. One of the premises in constructing this type of theories was the existence of an inverse for the relation (9), which can be studied through the Jacobian of the mapping Zumalacárregui and GarcíaBellido (2014). If this assumption is broken the resulting theory is mimetic gravity Chamseddine and Mukhanov (2013), a gravitational alternative to dark matter. Interestingly, the conformal relation between kinetic conformal gravity (10) and GR ensures that this is one of the theories in which the speed of GWs is nontrivially equivalent to the speed of light Ezquiaga and Zumalacárregui (2017); Creminelli and Vernizzi (2017).
The best known beyond Horndeski theory is given by the GleyzesLangloisPiazzaVernizzi (GLPV) action Gleyzes et al. (2015a), which consists of Horndeski plus the additional Lagrangian terms:
(11)  
(12) 
Horndeski and GLPV Lagrangians of the same order, i.e. (7+11) or (8+12), can be mapped to Horndeski via showing the viability of these combinations Gleyzes et al. (2015a, b). For generic combinations of Horndeski and GLPV, viability arguments were first based on a special gauge (unitary gauge) that assumed that the scalar field derivative is timelike. Subsequent analyses eventually lead to covariant techniques to study the degeneracy conditions Langlois and Noui (2016) (see Deffayet et al. (2015) for earlier criticism). These techniques later showed that not all Horndeski and GLPV combinations met the degeneracy condition on a covariant level Crisostomi et al. (2016a).
The study of degeneracy conditions for scalartensor theories ultimately lead to the degenerate higherorder scalartensor (DHOST) Langlois and Noui (2016) paradigm classification of theories with the right number of degrees of freedom (also known as Extended ScalarTensor or EST) Crisostomi et al. (2016b). DHOST theories include cases beyond conformal kinetic gravity (10) and GLPV theories (11,12). DHOST theories are invariant under general disformal transformations (9), which can in turn be used to classify them Ben Achour et al. (2016a) (see also Crisostomi et al. (2017)). DHOST theories have been fully identified including terms with up to cubic secondfield derivatives in the action, e.g. Ben Achour et al. (2016b). Demanding the existence of a Poissonlike equation for the gravitational potential restricts the space of DHOST theories to those that are related to Horndeski via disformal transformations (9) Langlois et al. (2017).
Vector field
Theories with vector fields have been proposed as modifications to GR and in the context of dark energy. A background vector field does not satisfy the isotropy requirements of the cosmological background, unless it points in the time direction and only depends on time . Isotropy can also happen on average, if a vector with a spacelike projection oscillates much faster than the Hubble time Cembranos et al. (2012). In that case the background is isotropic on average but the perturbations (including gravitational waves) inherit a residual anisotropy Cembranos et al. (2017). Finally, theories with multiple vectors can satisfy isotropy, for instance, if they are in a triad configuration ArmendarizPicon (2004).^{5}^{5}5Technically speaking, multiple vectors can lead to isotropic solutions if they have an internal symmetry that together with the broken spacetime symmetries leaves a residual Beltran Jimenez and Heisenberg (2018). For the case of the triad, the symmetry group is . A large number of vectors can also lead to statistical isotropy (e.g. if the orientations are random) Golovnev et al. (2008). The kinetic term for a vector field, , is defined by the gauge invariant field strength and the addition of a mass term is known as Proca theory Proca (1936).
Proca theories have been generalized to include explicit gravitational interactions of a massive vector field Tasinato (2014); Heisenberg (2014); Allys et al. (2016a); Beltran Jimenez and Heisenberg (2016). The vector field Lagrangian is built so that precisely one extra (longitudinal) scalar mode propagates in addition to the two usual Maxwelllike transverse polarizations. Its full generalization contains terms with direct couplings between the vector and spacetime curvature, whose structure closely resembles those of Horndeski’s theory (7,8). In analogy to beyond Horndeski, there are also beyond generalized Proca interactions Heisenberg et al. (2016); Kimura et al. (2017). Further extensions to multiple vector fields known as generalized multiProca/YangMills theories are able to incorporate new couplings Allys et al. (2016b) and configurations Beltran Jimenez and Heisenberg (2017), e.g. the extended triad , as do theories with a vector and a scalar (ScalarVectorTensor) Heisenberg (2018b). For more details about these theories we recommend Ref. Heisenberg (2018a).
An iconic theory containing a vector is the TensorVectorScalar (TeVeS) theory by Bekenstein Bekenstein (2004). TeVeS emerged as a relativistic theory able to describe Modified Newtonian Dynamics (MOND), and thus as an alternative to dark matter. For an overview of fieldtheoretical aspects of TeVeS and related theories, including other relativistic MOND candidates, see Ref. Bruneton and EspositoFarese (2007). TeVeS theory introduces several nonminimal ingredients. In addition to the gravitational metric matter is minimally coupled to an effective metric
(13) 
which generalizes the scalar disformal relation (9), incorporating the vector. Here is the gravitational metric, is the scalar. The vector is enforced to be timelike and normalized with respect to the gravitational metric . TeVeS has a very rich phenomenology, including effects in GW propagation Sagi (2010). At the level of cosmology it is partially able to mimic DM, although the oscillations of the fields make it hard for the theory to reproduce the peaks in the CMB Skordis et al. (2006); Bourliot et al. (2007); Skordis (2009).
ii.1.3 Massive gravity and tensor fields
Giving a mass to the graviton is another means to extend GR, with gravity mediated by a particle with mass , spin and polarization states (see de Rham et al. (2017) for bounds on the graviton mass). Weinberg theorem on the structure of GR relies on the infrared properties of the interactions: a mass term changes this structure. Despite this clear loophole, constructing a selfconsistent theory of massive gravity, free of pathologies and with the right number of degrees of freedom proved an extremely hard endeavor that took nearly 70 years to complete. The linear theory of massive gravity was formulated in 1939 by Fierz & Pauli Fierz and Pauli (1939) as linearized GR plus a mass term
(14) 
It was later found that FierzPauli theory was discontinuous and gave different results from GR in the limit (vDVZ discontinuity) van Dam and Veltman (1970); Zakharov (1970). The discrepancy is due to the longitudinal polarization of the graviton (the helicityzero mode) not decoupling in that limit. Considering nonlinear interactions solved the apparent discontinuity by hiding the helicityzero mode, which is strongly coupled in regions surrounding massive bodies and effectively decouples, recovering the GR predictions when Vainshtein (1972). Despite this progress, massive gravity had another important flaw: all theories seemed to have an additional mode (known as BoulwareDeser (BD) ghost) that renders the theory unstable Boulware and Deser (1972); Creminelli et al. (2005).
Ghost Free Massive Gravity
The apparent difficulties were overcome in de RhamGabadadzeTolley theory (dRGT) de Rham et al. (2011), also known as ghostfree massive gravity (for current reviews on the theory see Hinterbichler (2012); de Rham (2014)). dRGT is a ghost free theory propagating the 5 polarizations corresponding to a spin2 massive particle, universally coupled to the energymomentum tensor of matter (cf. Fig. 6). The ghostfree property was initially shown in the decoupling limit (in which the helicity0 mode decouples from the other polarizations) and then in the full theory Hassan and Rosen (2012a, b) . The phenomenological deviations induced by massive gravity are primarily due to the helicity0 mode. On small enough scales the Vainshtein mechanism Vainshtein (1972) (see Babichev and Deffayet (2013) for a review) effectively suppresses these interactions, leading to predictions very similar to GR on Solar System scales (however, new classes of solutions for black holes do exist, in addition to the usual ones Babichev and Brito (2015)).
Massive gravity may offer a solution to the accelerating universe. A heuristic argument is that the force mediated by the massive graviton has a finite range , weakening over distances larger than the Compton wavelength of the graviton . Hence, if the mass of the graviton is then gravity weakens at late times and on cosmological scales, causing an acceleration of the cosmic expansion relative to the GR prediction. The program to apply massive gravity as a dark energy model has hit important barriers, as flat FLRW solutions do not exist in this theory D’Amico et al. (2011). Accelerating solutions without a cosmological constant (CC) do exist with open spatial hypersurfaces Gumrukcuoglu et al. (2011), but they are unstable De Felice et al. (2013). Proposed solutions include the graviton mass being generated by the vacuum expectation value of a scalar D’Amico et al. (2011) or deformations of the theory in which the BD ghost is introduced, which provides dynamical and accelerating but metastable solutions Könnig et al. (2016). Other ways to make massive gravity dynamical include the addition of a new field, such as a scalar field, e.g. quasidilaton D’Amico et al. (2013), or one (or several) dynamical tensors in bigravity (and multigravity).
Bigravity and Multigravity
In order to write a mass term for the metric, dRGT incorporates an additional, nondynamical tensor, akin to the occurrence of in eq. (14). Massive gravity can be extended by including a kinetic term to the auxiliary metric, which becomes fully dynamical. This leads to the theory of bigravity (or bimetric gravity) Hassan and Rosen (2012c), which contains two spin2 particles: one massive and one massless. The same procedure can be extended to more than two interacting metrics, leading to multigravity theories Hinterbichler and Rosen (2012). In these constructions there is always one massless excitation of the metric (a combination of the different tensor fields), with all other excitations being massive.
Bigravity solves the problem of cosmological evolution, at least at the background level. Flat FLRW solutions do exist, and many viable expansion histories have been found to be compatible with data Akrami et al. (2013). However, it was later found that these models had instabilities that affected the growth of linear perturbations Comelli et al. (2014), which were found to be quite generic across different branches of solutions Könnig (2015). In some cases the instabilities affect only scales sufficiently small for nonlinear effects to be important (i.e. the Vainshtein mechanism) which might render the theory stable Mortsell and Enander (2015). Another solution is to choose the parameters of the theory so instabilities occur at early times, when characteristic energies are high and bigravity is not a valid effective field theory. This happens by choosing a large hierarchy between the two Planck masses: the soobtained theory is practically indistinguishable from GR plus a (technically natural) CC Akrami et al. (2015).
ii.2 Descriptions of cosmological gravity
The immense variety of alternative theories has motivated the search for effective descriptions able to capture the phenomenology of generic dark energy models. The covariant actions approach reviewed in Sec. II.1 offers several advantages, including 1) full predictivity, as (classical) solutions can be found from microscopic scales, to strong gravity and all the way to cosmology, 2) selfconsistency, as different regimes can be computed for the same theory, leading to tighter constraints when the data is combined. For instance, following this approach, we discuss the cosmology of covariant Galileons in Sec. IV.1. Nonetheless, a great downside of this approach is that the predictions for every model/theory have to be obtained from scratch, which makes the exploration of the theory space a daunting task.
An alternative route is to constrain deviations from GR, without reference to any fundamental theory. The tradeoff is to keep the theory of gravity as general as possible at the expense of dealing with a very simple spacetime. The simplest situation is where the background spacetime is flat and maximally symmetric (Minkowski), a setup useful to model gravity in the Solar System. In this simple case one can define a series of quantities, known as Parameterized PostNewtonian (PPN) coefficients, that describe general modifications of gravity over Minkowski space (see Ref. Will (2014) for details, including constraints and additional assumptions). These PPN parameters that can be constrained by experiments (such as the deflection of light by massive bodies) and computed for any theory, and thus provide a very efficient phenomenological dictionary.
In cosmology we are interested in describing gravity over a slightly less symmetric background: a spatially homogeneous and isotropic, but time evolving, FriedmannLemaitreRobertsonWalker (FLRW) metric:
(15) 
where metric perturbations are in Newtonian gauge with the sign conventions of Ma & Bertschinger Ma and Bertschinger (1995). The tensor perturbation is symmetric, transverse and traceless () and we have ignored vector perturbations. The timeevolution of the cosmological background makes an extension of PPN approach to cosmology a difficult task, as instead of constant coefficients one needs to deal with functions of time due to the evolution of the universe.
The most important example of an effective description in cosmology is the parameterization of the cosmological background, often done in terms of the equation of state Chevallier and Polarski (2001); Linder (2003). Instead of computing the modifications to the Friedmann equations and the pressure and energy density contributed by the additional fields, a general approach to cosmological expansion is to specify so that
(16)  
(17) 
This is sufficient to describe any cosmological expansion history and in any theory (as long as matter is minimally coupled) just by using the Friedmann equation (16) as a definition for .
Describing the perturbations requires more functional freedom. Here we will review two common procedures, namely the effective theory of dark energy and the modified gravitational “constants”. The different approaches (including the covariant theory approach), their features and connections are outlined in Fig. 4. Consistency checks between the background and perturbations can also be used to test the underlying gravity theory Ruiz and Huterer (2015); Bernal et al. (2016a).
ii.2.1 Effective theory of Dark Energy
The effective (field) theory of dark energy (EFTDE) Gubitosi et al. (2013); Bloomfield et al. (2013); Gleyzes et al. (2013) can be used to systematically describe general theories of gravity over a cosmological background (see Ref. Gleyzes et al. (2015c) for a review). The original formulation applies to theories with a scalar field and uses the unitary “gauge”: a redefinition of the time coordinate as the constant hypersurfaces (this is always possible if is timelike and nondegenerate, as in perturbed cosmological backgrounds, but not in general). One then constructs all the operators compatible with the symmetries of the background (recalling that the time translation invariance is broken by the cosmological evolution).
A very convenient basis for the EFT functions was proposed by Bellini & Sawicki Bellini and Sawicki (2014), when restricted to Horndeski’s theory. In their approach the EFT functions are defined by the kinetic term of the propagating degrees of freedom in the equations of motion. The dynamical equation for tensor perturbations
(18) 
introduces two dimensionless functions

tensor speed excess describes the modification in the GW propagation speed . This modification is frequency independent (see Sec. V).

Planckmass run rate enters as a friction term. It is related to the cosmological strength of gravity (the kinetic term of tensor perturbations) by (see Sec. VI).
The equations in the scalar sector (eqs. (3.20), (3.21) of Bellini and Sawicki (2014)) can be used to define the remaining functions. If we look only at the second time derivatives (that is, the kinetic terms)
(19)  
(20) 
(note the ellipsis denote terms without second time derivatives) one can define

kineticity modulates the “stiffness” of the scalar field (how hard it is to excite perturbations in ). The kineticity is intimately related to the propagation speed of scalar perturbations, which satisfies : higher kineticity values lead to slower scalar waves and vice versa.
These functions can be computed from the Lagrangian functions in (4), and for a given theory will depend on the value of the scalar field and its time derivative. Constraints on the functions can also be used to reconstruct the terms in a fundamental theory, as shown in Tab. 1. Systematic reconstructions of the Lagrangian from the functions have been also explored Kennedy et al. (2017, 2018).
Horndeski  DHOST  
GLPV  
zero, nonzero (arbitrary), nonzero (constrained)
Increasingly complex theories of gravity lead to a larger number of EFT functions. In beyond Horndeski theories of the GLPV type, e.g. (11,12), a new function is introduced Gleyzes et al. (2015b) which phenomenologically produces a weakening on gravity on small but linear cosmological scales D’Amico et al. (2017). In DHOST theories including (10) the situation is more involved, as the new EFT functions () need to be related to each other and by the degeneracy conditions that prevent the introduction of additional degrees of freedom Langlois et al. (2017): this leads to two classes of theories with one free function, which is either or one among . New EFT functions appear beyond scalartensor theories, as has been explicitly derived for vectortensor Lagos et al. (2016) and bimetric Lagos and Ferreira (2017) theories (including bimetric gravity), with a unified treatment of theories with different degrees of freedom Lagos et al. (2018).
Different versions of the linear EFTDE approach has been implemented in numerical codes able to obtain predictions based on linear perturbation theory. Publicly available implementations exist in EFTCAMB Hu et al. (2014), hi_class Zumalacárregui et al. (2017) and COOP Huang (2016), with the first two based on the CAMB and CLASS Boltzmann codes Lewis et al. (2000); Blas et al. (2011b). In addition, the CLASSGal code (integrated into CLASS) can be used to compute relativistic corrections to cosmological observables Di Dio et al. (2013). These and other codes have been tested against a large class of models at a level of precision sufficient for current and nextgeneration cosmological experiments Bellini et al. (2018).
The EFT framework has been tested using linear observables. Horndeski theories were tested against current experiments, leading to constraints on the functions varying over Bellini et al. (2016), with Ade et al. (2016) and setting to reflect the strong bounds on the GW speed Kreisch and Komatsu (2017) ( is very weakly constrained by current data). Future experiments have great potential to improve on these bounds, and are expected to improve the sensitivity to Gleyzes et al. (2016); Alonso et al. (2017); Lorenz et al. (2018); A. et al. (2018); Reischke et al. (2018). EFTbased modifications of gravity might be observable through relativistic effects on ultralarge scales Renk et al. (2016); Lorenz et al. (2018) (see also the discussion in Sec. II.2.2): these techniques might improve significantly our ability to constrain , although it will remain the hardest to measure Alonso et al. (2017). Those works used simple functional dependence of the EFT functions. It has been nonetheless shown that simple parameterizations are indistinguishable from more complex models in most cases, even for nextgeneration cosmology experiments Gleyzes (2017).
The EFT approach has been generalized beyond linear perturbations for Horndeski theories. Including nonlinear cosmological perturbations in general introduces new functions at every order in perturbation theory (e.g. to compute the bispectrum Bellini et al. (2015)). However, a restriction to cubic and quartic operators (in the unitary gauge) leads to only 3 new operators on quasistatic scales Cusin et al. (2018a). Some applications of nonlinear EFTDE include corrections to the power spectrum (e.g. Cusin et al. (2018b)), the use of higherorder correlations as a test of gravity, such as the bispectrum of mattter Bellini et al. (2015), galaxies Yamauchi et al. (2017) and CMB lensing Namikawa et al. (2018) or the the nonlinear shift of the BAO scale Bellini and Zumalacarregui (2015).
ii.2.2 Modified Gravitational “constants”
A very commonly used approach employs general modifications of the equations relating the gravitational potentials to the matter density contrast
(21)  
(22) 
(note that different conventions exist in the literature). Here is the density contrast in Newtonian gauge (15) and the functions parameterize the evolution of the gravitational potentials as a function of time and scale . The functions are often referred to as , because gradients of determines the force felt by nonrelativistic particles and those of the geodesics of massless particles (and thus the lensing potential). The ratio of the gravitational potentials,
(23) 
is of particular interest, since GR predicts that it is exactly one in the absence of radiation and any sizable deviation could be an indication of modified gravity.
This approach has numerous advantages as a test of gravity against data. It is completely theory agnostic, not requiring any information on the ingredients or laws of the theories being tested. Most importantly, it is completely general for universally coupled theories: given any solution it is possible to obtain through (21,22). In this sense, any finding of might point towards deviations from GR and warrant further investigation.
The main shortcoming of this approach is its great generality: any practical attempt to implement (21,22) requires a discretization of the functional space, introducing free parameters for a homogeneous binning. In contrast, the EFT approach for Horndeski theories (18,19) requires only parameters, making it a more economic parameterization for all but the simplest scaledependencies (). Capturing the full scale dependence of requires either a large parameter space or assumptions about the dependence.
A common practice to overcome this limitation is to choose a functional form for as a function of scale. For Horndeski theories the functional form is a ratio of quadratic polynomials in De Felice et al. (2011); Amendola et al. (2013)
(24) 
for functions that depend on redshift through the theory (4) and the scalar field evolution. The mapping is exact on small scales in which the field dynamics can be neglected, below scalar sound horizon Sawicki and Bellini (2015). A dependence as the ratio of polynomials is generic in local theories at quasistatic scales Silvestri et al. (2013), with higher order polynomials possible in Lorentzviolating Baker et al. (2014), multifield Vardanyan and Amendola (2015) theories. Studies with current data have tested rather simple parameterizations of : for instance the Planck survey tested the case of independent in addition to the theorymotivated (24) Ade et al. (2016). Future surveys will improve the resolution on the scaledependence: 3 bins are the minimum to constraint all the parameters in eq. (24), with 6 bins in Amendola et al. (2014a); Taddei et al. (2016). A limited handle on scaledependence on ultralarge scales might be achievable Baker and Bull (2015); Villa et al. (2018) (see also Raccanelli et al. (2014); Lombriser et al. (2013); Bonvin and Fleury (2018) for related parameterizations).
Another main shortcoming of the completely general approach is that there is no information from other regimes. The major setback with respect to EFT is the lack of information from gravitational wave observables, while in EFT the tensor and scalar sectors are modified accordingly i.e. GW data restrict the modifications available to scalar perturbations, for instance, theories with require either or to be nonzero Saltas et al. (2014). Attempts to explore the connections between and the EFT approach in Horndeskilike theories have used very general parameterizations: connecting theoretical viability conditions of the theory with the behavior of Perenon et al. (2015), including the case with to address the impact of the GW speed measurement Peirone et al. (2018a). General properties of Horndeski theories could be inferred from detailed measurements of Pogosian and Silvestri (2016). Similarly to the EFT approach, the background evolution is unknown and the equation of state (17) is in principle arbitrary. However, theoretical priors on can be obtained for broad classes of Lagrangians (e.g. quintessence Marsh et al. (2014)) or from stability conditions in general realizations of the EFT functions Raveri et al. (2017).
Iii Basics of Gravitational Waves
Gravity is a universal, longrange force. This, in field theory language, implies that it must be described by a metric field in order to manifestly preserve locality and Lorentz invariance. At low energies, the leading derivative interactions are second order. Therefore, gravity theories generically predict the existence of propagating perturbations or, in other words, the existence of GWs. One can define a metric perturbation as a small difference between the metric field and the background metric
(25) 
where . However, in curved space it is nontrivial to distinguish the perturbation from the background unless the latter posses some degree of symmetry, e.g. flat space or FLRW. A way out is to define GWs via geometric optics Misner et al. (1973). In this context, the key element to distinguish the GW from the background is the size of the fluctuations with respect to the typical size of the background variation . One could associate the typical variation scale in the background with the minimum value of the components of the background Riemann tensor
(26) 
For astrophysical sources, we will see later that the wavelength of the GW is orders of magnitude smaller than the typical variations of for cosmological setups. The fact that implies that there is a clear hierarchy between background and perturbations, allowing to solve the problem using an adiabatic (or WKB) expansion.
In the following, we describe the basics of GWs. We begin by introducing GWs in GR. Then, we explore the propagation in cosmological backgrounds. Subsequently, we describe how this picture is changed beyond GR. Finally, we discuss the status of present and future GW detectors. We recommend the reader Ref. Misner et al. (1973); Maggiore (2008, 2018); Flanagan and Hughes (2005); Carroll (2004) for more details.
iii.1 GWs in General Relativity
General Relativity is a universal, infiniterange force. As we have seen in the previous section, this implies that it is described by a massless, spin2 field. The dynamics is described by Einstein’s equations (2). Importantly, not all the components of the Einstein tensor contain second order time derivatives of the metric . This implies that not all of the 10 components of will propagate. In particular, the equations act as 4 constraint equations. This, together with the 4 unphysical modes reduced by the gauge choice, leaves only 2 propagating degrees of freedom. This is precisely what one would expect for a massless spin2 particle.
In order to study GWs, the next step is to study the linearized Einstein’s equations. To diagonalize the equations for the tensor perturbations, one has to introduce the tracereversed perturbation
(27) 
whose name comes from the fact that where and are the traces of and respectively. Fixing the Lorenz gauge for this new variable , the linearized Einstein equations in curved spacetime read
(28) 
where covariant derivatives are built with the background metric . Here, we have introduced the perturbed energymomentum tensor as the difference of the total energy momentum tensor with respect to the background solution . One should note that, in vacuum, all the Ricci tensors vanish in the second line. Moreover, for shortwave GWs , the Riemann tensor in the first line has a subdominant contribution.
To deal with the two GW polarizations, it is convenient to work in the transversetraceless (TT) gauge, which is defined by
(29) 
Note that in the TT gauge, the tracereversed perturbation is equal to the original perturbation . If the GW is propagating in the direction, the spatial components become
(30) 
with and being the two polarizations of GR.
iii.1.1 Generation
A first question to address is how GWs are produced. Let us consider a GW source in vacuum within the shortwave approximation. Then, the general propagation equation (28) reduces to
(31) 
This wave equation can be solved in analogy to electromagnetism using a Green’s function. In terms of the retarded time , the solution is
(32) 
For an isolated, far away, nonrelativistic source, this solution can be simplified. In fact, one can make a multipole expansion. The zeroth moment corresponds to the massenergy of the source . However, conservation of energy for an isolated source tells us that cannot vary in time. Next, the mass dipole moment is associated to the motion of the center of mass. Nevertheless, its time derivate is the momentum of the source that also has to be conserved^{6}^{6}6Similar arguments apply for the spin angular momentum in case the source exhibit some internal motion.. Consequently, the leading contribution is the mass quadrupole moment , which generates GWs through its second time derivatives
(33) 
For a binary system of masses and , the quadrupole radiation is
(34) 
where is a function of the orientation of the binary that depends on the polarization or (recall (30)), is the phase and we have introduced the chirp mass
(35) 
As the masses orbit one around the other, they will lose energy with the emission of GWs. They will begin getting closer and orbiting faster until they eventually merge. Thus, the frequency of GWs will increase with a characteristic chirp signal following
(36) 
Note that to consider the energy loss due to GWs emission one has to go to second order in perturbation theory. An example of the typical GW strain and frequency of a compact binary coalescence is presented in Fig. 5.
Typical binary compact objects emitting detectable GWs are binary neutron stars (BNS) and binary blackholes (BBH). The order of magnitude of the frequency of the GWs of these systems is
(37) 
where is equal to one solar mass. This implies that higher masses lead to lower frequencies. In terms of the wavelength one finds
(38) 
This allows us to compare the size of the wavelength with the typical size of the background curvature variation . For cosmology, the size of the curvature is related to the Hubble horizon m. For our galaxy one can estimate m and for the Solar System m. As it can be observed, the geometric optics expansion is an excellent approximation due to the great hierarchy between and . Only GWs passing near a very dense object such as a BH, km, would break this shortwave approximation.
The typical amplitude of a GW from a compact binary can be estimated using (34), leading to
(39) 
Contrary to EM waves, GW detectors are directly sensitive to the amplitude of the wave, which falls like and not as the luminosity . This means that even if the amplitudes are very small, GW detectors are more sensitive to distant sources.
iii.1.2 Propagation
Once the GW is generated, it will propagate in vacuum following
(40) 
A general solution of this wave equation can be written as the sum of plane waves
(41) 
where Re denotes the real part. By plugging this expression in the wave equation and expanding in powers of , one finds at leading order that
(42) 
Therefore, GWs propagate in null geodesics determined by the background metric. This means that the GWcone is the same as the lightcone and both waves propagate at the same speed. Moreover, the wave is transverse to the propagation direction
(43) 
similarly to electromagnetic waves. Finally, by defining the scalar amplitude one realizes that
(44) 
which can be interpreted as the conservation of gravitons. One should note that in the wave equation only modifies the amplitude at second order. Consequently, at first order in geometricoptics, the wave equation can be rewritten as
(45) 
This expression could be used as a gauge invariant, coordinate independent definition of the propagation of GWs in vacuum.
iii.1.3 Detection
To see the effect of a GW passing by, one has to study the deviation of nearby geodesics. Given two particles with fourvelocity separated by , their separation evolves as
(46) 
where is the proper time. At leading order, the four velocity is just the unit vector , and we only have to compute the Riemann tensor in the TT gauge. The result is
(47) 
where we have also used that at leading order the proper time and the coordinate time coincide. Accordingly, only the components of the separation vector transverse to the propagation vector will feel the effect of the GW. In these directions, the separation between the test particles will oscillate as the GW travels perpendicular to them. In Fig. 6, we plot the effect of the different GW polarizations crossing a circle of test masses.
GW detectors precisely rely on this principle that GWs can alter the separation between test masses. Modern detectors are interferometers. In brief, they are constituted by two perpendicular arms of the same length with two mirrors in free fall at their ends (acting as test particles). A laser beam is split in the two arms so that the beams reflect in each mirror and come back to the splitting point. In the absence of a GW, both laser beams returning will interfere destructively and no signal would arrive to the detector. However, if a GW crosses the interferometer, it will change the length of the arms differently. This means that the laser beams will take different times to travel the arms, arriving at the splitting point with different phases. Then, the destructive interference is lost and some signal gets to the detector.
Note that the typical distance variation of two test masses separated by is approximately . For compact binaries, we have seen that the strain amplitude is . Therefore, LIGOtype detector with arms of the order of kilometers have to measure distance variations
(48) 
a thousand times smaller than the nucleus of an atom. To achieve that, each arm has a resonant cavity in which the laser beams bounce back and forth about 300 times. This effectively makes groundbased interferometer arms to be km long (since the variation time of the GW is much longer the travel time of the laser in the cavity). Accordingly, LIGO is sensitive to frequencies of Hz. For the future spacebased interferometer LISA, the working principle will be the same but with longer arms km, being thus sensitive to much smaller frequencies, Hz.
iii.2 GWs in cosmology
At large scales, the universe is homogeneous and isotropic to very high accuracy. The background geometry is then described by a (flat) FriedmannLemaitreRobertsonWalker (FLRW) metric
(49) 
where is the scale factor and we are timing in conformal time . The propagation equation (40) becomes in Fourier space
(50) 
where is the Hubble parameter and primes denote derivatives with respect to conformal time. This is nothing but a wave equation with a friction term due to the cosmic expansion. This Hubble friction will produce a redshift of the frequencies and a rescaling of the GW amplitude . The previous formulae for a compact binary (3436) written in terms of the observed frequency are thus valid if we replace the chirp mass by the redshifted chirp mass
(51) 
and the physical distance by the GW luminosity distance
(52) 
where is the speed of light and the redshift. In this way, all the terms cancel each other. Note that there is an intrinsic degeneracy between the redshift and the Hubble parameter in the GW luminosity distance. Therefore, the expansion history can only be obtained from the GW amplitude if the redshift is known. For near by sources , the Hubble constant can be obtained
(53) 
showing the power of GW astronomy to do cosmology. We will review this topic in more detail in Sec. IV.
Finally, let us mention that we have only focused on GWs from binary sources in the late universe. However, there could be other sources of GWs in the early universe leading to stochastic, cosmological backgrounds. For a nice review in the subject one can follow Caprini and Figueroa (2018). One may wonder if there could be an effect in the GW propagation when traveling through the cold dark matter. This question has been addressed recently and the answer is that the effect is too small Baym et al. (2017); Flauger and Weinberg (2018).
iii.3 GWs beyond GR
As we have emphasized at the beginning of this section, the existence of wave solutions for metric perturbations is generic of second order gravity theories. However, the behavior of these GWs can be very different depending on the gravity theory. The differences can arise either at the production or the propagation.
iii.3.1 Additional polarizations
During the generation of GWs, the main differences in theories beyond GR is that there could be another polarizations excited. We have seen that in GR only the 2 tensor polarizations propagate (recall (30)). Nevertheless, modifications of gravity might introduce new degrees of freedom. For instance, in scalartensor theories there will be an additional scalar mode. Or in Massive Gravity, where there will be in addition 2 vector modes and a scalar one. For a GW propagating in the direction, one could decompose the amplitude in the different polarizations
(54) 
where and are the two tensor modes, the two vector polarizations, the transverse (breathing) scalar and the longitudinal scalar mode. One should note that these other types of polarizations will also leave an imprint in the detectors. Each polarization will have a different effect as we exemplify in Fig. 6. In principle, with a set of 6 detectors one could distinguish all possible polarizations.
Before continuing, it is important to remark that if a source is emitting additional polarizations, it will lose energy more rapidly. For binary pulsar, if additional modes were emitted, the orbit would shrink faster due to the higher energy loss. For PSR B1913+16 (better known as HulseTaylor pulsar) Hulse and Taylor (1975), the orbit has been tracked for more than four decades now, showing an impressive agreement with GR Weisberg et al. (2010). Binary pulsars have been intensively used to constrain alternative theories of gravity, placing severe bound on dipolar radiation as reviewed in Stairs (2003); Wex (2014). An example of this are EinsteinAether propagating waves Jacobson and Mattingly (2004), which have been constrained from pulsars due to dipolar GW emission Yagi et al. (2014a, b). Another would be the constraints on BransDicke from a pulsarwhite dwarf binary Freire et al. (2012).
Due to these constraints on the emission of additional polarizations, it is usually invoked a screening mechanism around the source to evade them. If this is the case, deviations of GR could only be measured in the propagation of GWs. We will discuss more about the emission of extra modes and screening mechanisms in Sec. VIII.2.
iii.3.2 Modified propagation
The propagation of GWs in gravity theories beyond GR can be very complicated. The additional fields might modify the background over which GWs propagate and their perturbations could even mix with the metric ones. For simplicity, we will restrict here to cosmological backgrounds. In that case, due to the symmetries of FLRW, tensor perturbations can only mix with other tensor perturbations. Possible deviations from the cosmological wave equation in GR (50) can be parametrized by Nishizawa (2017)
(55) 
where is an additional friction term, accounts for an anomalous propagation speed, is an effective mass and is a source term originated by the additional fields. For instance, the scalartensor analogue of this equation is (18). It is interesting that the modified GW propagation can also be understood in analogy with optics as GWs propagating in a diagravitational medium Cembranos et al. (2018).
Focusing in the case without sources, , the original GR waveform , given by (34) for instance, will be modified by
(56) 
where we have introduced . Mainly, the additional friction will modify the amplitude while the anomalous speed and the effective mass change the phase. The modified luminosity distance is then
(57) 
We will discuss how to test the GW phase in Sec. V and the damping of the strain in Sec. VI.
For GWs propagating in FLRW backgrounds, a source is present when there are additional tensor modes propagating. A paradigmatic example of this is bigravity, where there are two dynamical metrics. In that case, we have to track the evolution of both metric perturbations Narikawa et al. (2015); Max et al. (2017, 2018)
(58) 
where for shortness we have absorbed the Hubble friction in the definition of the perturbation and we do not show the spatial indices. Here is the effective mass (one of the tensor fields is massive) and is the mixing angle. Since there are interactions between and , this means that the mass eigenstates are not the same as the propagation eigenstates. In analogy with the propagation of neutrinos, there can be GW oscillations. In Sec. VII.1 we will see how GW oscillations can be tested. One should note that the possibility of having GW oscillations is not restricted to bigravity. Any gravity theory in which the additional degrees of freedom can arrange to form a tensor perturbation over FLRW background could display the same phenomenology. In particular, this is what happens with gauge fields in a SU(2) group Caldwell et al. (2016); Beltran Jimenez and Heisenberg (2018).
iii.4 Present and future GW detectors
Before presenting the different tests of gravity with multimessenger GW astronomy, let us outline briefly the status of present and future GW detectors. We summarize the different sensitivities of each detector and the typical sources in Fig. 7. The capabilities of multimessenger GW astronomy depend mainly on two aspects:

Number of detections: this is most sensitive to the size of the volume of the Universe covered by the GW detector. However, there is a large uncertainty in the actual population of the sources, e.g. BNS.

Sky localization: this is most sensitive to the number of detectors that allow for a better triangulation of the source. A better localization of the GW events simplifies the search for a counterpart.
We draft a summary of present expectations for the range of detection and localization angle of different GW detectors in Fig. 2. The reader should be aware that these expectations, specially the ones far in the future, might be subject to important modifications.
At present, we are in the second generation (2G) of groundbased detectors. There have been already two operation runs. In the first one, only the two aLIGO detectors were online with a detection range for BNS of the order of 80 Mpc. In the second one, aVirgo joined. Although its sensitivity was still lower, aVirgo helped to reduce the localization area an order of magnitude, from to . For illustration, we plot in Fig. 7 the strain of the first event GW150914 Vallisneri et al. (2015).
However, neither aLIGO nor aVirgo has reached their designed sensitivity yet. Moreover, other two 2G detectors are on the way. KAGRA Somiya (2012) in Japan is under construction and it is expected to start operating in 2020. On the other hand IndIGO et al. (2011), a replica of LIGO located in India has been approved. This means that in the coming years two main improvements are expected: a larger event rate and a more precise localization Abbott et al. (2018d). The range of detection is expected to improve by a factor of 3 implying a factor 27 in the detection rate. The localization is expected to reduce to areas of with KAGRA and to a few with IndIGO. Note that this is a key point in order to associate any counterpart with a GW event.
A third generation (3G) of groundbased detectors is being planned. The European 3G proposal is the Einstein telescope (ET) Sathyaprakash et al. (2012), an underground, three 10kmarms detector. Its current design aims at improving by a factor of 10 present sensitivity. The US 3G proposal, Cosmic Explorer (CE) Abbott et al. (2017h), is more ambitious with two 40km arms further improving the sensitivity of ET. In any case, 3G detectors imply a substantial change in GW astronomy. While 2G detectors will only be able to reach up to for BNS and for BBHs, 3G detectors might reach for BNS and for BBHs. In terms of multimessenger events, this corresponds to thousands or tens of thousands standard sirens.
The sky localization of events in 3G will vary depending on the available network of detectors Mills et al. (2018). In this sense, there are already plans to upgrade advanced LIGO detectors. This envisioned upgrade is known as LIGO Voyager et al. (2017). Voyager could reach sensitivities between 2G and 3G. The localization will thus vary depending on the redshift of the source since the sensitivity of the network will not be homogeneous. A network of three Voyager detectors plus ET would localize of the events within , while a setup with three ET detectors would localize of the events within Mills et al. (2018).
Moreover, spacebased GW detectors have been also projected. The European space agency has approved LISA Audley et al. (2017). Being in space and with million kilometer arms, the frequency band and targets of LISA are very different from groundbased detectors (see Fig. 7). Expected sources include supermassive BHs, extreme mass ratio inspirals (EMRI) and some already identified white dwarf binaries (known as verification binaries). It is presumed that these sources could be observed with counterparts, enlarging the reach of multimessenger GW astronomy. For reference, we have included in Fig. 7 the expected background of massive BBH () and unresolved galactic whitedwarf binaries Moore et al. (2015) (see more details about the different sources in Fig. 1 of Audley et al. (2017)).
Finally, there are other proposals to detect GWs at even lower frequencies, in the band of 1100 nHz. Sources in this regime could be binary SMBH in early inspiral or stochastic, cosmological backgrounds. These GWs could be observed using a network of millisecond pulsar, in which the pulsation is extremely wellknown, for instance with PPTA Zhu et al. (2014). Other proposals are to use astrometry with GAIA, which is capable of tracking the motion of a billion stars Moore et al. (2017), or to use radio galaxy surveys Raccanelli (2017).
Iv Standard sirens
GWs coming from distant sources can feel the cosmic expansion in the same way as EM radiation does. In fact, we have seen in Sec. III.2 that the amplitude of the GWs is inversely proportional to the GW luminosity distance . In GR the GW luminosity distance is equal to EM luminosity distance, with the standard formula given by (52). However, this is not a universal relation in theories beyond GR as we will discuss in Sec. VI. For the moment, we will restrict to Einstein’s theory only.
In order to measure distances in cosmology one needs both a time scale and a proper ruler. The inverse dependence of the strain with makes GWs natural cosmic rulers. Introducing the full cosmological dependence^{7}^{7}7In (52) we had assumed a flat universe., the GW luminosity distance is given by
(59) 
where for a positive, zero and negative spatial curvature respectively. Assuming a CDM cosmology, the Hubble parameter is a function of the matter content , the curvature and the amount of DE (radiation at present time is negligible)
(60) 
On the contrary, GWs alone do not provide information about the source redshift. This is because gravity cannot distinguish a massive source at large distances with a light source at short distances. Nevertheless, when GWs events are complemented with other signals that allow a redshift identification, these events become standard sirens Schutz (1986). Standard sirens are complementary to already wellestablished standard candles, SN events in which the intrinsic luminosity can be calibrated allowing for a measurement of the EM luminosity distance. There are also standard rulers, such as the one determined by baryon acoustic oscillations (BAO) which provides the angular diameter distance. For binary blackholes (BBH) it is not expected to observe any counterpart unless there is matter around the BHs Perna et al. (2016). Fortunately, binary neutron stars (BNS) and blackhole neutron star systems (BHNS) are expected to emit short gammaray bursts (sGRB) and other EM counterparts, becoming clear standard siren targets.
The first ingredient for a standard siren is the measurement of the GW luminosity distance. However, is degenerate with the inclination of the binary. More precisely, showing the explicit angular dependence of the waveform (34) one finds that the two polarizations scale as
(61) 
where is the inclination angle. This distanceinclination degeneracy is the main source of uncertainty of present measurements of Abbott et al. (2016c). One possibility to break this degeneracy is to have an identification of both polarizations. This requires at least a three detector network and a good sky localization. Another possibility to break this degeneracy is when the binary has a precessing spin. Then, there is a characteristic modulation of the amplitude that can disentangle the inclination angle. Orbital precession is more significant for large effective spin ^{8}^{8}8The effective spin is the mass weighted projection of the two spins of the binaries into the orbital angular momentum. and/or small mass ratios since there is also an effective spinmass ratio degeneracy. Possibly good candidates for this would be BHNS binaries since BNS typically have a mass ratio close to 1.
The other ingredient for a standard siren is the identification of the redshift. This can be achieved by different means. The simplest consists in finding an EM counterpart of the GWs from the binary coalescence Schutz (1986). Then, the redshift could be extracted from the EM counterpart or from the host galaxy depending on the case. BNS will produce a sGRB after the merger. This sGRB is characterized by a beaming angle , which is typically expected to be . This means that depending on the orientation of the source we will be able to detect both signals only in a small fraction of the events. Observing a bright afterglow or kilonovae Metzger (2017) might increase the changes of detecting a counterpart. BNS will be the primary source for LIGO Dalal et al. (2006), although BHNS could also play an important role Vitale and Chen (2018). SMBHs might be good standard sirens for LISA as well Holz and Hughes (2005). Several multimessenger observations will lead to a precise measurement of the cosmic expansion either for second generation detectors Nissanke et al. (2010, 2013) or for third generation Sathyaprakash et al. (2010).
There are alternative proposals to identify the redshift without observing a counterpart. Based on statistical methods, one could associate every GW event with all the galaxies within the error in the localization and compute the cosmology Schutz (1986); Del Pozzo (2012). For a large number of events, the true cosmology will statistically prevail. Conveniently, this method applies to any type of source, including BBH which is the most common observation. Moreover, for very loud (golden) events there might be only few galaxies in the localization box Chen and Holz (2016). On the con side, this method relies on a complete galaxy catalogue.
For events involving a NS there are other possibilities. If the EoS of the NS is known, one could compute the tidal effects in the GW phase, which breaks the degeneracy between the source masses and the redshift Messenger and Read (2012). A good sensitivity could be achieved with Einstein Telescope Del Pozzo et al. (2017). Since this method relies on the knowledge of the EoS, which most probably will be uncover through GW observations also, an iterative approach could be performed. In addition, one could benefit from the narrow mass distribution of NS to statistically infer the redshift Taylor et al. (2012). Finally, numerical simulations suggests that in BNS a short burst of GWs with a characteristic frequency will be emitted after the merger. If this burst was observed, a redshift measurement could be obtained Messenger et al. (2014). The main challenge of this method is possibly the low SNR of the GW burst.
GW170817 has become the first standard siren detected. The redshift, , was obtained identifying the host galaxy NGC4993 through the different EM counterparts Abbott et al. (2017f). For such a close event, only the leading term in the cosmic expansion could be obtained following (53). The precise value obtained was Abbott et al. (2017g),
(62) 
This result has the relevance of being the first independent measurement of using GWs. Still, since it is only one event, the relative error is large, of the order of . From this error budget, arises from the uncertainty in the measurement of the distance due to present detector sensitivity and the previously mentioned degeneracy with the inclination angle. The rest of the error comes from the uncertainty in the estimation of the peculiar velocity of the host galaxy. Observations of the afterglow in different frequencies can help in reducing the inclination uncertainty Guidorzi et al. (2017); Hotokezaka et al. (2018). One could also use the statistical method to obtain without information of the counterpart, although the error is significantly larger Fishbach et al. (2018). Recent studies have shown that with order BNS standard sirens events could be measured at the level of Chen et al. (2017); Feeney et al. (2018a). Depending on the actual population of BNS this might be achieved with second generation detectors. LISA will detect mergers of SMBHs (with EM counterparts), providing measurements of cosmic expansion up to and potentially measuring with precision Tamanini et al. (2016).
iv.1 The Hubble rate tension
Standard siren observations of the cosmic expansion can also explore the tension on the Hubble parameter: where a distance ladder measurement gives a value