Unveiling the Dynamics of the Universe

We explore the dynamics and evolution of the Universe at early and late times, focusing on both dark energy and extended gravity models and their astrophysical and cosmological consequences. Modified theories of gravity not only provide an alternative explanation for the recent expansion history of the universe, but they also offer a paradigm fundamentally distinct from the simplest dark energy models of cosmic acceleration. In this review, we perform a detailed theoretical and phenomenological analysis of different modified gravity models and investigate their consistency. We also consider the cosmological implications of well motivated physical models of the early universe with a particular emphasis on inflation and topological defects. Astrophysical and cosmological tests over a wide range of scales, from the solar system to the observable horizon, severely restrict the allowed models of the Universe. Here, we review several observational probes -- including gravitational lensing, galaxy clusters, cosmic microwave background temperature and polarization, supernova and baryon acoustic oscillations measurements -- and their relevance in constraining our cosmological description of the Universe.


Introduction
During the last few decades Cosmology has evolved from being mainly a theoretical area of physics to become a field supported by high precision observational data. Recent experiments call upon state of the art technology in Astronomy and Astrophysics to provide detailed information on the contents and history of the Universe, which has led to the measurement of parameters that describe our Universe with increasing precision. The standard model of Cosmology is remarkably successful in accounting for the observed features of the Universe. However, a number of fundamental open questions remains at the foundation of the standard model. In particular, we lack a fundamental understanding of the recent acceleration of the Universe [1,2]. What is the so-called "dark energy" that is driving the cosmic acceleration? Is it vacuum energy or a dynamical field? Or is the acceleration due to infra-red modifications of Einstein's theory of General Relativity (GR)? How is structure formation affected in these alternative scenarios? What are the implications of this acceleration for the future of the Universe?
The resolution of these fundamental questions is extremely important for theoretical cosmology. Dark energy models are usually assumed to be responsible for the acceleration of the cosmic expansion in most cosmological studies. However, it is clear that these questions involve not only gravity, but also particle physics. String theory provides a synthesis of these two branches of physics and is widely believed to be moving towards a viable quantum gravity theory. One of the key predictions of string theory is the existence of extra spatial dimensions. In the brane-world scenario, motivated by recent developments in string theory, the observed 3-dimensional universe is embedded in a higher-dimensional spacetime [3]. The new degrees of freedom belong to the gravitational sector, and can be responsible for the late-time cosmic acceleration [4,5]. On the other hand, generalizations of the Einstein-Hilbert Lagrangian, including quadratic Lagrangians which involve second order curvature invariants have also been extensively explored [6][7][8][9][10]. These modified theories of gravity not only provide an alternative explanation for the expansion history of the Universe [11][12][13], but they also offer a paradigm fundamentally distinct from the simplest dark energy models of cosmic acceleration [14], even from those that perfectly mimic the same expansion history. Nevertheless, it has been realized that a large number of modified gravity theories are amenable to a scalar-tensor formulation by means of appropriate metric re-scalings and field redefinitions.
It is, therefore, not surprising that we can think of scalar-tensor gravity theories as a first stepping stone to explore modifications of GR. They have the advantage of apparent simplicity and of a long history of examination. First proposed in its present form by Brans and Dicke for a single scalar field [15], they have been extensively generalised and have maintained the interest of researchers until the present date. For instance, an extensive field of work has been developed in the cosmological dynamics of scalar-tensor theories. This can be elegantly summarised by carrying out a unified qualitative analysis of the dynamical system for a single scalar field. We also have an established understanding of the observational bounds in these models, where we can use the Parametrized Post-Newtonian formalism to constrain the model parameters. Finally, with a conformal transformation, these theories can be recast as matter interacting scalar fields in General Relativity. In this format, they can still play an important role in dark energy modelling, such as in coupled quintessence models. The consideration of a multi-scalar fields scenario, which can be perceived as the possible reflection of a multi-scalar tensor gravity theory allows for cooperative effects between the fields yielding assisted quintessence.
Indeed, scalar fields are popular building blocks used to construct models of present-day cosmological acceleration. They are appealing because such fields are ubiquitous in theories of high energy physics beyond the standard model and, in particular, are present in theories which include extra spatial dimensions, such as those derived from string theories. Recently, relative to scalar-tensor theory, much work has been invested in the Galileon models and their generalizations [16]. The latter models allow nonlinear derivative interactions of the scalar field in the Lagrangian and lead to second order field equations, thus removing any ghost-like instabilities. The Lagrangian was first written down by Horndeski in 1974 [17], which contains four arbitrary functions of the scalar field and its kinetic energy. The form of the Lagrangian is significantly simplified by requiring specific self-tuning properties (though it still has four arbitrary functions), however, the screening is too effective, and will screen curvature from other matter sources as well as from the vacuum energy [18]. An alternative approach consists of searching for a de Sitter critical point for any kind of material content [19]. These models might alleviate the cosmological constant problem and can deliver a background dynamics compatible with the latest observational data.
A promising alternative to explain the late-time cosmic acceleration is to assume that at large scales Einstein's theory of GR breaks down, and a more general action describes the gravitational field. Thus, one may generalize the Einstein-Hilbert action by including second order curvature invariants such as R 2 , R µν R µν , R µναβ R µναβ , C µναβ C µναβ , etc. Some of the physical motivations for these modifications of gravity were inspired on effective models raised in string theory, which indeed may lead to the possibility of a more realistic representation of the gravitational fields near curvature singularities [20]. Moreover, the quantization of fields in curved space-times tell us that the high-energy completion of the Einstein-Hilbert Lagrangian of GR involves higher-order terms on the curvature invariants above [21]. This is in agreement with the results provided from the point of view of treating GR as an effective field theory [22]. Among these extensions of GR the so-called f (R) gravity has drawn much attention over the last years, since it can reproduce late-time acceleration and in spite of containing higher order derivatives, it is free of the Ostrogradsky instability, as can be shown by its equivalence with scalar-tensor theories (for a review on f (R) gravity see Refs. [6][7][8][9]). Moreover, f (R) gravities have been also proposed as solutions for the inflationary paradigm [23], where the so-called Starobinsky model is a successful proposal, since it satisfies the latest constraints released by Planck [24]. In addition, the equivalence of f (R) gravities to some class of scalar-tensor theories has provided an extension of the so-called chameleon mechanism to f (R) gravity, leading to some viable extensions of GR that pass the solar system constraints [25,26]. Other alternative formulations for these extensions of GR have been considered in the literature, namely, the Palatini formalism, where metric and affine connection are regarded as independent degrees of freedom, which yields an interesting phenomenology for Cosmology [27]; and the metric-affine formalism, where the matter part of the action now depends and is varied with respect to the connection [28]. Recently, a novel approach to modified theories of gravity was proposed that consists of adding to the Einstein-Hilbert Lagrangian an f (R) term constructed a la Palatini [29]. It was shown that the theory can pass the Solar System observational constraints even if the scalar field is very light. This implies the existence of a long-range scalar field, which is able to modify the cosmological and galactic dynamics, but leaves the Solar System unaffected.
Note that these modified theories of gravity are focussed on extensions of the curvature-based Einstein-Hilbert action. Nevertheless, one could equally well modify gravity starting from its torsion-based formulation and, in particular, from the Teleparallel Equivalent of General Relativity (TEGR) [30]. The interesting point is that although GR is completely equivalent with TEGR at the level of equations, their modifications (for instance f (R) and f (T) gravities, where T is the torsion) are not equivalent and they correspond to different classes of gravitational modifications. Hence, f (T) gravity has novel and interesting cosmological implications, capable of describing inflation, the late-time acceleration, large scale structure, bouncing solutions, non-minimal couplings to matter, etc [31][32][33].
Another gravitational modification that has recently attracted much interest is the massive gravity paradigm, where instead of introducing new scalar degrees of freedom, such as in f (R) gravity, it modifies the graviton itself. Massive gravity is a well-defined theoretical problem on its own and has important cosmological motivations, namely, if gravity is massive it will be weaker at large scales and thus one can obtain the late-time cosmic acceleration. Fierz and Pauli presented the first linear massive gravity. However, it was shown to suffer from the van Dam-Veltman-Zakharov (vDVZ) discontinuity [34,35], namely the massless limit of the results do not yield the massless theory, namely, GR. The incorporation of nonlinear terms cured the problem but introduced the Boulware-Deser (BD) ghost. This fundamental problem puzzled physicists until recently, where a specific nonlinear extension of massive gravity was proposed by de Rham, Gabadadze and Tolley (dRGT), in which the BD ghost is eliminated by a Hamiltonian constraint [16]. This new nonlinear massive gravity has interesting cosmological implications, for instance, it can give rise to inflation, late-time acceleration [16]. However, the basic versions of this theory exhibit instabilities at the perturbative level, and thus suitable extensions are necessary. These could be anisotropic versions, f (R) extensions, bigravity generalizations, partially-massive constructions. The crucial issue is whether one can construct a massive gravity and cosmology that can be consistent as an alternative to dark energy or other models of modified gravity, and whether this theory is in agreement with high-precision cosmological data, such as the growth-index or the tensor-to-scalar ratio, remains to be explored in detail.
Quantum field theory predicts that the universe underwent, in its early stages, a series of symmetry breaking phase transitions, some of which may have led to the formation of topological defects. Different types of defects may be formed depending on the (non-trivial) topology of the vacuum manifold or the type of symmetry being broken. For instance, Domain Walls-which are surfaces that separate domains with different vacuum expectation values-may arise due to the breaking of a discrete symmetry, whenever the vacuum manifold is disconnected. Line-like defects, or Cosmic strings, are formed if the vacuum is not simply connected or, equivalently, if it contains unshrinkable loops. This type of vacuum manifold results, in general, from the breaking of an axial symmetry. Moreover, if the vacuum manifold contains unshrinkable surfaces, the field might develop non-trivial configurations corresponding to point-like defects, known as Monopoles. The spontaneous symmetry breaking of more complex symmetry groups may lead to the formation of textures, delocalized topological defects which are unstable to collapse. The production of topological defects networks as remnants of symmetry breaking phase transitions is thus predicted in several grand unified scenarios and in several models of inflation. Moreover, recent developments in the braneworld realization of string theory suggest that its fundamental objects-p-dimensional D-branes and fundamental strings-may play the cosmological role of topological defects.
Topological defects networks, although formed in the early universe, may in most instances survive throughout the cosmological history and leave a variety of imprints on different observational probes. The observational consequences of topological defect networks can be very diverse, depending both on the type of defects formed and on the evolution of the universe after they are generated. Although the possibility of a significant contribution to the dark energy budget has been ruled out both dynamically [36] and observationally [37], light domain walls may leave behind interesting astrophysical and cosmological signatures. For instance, they may be associated to spatial variations of the fundamental couplings of nature (see, e.g., [38]). On the other hand, cosmic strings may contribute significantly to small-scale cosmological perturbations and have consequently been suggested to have significant impact on the formation of ultracompact minihalos [39], globular clusters [40], super-massive black holes [41] and to provide a significant contribution to the reionization history of the Universe [42]. Both cosmic strings and domain walls may be responsible for significant contributions to two of the most significant observational probes: the temperature and polarization anisotropies of the cosmic microwave background and the stochastic gravitational wave background. This fact-alongside the possibility of testing string theory through the study of topological defects-greatly motivates the interest on the astrophysical and cosmological signatures of topological defects.
An extremely important aspect of modern cosmology is the synergy between theory and observations. Dark energy models and modified gravity affect the geometry of the universe and cosmological structure formation, impacting the background expansion and leaving an imprint on the statistical properties of the large-scale structure. There are a number of well-established probes of cosmic evolution, such as type Ia supernovae, baryon acoustic oscillations (BAO), weak gravitational lensing, galaxy clustering and galaxy clusters properties [43]. Different methods measure different observables, probing expansion and structure formation in different and often complementary ways and have different systematic effects. In particular, joint analyses with Cosmic Microwave Background (CMB) data are helpful in breaking degeneracies by constraining the standard cosmological parameters. Indeed, CMB has revolutionized the way we perceive the Universe. The information encoded in its temperature and polarization maps provides one of the strongest evidences in favour of the hot Big-Bang theory and has enabled ways to constrain cosmological models with unprecedented accuracy [44].
The CMB also encodes additional information about the growth of cosmological structure and the dynamics of the Universe through the secondary CMB anisotropies. These are originated by physical effects acting on the CMB after decoupling [45], such as the integrated Sachs-Wolfe effect and the Sunyaev-Zel'dovich (SZ) effect, manifest respectively on the largest and arc-minute scales of the CMB. In this review, we will discuss in some detail both a well-established acceleration probe (weak lensing) and a few promising ones related to galaxy cluster properties and the SZ effect. Galaxy clusters are the largest gravitationally bound objects in the Universe and are among the latest bound structures forming in the Universe. For this reason, their number density is highly sensitive to the details of structure formation as well as to cosmological background parameters and cluster abundance is a well established cosmological probe. The Sunyaev-Zel'dovich effect [46,47] is the scattering of CMB photons by electrons in hot reservoirs of ionized gas in the Universe, such as galaxy clusters. In particular, SZ galaxy cluster counts, profiles, scaling relations, angular power spectra and induced spectral distortions are promising probes to confront model predictions with observations. Weak gravitational lensing [48] describes the deflection of light by gravity in the weak regime. Its angular power spectrum is a direct measure of the statistical properties of gravity and matter on cosmological scales. Weak lensing, together with galaxy clustering, is the core method of the forthcoming Euclid mission to map the dark universe. Euclid [49] will provide us with weak lensing measurements of unprecedented precision. To obtain high-precision and high-accuracy constraints on dark energy or modified gravity properties, with both weak lensing and SZ clusters, non-linearities on structure formation must be taken into account. While linear covariant perturbations equations may be evolved with Boltzmann codes, non-linearities require dedicated N-body cosmological simulations [50]. There currently exist a number of simulations for various modified gravity and dark energy models, together with a set of formulas that fit a non-linear power spectrum from a linear one. Hydrodynamic simulations, commonly used in cluster studies, are increasingly needed in weak lensing applications to model various baryonic effects on lensing observables, such as supernova and AGN feedback, star formation or radiative cooling [51].
The ΛCDM framework provides a very good fit to various datasets, but it contains some open issues [52]. As an example, there are inconsistencies between probes, such as the tension between CMB primary signal (Planck) and weak lensing (CFHTLenS) [53], as well as problems with the interpretation of large-scale CMB measurements (the so-called CMB anomalies) [54]. Alternatives to ΛCDM or deviations to General Relativity are usually confronted with data using one of two approaches: model selection or parameterizations. In model selection a specific model is analyzed and its parameters constrained. Such analyses have a narrower scope but may be better physically motivated. Parameterizations are good working tools and are helpful in highlighting a particular feature and in ruling out larger classes of models, however they must be carefully defined in a consistent way. Parameterizations are commonly applied to the dark energy equation of state and to deviations from General Relativity. An example of the latter is the gravitational slip, which provides an unambiguous signature of modified gravity, and can be estimated combining weak lensing measurements of the lensing potential with galaxy clustering measurements of the Newtonian potential. Model parameters in both of the approaches are usually estimated with Monte Carlo techniques, while the viability of different models may be compared using various information criteria.
Besides model testing, cosmological data is also useful to test foundational assumptions, such as the (statistical) cosmological principle and the inflationary paradigm. The common understanding is that cosmological structures are the result of primordial density fluctuations that grew under gravitational instability collapse. These primordial density perturbations would have originated during the inflationary phase in the early universe. Most single field slow-roll inflationary models produce nearly Gaussian distributed perturbations, with very weak possible deviations from Gaussianity at a level beyond detection [55][56][57]. However, non-Gaussianities may arise in inflationary models in which any of the conditions leading to the slow-roll dynamics fail [58], such as the curvaton scenario [59][60][61], the ekpyrotic inflacionary scenario [62,63], vector field populated inflation [64][65][66] and multi-field inflation [67][68][69][70]. Tests of non-Gaussianity are thus a way to discriminate between inflation models and to test the different proposed mechanisms for the generation of primordial density perturbations. Likewise, the assumption of statistical homogeneity may be tested. Locally, matter is distributed according to a pattern of alternate overdense regions and underdense regions. Since averaging inhomogeneities in the matter density distribution yields a homogeneous description of the Universe, then the apparent homogeneity of the cosmological parameters could also result from the averaging of inhomogeneities in the cosmological parameters, which would reflect the inhomogeneities in the density distribution. The theoretical setup closest to this reasoning is that of backreaction models, where the angular variations in the parameters could also source a repulsive force and potentially emulate cosmic acceleration. Hence the reasoning is to look for these inhomogeneities, not in depth, but rather across the sky and then to use an adequate toy model to compute the magnitude of the acceleration derived from angular variations of the parameters compared to the acceleration driven by a cosmological constant [71].
This work will focus on all of the above-mentioned topics. More specifically, this article is outlined in the following manner: In Section 2, we present scalar-tensor theories, and in Section 3, we consider Horndeski theories and the self-tuning properties. In Section 4, f (R) modified theories of gravity and extensions are reviewed. In Section 5, an extensive review on topological defects is carried out. The following sections are dedicated to observational cosmology. In particular, in Section 6, cosmological tests with galaxy clusters at CMB frequencies are presented and in Section 7, gravitational lensing will be explored. In Section 8, the angular distribution of cosmological parameters as a measurement of spacetime inhomogeneities will be presented. Finally, in Section 9 we conclude.

Scalar-tensor Theories
Scalar-tensor theories of gravity stem from the original proposal of Brans-Dicke theory (BD) [15], although theories involving a scalar field with a gravitational role in addition to the tensor metric fields can be traced back to the earlier Kaluza-Klein type theories upon dimensional reduction [72]. One notorious example was Jordan's proposal for a field theory realization of Dirac's Large Number Hypothesis [73]. In the course of time, it has been realized that a large number of modified gravity theories are amenable to a scalar-tensor formulation [74] by means of appropriate metric re-scalings and field redefinitions. Interestingly, this possibility of a scalar-tensor description occurs regardless of the field equations being derived through the metric or, alternatively, the Palatini variational prescription [29]. For instance, it is well known that the f (R) higher-order gravity theories can be cast as a scalar-tensor theory [74][75][76]. The same happens for the theories based on Lagrangeans with non-linear kinetic terms of a scalar field dubbed k-essence [77].

General Formalism
The general scalar-tensor class of gravity theories encompassing BD's and a plethora of modified theories that we just alluded to, were first examined by Bergmann [78], Wagoner [79], and Nordvedt [80] (BWN). Their fundamental action can be cast in the following form [81] where R is the usual Ricci curvature scalar of the spacetime, φ is a scalar field, ω(φ) is a dimensionless function of φ which calibrates the coupling between the scalar field and gravity (in what follows we shall refer to it as the coupling parameter), λ(φ) is another function of the scalar field which can be interpreted both as a potential for φ and as a cosmological parameter and, finally, S M represents the action for the matter fields. It is worth mentioning at this point that the action should include an additional term involving the extrinsic curvature, as in GR [82], to account for a boundary term associated with the variation of R αβ with respect to g αβ . A discussion of this term in the scalar-tensor theories is provided in [74]. It becomes immediately clear that the scalar field φ plays the role which in Einstein's GR is confined to the gravitational constant, and has dimensions of a mass squared, φ ∼ 1/G. Since φ is now a dynamical variable, we thus see that the trademark of these theories is the fact that they exhibit a varying gravitational "constant". In the archetypal Brans-Dicke theory the coupling ω(φ) is constrained to be a constant. The consideration of couplings between a scalar field and the spacetime geometry envisaged by this class of theories allows a generalization to multiple scalar fields so that we may consider multi-scalar-tensor theories. The Lagrangian can then be cast as where R is the usual Ricci curvature scalar of a space-time endowed with the metric g αβ , φ µ (µ = 1, . . . , n) represents n scalar fields, f (φ µ ) is a function of those scalar fields, ω µν (φ β ) is the metric of the internal space of the φ µ fields, V(φ µ ) is the scalar potential. The Lagrangian of the matter fields L m does not carry any explicit dependence on φ µ which guarantees that free-falling test particles follow the geodesics of the spacetime (this amounts to keeping valid the conservation of the energy-momentum tensor ∇ β T αβ = 0). We leave aside a detailed consideration of this extended class of theories, which have been addressed in [83,84]. However, in subsection 2.4 we shall analyse a model of quintessence stemming from the cooperative behaviour of several scalar fields coupled to several matter components which as it will be argued can be interpreted as a particular realisation of the models characterised by Equation (2). Taking the variational derivatives of the action (1) with respect to the two dynamical variables g αβ and φ , yields the field equations where T ≡ T γ γ is the trace of the energy-momentum tensor, T α β , of the matter content of spacetime, and G N is the gravitational constant normalized to its present value (in what follows in this section, we set G N = 1 ).
From these equations the role of λ(φ) becomes more apparent. As pointed out in [78], this term enters equations (3) as a cosmological function, but intervenes in the second equation (4) as a mass term. A consequence of this latter feature is that the scalar waves propagate with a speed different from the speed of the electromagnetic waves. Another interesting aspect to point out in connection with this term is that the inverse square law for the static gravitational field is modified in the presence of a non-vanishing λ(φ) . This manifests itself in that a Yukawa-type term arises [79,81].
It is important to notice two additional points. First, we assume that the usual relation ∇ β T αβ = 0 establishing the conservation laws satisfied by the matter fields holds. This is ensured by requiring that the matter fields Lagrangian L m be independent of φ [78] The role of the scalar field is then that of helping to generate the spacetime curvature associated with the metric. Matter may create this field, but the latter cannot act back directly on the matter, which thus responds only to the metric [85].
The second point follows from Equation (4). By allowing the coupling parameter ω(φ) to tend to infinity, it can be seen that the right-hand side of the equation may vanish. This will be the case when ω → ∞ dominates over ω (φ) , λ(φ) and λ (φ) [86]. Indeed, in the absence of λ(φ), one requires ω /w 3 → 0 [80]. A solution with constant φ is then asymptotically admitted by the Klein-Gordon equation, and thus GR is recovered in the ω → ∞ limit of scalar-tensor theories (see also [83]). However, it is important, on the one hand, to separate the cases where the trace of the energy-momentum tensor is trivially vanishing, i.e., vacuum and radiation. In the latter cases, and again when there is no cosmological potential λ(φ), it has been argued that there are asymptotic behaviors that differ from GR [83,[86][87][88]. On the other hand, in the presence of the potential λ(φ), the asymptotic behaviour is very much dependent on it. In Ref. [86], it was shown that when the potential is of a power-law nature, in particular, when φ λ(φ) ∝ φ 2 , it governs the approach to GR, superseding the mechanism due to the divergence of ω(φ). Morevover, as shown below, BD theory can also be an attractor [89].

Conformal Picture
A convenient device to deal with modified gravity theories is to make use of an appropriate conformal transformation [74,75,82] to bring the theory into the form of GR plus a minimally coupled field. Consider the transformationg where φ 0 is a constant with the dimensions of mass squared which we introduce to make the scaling factor dimensionless. In order to keep the next results as simple as possible, we shall assume that φ 0 = 1 and that φ in the following expressions is dimensionless (i.e., whenever we write φ it corresponds to φ/φ 0 where φ 0 is some value of reference, for instance φ 0 = 1/G). Using the transformed metric g αβ and redefining the scalar field as we find that the action adopts the form where V(ϕ) = U (φ(ϕ)) /φ 2 (ϕ), and alsoL m = φ −2 L m . Thus, we realize that the action of the STT becomes the action of GR with matter and a massive scalar field. Writing the Einstein field equations in the transformed frame (dubbed the Einstein frame), we obtaiñ It becomes immediately apparent that, when the trace of the energy-momentum tensor of the matter fields in the original frame is vanishing, there is an equivalence between the two pictures related via the conformal transformation. This happens trivially for a vacuum model, and also in the radiation case. In all other cases, the GR-like frame exhibits a coupling between matter and the redefined scalar field. One is then drawn to conclude that the equivalence principle is violated. However, since this variation is induced by the presence of the scalar field which also acts as a source of the Einstein equations (8) the latter conclusion seems arguable and has been much debated [90][91][92]. If one considers the scalar field as a component of the total energy-momentum tensor on the right-hand side of the field equations, then the realization that the original matter does not exhibit the geodesic behaviour is a mere consequence of the fact that it couples to the other source component.
At this stage we should mention that the conformal transformation mixes the geometric and matter degrees of freedom, which results in many interpretational ambiguities [93]. One may also mention that if one only restricts attention to the Einstein frame, one may also lose sight of the original motivations and modifications of gravity in the geometrical sector [94,95]. However, the Einstein and Jordan frames are physically equivalent, as can be traced back to Dicke's original paper [96], where the conformal transformation technique was introduced. In fact, both conformal frames are equivalent, in the spirit of Dicke's paper, provided that in the Einstein frame the units of mass, time and space scale as appropriate powers of the scalar field, and are thus varying. We will discuss the issues of discriminating the geometric and matter sectors in Section 4.1.2.

Cosmological Dynamics
Here we consider the dynamics of cosmological models in general scalar-tensor gravity theories. A number of important, exact solutions of Brans-Dicke models [97,98], and of more general class of theories [99][100][101] have been derived in the literature. In what follows we resort to a unified qualitative analysis of the major dynamical features of these scalar field cosmological models [86,89,[102][103][104][105] which has the advantage of relying on a closed, autonomous system of equations where the functional dependence on φ of both the coupling ω(φ) and cosmological potential λ(φ) is left unspecified. The analysis of the relevant features of the dynamics, namely, the set of possible asymptotic behaviors, leads to a classification of the classes of couplings and potentials vis-à-vis to their cosmological impact. Alternative dynamical system analyses that are found in the literature usually envisage specific classes of scalar-tensor theories [106][107][108] or, when more general [109][110][111][112], are not so transparent as the treatment that follows.
In order to investigate the cosmological dynamics, it is most convenient to introduce dimensionless expansion normalized variables akin to the square roots of density parameters. This associates to the equilibrium points of the dynamical system (also dubbed singular, critical, or fixed points in the jargon of the theory of dynamical systems) non-static solutions of the cosmological evolution. An alternative qualitative approach based on the bare dimensional variables of the cosmological model, namely, on the scale factor of the universe, Hubble factor and/or on the energy densities of the various source components is also possible, but is of more limited scope in the scrutiny of non-static asymptotic regimes. It is however required whenever the coupling function ω(φ) or the cosmological potential λ(φ) have vanishing minima.
In what follows, we shall restrict to the homogeneous and isotropic universes given by the Friedmann-Robertson-Walker (FRW) metric where k = 0, ±1 distinguishes the curvature of the spatial hypersurfaces.
In the Einstein frame the field equations with a matter content given by a perfect fluid with p = (γ − 1)ρ can be cast as [86] whereã = √ Φ/Φ * a, the overdots stand for the derivatives with respect to the conformally , where µ 0 = 8πρ 0 sets the arbitrary initial condition for the energy density of the perfect fluid, and, more importantly, ∂ ϕ ln m(ϕ) ∝ (4 − 3γ)α(ϕ) where α = ( √ 2Ω + 3) −1 is the PPN function that translates the coupling between the scalar field and matter in the Einstein frame as used in Damour and Nordtvedt [83].
These equations are those of GR with a scalar field subject both to a time-independent potential, V(ϕ), and to a time-dependent potentialm(ϕ) a −3γ illustrating the fact that in the Einstein frame the masses of particles vary. (Hereafter we shall set 8π/φ * = 1).
Introducing the new time variable N = ln a and the dimensionless density parameters as well as the expansion normalized curvature term K = k/(aH) 2 , the Einstein field equations become a fourth order, autonomous dynamical system where we have used ρ/3H 2 = Ω m = 1 − x 2 − y 2 + K. GR models correspond to the case where ∂ ϕ ln m(ϕ) = 0, and Brans-Dicke models are characterised by an exponential coupling m, i.e., by ∂ ϕ ln m(ϕ) = α 0 (as well as by V(ϕ) = 0 in BD original version, but we do not require it here). The crucial point regarding the qualitative study of general models with scalar fields lies in the ϕ -equation, since it allows the consideration of arbitrary choices of V(ϕ) and of m(ϕ) [102]. We compactify the phase space of (15-17) by considering ϕ ∈ ∪ {∞}, and distinguish the fixed points arising at finite values of ϕ, which require x = 0, from those at ϕ = ∞ (which we shall denote ϕ ∞ ) that require either x = 0 or ψ = ϕ −1 = 0. The K = 0 and the y = 0 subspaces are invariant manifolds.
For K = 0 the system reduces to the Equations (15), (16) and 18) and we find the following fixed points, namely, at finite ϕ: This case corresponds to de Sitter solutions, dominated by the scalar field, and arise at non-vanishing maxima or minima of the potential V. Notice that ∂ ϕ m may be different from 0. These solutions are attractors at minima of V and repellors at maxima of V. • x = y = 0, ϕ 0 where∂ ' m = 0. The second case corresponds to matter dominated solutions with V = 0 at maxima or minima of m. Note however that the system (15)-(18) is singular when V = 0, translating the fact that the variables in use are not regular. In the original variables, this second class of solutions exists provided that V has also an extremum, or is identically zero.
At ϕ ∞ , the finite ϕ solutions may exist also at ϕ = ∞ provided that V or m are asymptotically flat. However, other solutions may show up at ϕ = ∞, if both V and m are asymptotically exponential. In fact, in this case the equilibrium in ψ = 1/ϕ is trivially satisfied at ψ = 0 [102]. Defining W = • x V ± ∞ = ±1, y V ± = 0, that lie on the intersection of the invariant lines x 2 + y 2 = 1 and y = 0. These points correspond to the vacuum solutions of Brans-Dicke theory found by [98,99]. • x BM ∞ = −W/3, which lies on the invariant line x 2 + y 2 = 1, provided |W| < 3. These solutions (BM) correspond to those found in [114].
Z−W , that lies in the interior of the phase space domain, provided that These are scaling solutions (S) (see [89,113,115] and references therein)  For the K = 0 case, as K = 0 for K = 0, we see from the linearization of the system that which shows that the K = 0 subspace is stable whenever Q < 1 which is precisely the condition for inflationary behaviour. Thus inflation is required to guarantee both the stability of the K = 0 asymptotic solutions and the attraction to GR, whenever the latter applies. As previously referred, when the potential V(φ) yields vanishing minima in the Einstein frame the expansion normalized variables adopted here are not suitable, as they yield a dynamical system that is not defined at those minima. In such a case it is preferable to use variables such as the scale factor, the kinetic energy of the scalar field and its potential, as done in [86], or in some of the references in [106,109].
Summing up the qualitative analysis of the dynamical system (15)- (17) associates exact solutions to a classification of the fixed points, and also allows the consideration of how extended gravity theories dynamically relate to GR [89]. It reveals the interplay of the two main mechanisms that yield general relativity as a cosmological attractor of scalar-tensor gravity theories: (i) the vanishing of the coupling function α, and (ii) the existence of a minimum of the scalar field potential V. The latter mechanism is shown to supersede the former one. Moreover the approach to GR is then characterized by a de Sitter inflationary behavior. This guarantees the stability of the k = 0 GR solutions. We have also shown that at ϕ → ∞ there is a variety of power-law solutions that correspond to Brans-Dicke behavior and that can be attractors of the cosmological dynamics.
In what follows we will address a multifield scenario, which can be understood as a particular case of the multi-scalar tensor gravity models of Equation (2).

Assisted Quintessence
Assisted inflation is a model of early universe inflation where multiple scalar fields cooperate to fuel the accelerated inflationary expansion [116][117][118]. Scalar potentials that are not viable to provide inflation as a single field can be seen to work when several such fields are present in parallel. A similar idea can be implemented for our present observed accelerated expansion [119][120][121][122].
In the usual quintessence scenario, a single scalar field, uncoupled from dark matter, is responsible for the dark energy component in our universe [123]. In assisted quintessence, we have a scenario where the dark energy sector is composed of multiple scalar fields, with different scalar potentials. Since so little is known about the structure of the dark sector, there is the possibility of a coupling between the scalar fields responsible for inflation and dark matter. Coupled quintessence, where a scalar field interacts with dark matter was introduced in Refs. [124][125][126]. In [127,128] it was further suggested that the scalar field can couple in a different way to several separate dark matter species. To keep the model as general as possible, we assume the dark matter sector to consist of several matter components with different couplings to the dark energy fields [129].

General Equations
We consider an ensemble of n scalar fields φ i cross-coupled to an ensemble of m dark matter components ρ α . The cross-couplings are described by the matrix C iα where indexes i, j identify scalar field indexes and indexes α, β identify the dark matter components. The equation of motion for the fields and the various dark matter components, in a spatially flat Friedmann-Robertson-Walker metric with scale factor a(t), are then written as Notice that the ensemble of n scalar fields φ i cross-coupled to m dark matter components ρ α can be perceived as a particular case of the multi-scalar-tensor models of Equation (2). Indeed, if we set the metric of the internal space of the scalar fields to be flat, i.e., ω µν = δ µν and allow for a combination of dark matter components in the matter sector, the C iα coupling coefficients are related to the couplings α i = ∂ ln f /∂φ i of the multi-scalar-tensor theory in the Einstein frame, and we obtain the latter equations.
The solution for the dark matter component evolution can be given immediately in terms of the values of the fields as The rate of change of the Hubble function iṡ which is subject to the Friedmann constraint where It is possible to define an effective equation of state parameter for the dark sector, where .., φ n ) is the total dark energy pressure and ρ m = ∑ α ρ α the total dark matter energy density.
We are interested in looking at the critical points of the evolution to look for viable attractors of late time evolution. For this, we will consider two possible forms for the scalar potential, both leading to scaling solutions: a sum of exponential terms with V 1 (φ 1 , ..., φ n ) = M 4 ∑ i e −κλ i φ i and an exponential of a sum of terms with

Scalar Field Dominated Solution
We start our investigation with the case of a negligible matter contribution (ρ α = 0). We can show that the effective equation of state for the attractor is where λ eff can be interpreted as an effective logarithmic slope for an dynamically equivalent single field exponential potential. For the first potential, V 1 , we get that λ eff has to obey In this case incre asing the number of fields has the effect of decreasing the value of λ eff , so that an accelerated expansion can be achieved even in the case of each single field having a slope too large to fuel acceleration. For the second potential, V 2 , λ eff is given by Here we see that increasing the number of fields actually increases the value of λ eff , making an accelerated expansion more difficult.
Both of these results are similar to what is obtained for assisted inflation [116][117][118]. This was to be expected since the matter sector is being neglected.

Scaling Solution
Here we will look for fixed points where all the sectors, kinetic and potential energy densities for the scalar energy and the matter energy density, have non-vanishing contributions. We can again define an effective scalar potential logarithmic slope, λ eff . But now we also have to take into account the coupling of the scalar fields to the matter sector. This can also be represented by an effective coupling C eff . In terms of these parameters, we can obtain the effective equation of state and for the dark energy sector For the first potential, the value of λ eff is again obtained from Equation (27), whereas the value of C eff is where α can be taken to be any of the matter species. This is because the sum is constrained to yield the same result for all the matter species, otherwise there is no exact scaling solution. It is instructive to consider these quantities for the simple case of all the potentials being replicas, that is, having the same slopes and couplings. Then λ i = λ 1 and the couplings C are a diagonal matrix with C ii = C 11 . We then get where n is the number of scalar fields. The equation of state obtained is the usual result for a single coupled field, but in the scalar field energy we can see a mild dependence in the number of fields. For the second potential, the effective slope λ eff obtained in scalar dominance is no longer valid. Instead it is possible to show that we can combine the fields through a rotation so that the kinetic energy of all fields is zero except for the first field. In terms of this rotation matrix Q we get the effective equation of state as and the effective coupling to be The precise form of Q is hard to get in the general case. But the simple case of a diagonal coupling matrix C can give us an idea of what is to be expected. In this instance, we get and so Finally, if we consider again the case of n fields all with the same slope and couplings, we get λ eff = √ n λ 1 and C eff = C 11 / √ n. Here we have a strong dependence in the number of fields in the effective equation of state and the dark energy density, , This is clearly in contrast to what we have for the first potential type.

Scalar Potential Independent Solutions
Finally we look for solutions where the scalar potential energy is negligible. These solutions, therefore, are independent of the choice we make for the scalar potential. Note that the scalar field still contributes to the overall energy density through its kinetic energy. These solutions can be grouped into three different types: • First, we have the kinetic dominated solutions where all the energy density for the dark sector is due to the kinetic energy of the scalar fields. In this instance, we have w eff = Ω φ = 1.
• As a second possibility, we have a contribution from the kinetic energy and the matter energy densities for the dark sector. In this case we have w eff = Ω φ = 2 3 ∑ i C 2 iα . This will be true for any matter species α since the sum is constrained to yield the same for all matter species.
• Finally, as a third possibility we have a matter dominated scenario. This has been investigated previously in the literature for the case of a single field with two matter components [128,[130][131][132]. If we extend this to more than one field, we will see that there is a consistency condition imposed on the values of the couplings to allow a flat direction in the field space, otherwise this solution will be non existent. This is by its nature a transient solution since the matter contribution will eventually decay away and the scaling solution or the scalar field dominated solution will take over. However, this can be an interesting scenario for model building since it allows the dark matter sector to have a significant contribution before the acceleration of the Universe.
If we extend our analysis to the linear perturbations in these models we can show that there is a term proportional to ∑ i C iα C iβ that can act as a growth or damping term to the density contrast of the matter components. This yields a strong constraint on the types of models that can be compatible with observations [129].

Horndeski Theories and Self-tuning
In the previous sections we have considered the simplest scalar-tensor theories of gravity, which are given by the action (1). Those theories can be re-written as a scalar field, which is coupled to the geometry through a term f (φ)R in the action, with a canonical kinetic term and a potential. This is: where X = ∇ µ φ∇ µ φ. Nevertheless, there are more things in heaven and earth for scalar-tensor theories of gravity. K-essence models [133] assume a scalar field minimally coupled to gravity ( f (φ) = 1) with an action that is not necessarily the sum of a kinetic term and a potential but a more general function K(φ, X). One can go further and consider interaction terms in the Lagrangian containing second order derivatives of the scalar field, as in the kinetic braiding models [134] that have a scalar field Lagrangian given by K(φ, X) + G(φ, X) φ. These models can be stable because their field equations are only second order. However, theories with Lagrangians containing second order derivatives generically have field equations which contain higher than second order derivatives. Therefore, such theories usually propagate an extra (ghostly) degree of freedom which means that they are affected by the Ostrogradski instability [135]. In fact, in 1974, Horndeski found the most general scalar-tensor action leading to second order equations of motion [17]. Nevertheless, Horndeski work passed mostly unnoticed until Deffayet et al. [136] rediscovered his theory when generalizing the covariantized version [137] of the galileons models [138]. The Horndeski Lagrangian can be written as [17] where κ i (φ, X) and F(φ, X) are arbitrary functions, satisfying F ,X = κ 1,φ − κ 3 − 2Xκ 3,X . Brans-Dicke theory, k-essence, kinetic braiding and many other models can be seen as particular cases of the most general Horndeski family. Thus, although the Horndeski Lagrangian restricts the number of stable scalar-tensor theories. Actually it has been shown that there are still stable theories beyond the Horndeski Lagrangian [139][140][141]. Those theories have second order equations of motion only when one particular gauge is fixed. On the other hand, the vacuum energy also gravitates in modified theories of gravity. Thus, the cosmological constant problem [142][143][144] is still present if, as Weinberg suggested, the extra degree of freedom could screen only a given value of that constant [142]. In this context, the fab four models are based in the observation by Charmousis et al. [18,145] that Weinberg's no-go theorem can be avoided by relaxing one of the assumptions, that is allowing the field to have a non-trivial temporal dependence once the cosmological constant has been screened. Their dynamical screening is based in requiring Minkowski to be a critical point of the dynamics and, when the critical point is an attractor, it may alleviate the cosmological constant problem. Nevertheless, these models are forced by construction to decelerate its expansion while approaching the Minkowski final state. Thus, a late time accelerating cosmology does not naturally arise in this scenario. When considering the compatibility of dynamical screening the vacuum energy with a late time phase of accelerated expansion, the concept of self-adjustment was extended to non-Minkowskian final states [19]. As we will now show in detail, a scalar field self-tuning to de Sitter can lead to very promising scenarios from a phenomenological point of view [146,147]. In addition, it may alleviate the cosmological constant problem if the field is able to screen the vacuum energy before the material content in a particular model.

Dynamical Screening
In a cosmological context, one can consider a FLRW geometry to express the Horndeski Lagrangian in the minisuperspace. Once the dependence on higher derivatives is integrated by parts the point-like Lagrangian takes the simple form [18] V is the spatial integral of the volume element, H =ȧ/a is the Hubble expansion rate, and an over-dot represents a derivative with respect the cosmic time t. The functions Z i are given by where X i and Y i are written in terms of the Horndeski free functions [18]. Taking into account Equation (40), the Hamiltonian density can be written as Assuming that the matter content, described by ρ m (a), is minimally coupled with the geometry and non-interacting with the field, the modified Friedmann equation is given by Now, one can follow a similar argument as that presented in reference [18] for screening to Minkwoski, the only difference being that we require self-tuning to a more general given on-shell solution with H 2 = H 2 on = 0 (where the subscript "on" means evaluated on-shell-in-a), as done in reference [19]. Thus, assuming that the field is continuous in any phase transition that changes the value of the vacuum energy and noting that the screening of any value of the vacuum energy would lead to the screening of any material content, one has to require H on to be approached dynamical in order to have acceptable cosmological solutions. That is, we want H on to be an attractor solution, although this particular adjustment mechanism only ensures that it is a critical point. As it is subtly discussed in reference [18], this requirement leads to the following three conditions: • The field equation evaluated at the critical point has to be trivially satisfied leading the value of the field free to screen. This implies that the minisuperspace Lagrangian density at the critical point takes the form • The modified Freedman equation evaluated once screening has taken place has to depend onφ to absorb possible discontinuities of the cosmological constant appearing on the r. h. s. of this equation. Thus, taking into account Equations (43) and (44) (and itsφ-derivative), it leads to • The full scalar equation of motion must depend onä to allow for a non-trivial cosmological dynamics before screening. This leads to a condition equivalent to (45) if H on = 0.
A particular Lagrangian which satisfies conditions can be written as This Lagrangian (which is a trivial generalization of that presented in reference [18]) has a critical point at the solution H on by construction [19]. In order to obtian some information about the Z i 's, we note that, as it has been explicitly proven in [18] (it must be noted that such a proof is independent of the particular H on ), two Horndeski theories which self-tune to H on are related through a total derivative of a function µ (a, φ). Therefore, one can consider with L (a,ȧ, φ,φ) given by Equation (46). As this equation has to be valid during the whole evolution, we can equate equal powers of H to obtain which can be combined to yield [18,19] 3.2. Requiring the Existence of a de Sitter Critical Point: H 2 on = Λ In the first place, we consider k = 0. In this case the dependence on the scale factor a of the Z i 's vanishes, as can be easily noted from Equation (41), and the point-like Lagrangian (40) is independent of Y i 's. For k = 0, therefore, Equation (49) is equal to the L on (a on ) /a 3 on . In the second place, we demand H 2 on = Λ, this leads to [19] ∑ i=0..3 As the l. h. s. of this equation is independent of a, the r. h. s. should also be independent of a for any value ofφ, which allows us to obtain µ (a, φ). Thus, we have Therefore, there are three different kind of terms which can appear in the Lagrangian. These are: (i) X i -terms linear onφ; (ii) X i -terms with non-linear dependence onφ, which contribution has to vanish in the on-shell point-like Lagrangian; and, (iii) terms that are not able to self-tune because they contribute through total derivatives or terms multiplied by k in the Lagrangian. Thus, the point-like Lagrangian takes the form [19] with L linear and L nl to be given in the following discussion.

Linear Terms: "The Magnificent Seven"
In order to satisfy equation (51) considering only terms linear onφ, it is enough to have with the potentials U i and W i fulfilling the constraint The Lagrangian of the linear terms is with the potentials satisfying condition (54). As there are a total of eight functions and only one constraint, there are effectively only seven free functions which we coined 'the magnificent seven". The field equation for these models is where we have divided all the equations by a 3 for simplicity and we require that the term in the square bracket does not vanish. On the other hand, the Hamiltonian density is When only this kind of terms are present the modified Friedmann equation is H linear = −ρ m (a). These models are able to screen any material content dynamically as long as there is at least one W i = 0 with i = 0. The models with only W 0 and U i potentials and those with only U i 's would not spoil the screening of other models when combined with them, but they are not able to self-tune by themselves [19]. It is worthy to note that a non-self-tuning Einstein-Hilbert term is contained in these models and can be explicitly written redefining

Non-linear Terms
We now consider terms with an arbitrary dependence on φ andφ. These terms have a minisuperspace Lagrangian given by The functions X nl i are, however, not arbitrary. Taking into account Equation (51), the contribution of these terms has to vanish when L nl is evaluated at the critical point. Therefore, they have to satisfy the condition The Hamiltonian density of this models is Taking into account condition (59) and its derivatives, it can be seen that H nl, os depends onφ; thus, these models are always able to self-tune. If only these nonlinear terms are present, the modified On the other hand, taking into account Lagrangian (58) the field equation can be written as with ε linear + ε nl = 0.

Example of a Linear Model
The field equation and the Friedmann equation read [146] where a prime represents derivative with respect to ln a.
As an example of a linear model let us consider the three potentials U 2 , U 3 and W 2 . The constraint therefore we need: U 2,φ /W 2 = 1, during a matter domination epoch, and U 2,φ /W 2 = 4/3, for a radiation domination epoch. For example, the choice of the potentials, U 2 = e λφ + 4 3 e βφ , and W 2 = λe λφ + βe βφ , give us the desired behaviour, as shown in Figure 3. The de Sitter evolution is attained Unfortunately, the contribution of the field at early times is too large to satisfy current constraints.

Example of a Non-linear Model
We will restrict the analysis to the shift-symmetric cases, which means no dependence on φ, and make use of the convenient redefinition ψ =φ. Under these assumptions we obtain the equations of motion [147], where Q 0 , Q 1 , Q 2 , P 0 , P 1 , P 2 , are non-trivial functions of X i and H, and the average equation of state parameter of matter fluids is Let us consider a case involving the three potentials X 0 , X 1 and X 2 such that We can obtain a model with w ψ = w 0 + w a (1 − a), such that, w 0 = −0.98 and w a = 0.04, which is compatible with current limits and moreover, has a negligible dark energy contribution at early times. The evolution of the energy densities is illustrated in Figure 4.

Summary
In this section we have considered a subclass of the Horndeski Lagrangian that leads to a late-time de Sitter evolution of the Universe and that may provide a mechanism to alleviate the cosmological constant problem. We have presented examples of a linear and a non-linear model. In particular, the shift-symmetric non-linear models, the de Sitter critical point is indeed an attractor. We can, thus, understand the current accelerated expansion of our Universe as the result of the dynamical approach of the field to the critical point. The model presented satisfies current observational bounds on the background evolution and appear very promising.

f (R) Modified Theories of Gravity and Extensions
Nowadays, a great multidisciplinary effort has firmly established experimentally the standard ΛCDM model of cosmology, which is based on Einstein's General Theory of Relativity (GR), a cold dark matter source and a tiny cosmological constant. GR itself has been tested on scales ranging from the submilimeter to the solar system [81]. Nonetheless, a number of conceptual issues still refuses to naturally merge into this paradigm. Since the establishment of the theorems by Hawking and Penrose [148][149][150], space-time singularities both deep inside black holes and in the early universe have been a disturbing issue from a theoretical point of view since they represent regions where Physics looses predictability. On the other hand, introduction of a dark sector in the model allows to fit the data, but implies that up to 96% of the total energy-matter content of the universe is of unknown nature. Though several candidates for the microscopic origin of the dark matter component have been considered, like supersymmetric partners [151], this new physics has not been directly detected at particle accelerators so far. Moreover, the accumulated data in favour of a late-time accelerated expansion of the universe only hints at an exotic form of energy whose equation of state is nearly equal to that of a cosmological constant.
In view of these difficulties one may wonder whether this is a signal that the conceptual and mathematical framework represented by GR is reaching its limit of validity and that the increasing complexity of the additional ad hoc hypothesis needed to fit the data could be removed by finding a new framework for gravity where the difficulties (singularities and dark sources) are naturally encompassed. At the more fundamental level, great theoretical progress has been made within the unifying frameworks of both string theory [152] and loop quantum gravity [153] but, given the high energy/curvature scales at which the new effects would take place, they are troubled with large difficulties to make contact with experimental results. In this sense, many authors have focused their attention on phenomenological models where specific issues of the ΛCDM model (singularities, inflation, dark sources, etc) are dealt with. In the high-energy regime, these modified gravity models including higher-order powers of the curvature invariants (see e.g., [8] for a review) are supported by a high-energy completion of the Einstein-Hilbert Lagrangian of GR within the theory of quantized fields in curved space-times [21,154].

General Formalism
Recently, f (R) modified theories of gravity have been extensively explored in the literature, where the standard Einstein-Hilbert action is replaced by an arbitrary function of the Ricci scalar R [155,156]. These theories can be formulated through different approaches: (i) one could use the metric formalism, which consists in varying the action with respect to the metric g µν only; (ii) in the Palatini formalism [28,[157][158][159][160][161][162][163], where the metric and the affine connection are treated as independent variables (see Section (4.3) below for more details on this formulation); (iii) and the metric-affine formalism, where the matter part of the action now depends on and is varied with respect to the connection [28].
The action for the f (R) modified theories of gravity is given by where L m is the matter Lagrangian density, in which matter is minimally coupled to the metric g µν and ψ collectively denotes the matter fields.
Using the metric approach, by varying the action with respect to g µν yields the following field equation where Now, considering the contraction of Equation (63), provides the following relationship which shows that the Ricci scalar is a fully dynamical degree of freedom. f (R) gravity may be written as a scalar-tensor theory, by introducing a transformation {R, f } → {φ, V} defined as so that the action (62) can be written as Note that action is simply a Brans-Dicke type action (1) with the parameter ω = 0. The gravitational field equations of f (R) gravity, in the scalar-tensor representation, simply reduce to (3) and (4), with a redefinition of the potential V(φ) = 2φλ(φ).

Discriminating between Dark Energy and Modified Gravity Models
Due to the fact that subsets of scalar-tensor theories can be obtained from modifications to the gravitational sector through specific mappings, it is fundamental to understand how one may differentiate these modified theories of gravity from dark energy models. This ambiguity requires a practical classification, which is not a trivial issue as one cannot discriminate between the dark energy models and modified gravity, solely from the expansion rate of the Universe. However, as structure formation behaves differently by these two classes of theories, information on the growth of structure, at different scales and redshifts, will break the degeneracy and will serve to discriminate between both models of dark energy and modified gravity [164]. In this context, generic modifications of the dynamics of scalar perturbations, with respect to the ΛCDM background, can be represented by two new degrees of freedom in the Einstein constraint equations, namely, through the functions Q(a, k) and η(a, k), where a is the scale factor and k the perturbation scale. In modified theories of gravity, Q(a, k) results from a mass-screening effect due to local modifications of gravity, and effectively modifies Newton's constant. In the context of dynamical dark energy models, the function Q(a, k) incorporates the interaction with other fields, due to the perturbations. The function η, absent in ΛCDM, parameterizes the effective stresses due to the modification of gravity or specific dynamical dark energy models. Finally, the scale and time-dependence of both functions, Q and η, can be derived in the considered model and traced on a (Q, η) plane.
Following [164], a practical manner is to denote the term "modified gravity" if additional contributions to the Poisson equation are presented, which induces Q = 1, and if extra effective stresses arise, implying η = 1. In this context, "modified gravity" is related to models in which modifications are present in the gravitational sector and in which dark energy clusters or interacts with other fields. Following this practical classification [164], in the context of first order perturbation theory, models with Q = η = 1 are denoted standard dark energy models, such as, a minimally-coupled scalar field with standard kinetic energy [164]. On the other hand, models for Q = 1 and η = 1 are denoted "modified gravity", such as scalar-tensor theories, f (R) gravity, massive gravity and generalized galileons, Horndeski interactions, bi-(multi-) gravity, etc. Thus, in the context of the EUCLID mission [164], the definitions of the functions Q and η are extremely convenient, for instance, EUCLID can distinguish between standard dynamical dark energy and modified gravity by forecasting the errors on Q and η, and several combinations of these functions, such as Q/η.

Late-time Cosmic Acceleration
In this subsection, we show that f (R) gravity may lead to an effective dark energy, without the need to introduce a negative pressure ideal fluid. Consider the FLRW metric (11). Taking into account the perfect fluid description for matter, we verify that the gravitational field equation, Equation (63), provides the generalised Friedmann equations in the following form [165,166]: These modified Friedmann field equations may be rewritten in a more familiar form, as where ρ tot = ρ + ρ (c) and p tot = p + p (c) , and the curvature stress-energy components, ρ (c) and p (c) , are defined as respectively. The late-time cosmic acceleration is achieved if the condition ρ tot + 3p tot < 0 is obeyed, which follows from Equation (71). For simplicity, consider the absence of matter, ρ = p = 0. Now, taking into account the equation of state ω eff = p (c) /ρ (c) , with f (R) ∝ R n and a generic power law a(t) = a 0 (t/t 0 ) α [165], the parameters ω eff and α are given by respectively, for n = 1. Note that a suitable choice of n can lead to the desired value of ω eff < −1/3, achieving the late-time cosmic acceleration. It is interesting to note that observations might slightly favour the fact that the effective dark energy equation-of-state parameter lies in the phantom regime, i.e., ω eff < −1 [37], although the data are still far from being conclusive. In fact, recent fits to supernovae, CMB and weak gravitational lensing data indicate that an evolving equation of state crossing the phantom divide, is mildly favored, and several models have been proposed in the literature. In particular, models considering a redshift dependent equation of state, possibly provide better fits to the most recent and reliable SN Ia supernovae Gold dataset. In a cosmological setting, it has been shown that the transition into the phantom regime, for a single field is probably physically implausible [167], so that a mixture of various interacting non-ideal fluids is necessary. In this context, we refer the reader to [168] for a recent review on specific models in modified gravity. If confirmed in the future, this behavior has important implications for theoretical models of dark energy and modified gravity. For instance, this implies that dark energy is dynamical and excludes the cosmological constant and the models with a constant parameter, such as quintessence and phantom models, as possible candidates for dark energy. All of these models present an extremely fascinating aspect for future experiments focussing on supernovae, cosmic microwave background radiation and weak gravitational lensing and for future theoretical research.

Curvature-matter Couplings in f (R) Gravity
Recently, an extension of f (R) gravity was presented, by considering an explicit curvature-matter coupling. The latter coupling induces a non-vanishing covariant derivative of the energy-momentum tensor, ∇ µ T µν = 0, which potentially leads to a deviation from geodesic motion, and consequently the appearance of an extra force [169]. Implications, for instance, for stellar equilibrium have been studied in Ref. [170]. The equivalence with scalar-tensor theories with two scalar fields has been considered in Ref. [171], and a viability stability criterion was also analyzed in Ref. [172]. It is interesting to note that nonlinear couplings of matter with gravity were analyzed in the context of the accelerated expansion of the Universe [173][174][175], and in the study of the cosmological constant problem [176].
The action for f (R) curvature-matter coupling [169] takes the following form where f i (R) (with i = 1, 2) are arbitrary functions of the curvature scalar R. For notational simplicity we consider κ = 1 throughout this subsection.
Varying the action with respect to the metric g µν yields the field equations, given by where we have denoted F i (R) = f i (R), and the prime represents the derivative with respect to the scalar curvature. We refer the reader to [177] for the curvature-matter coupling in the Palatini formalism.
As mentioned above, an interesting feature of this theory is the non-conservation of the energy-momentum tensor. More specifically, by taking into account the generalized Bianchi identities [169,178,179], one deduces the following conservation equation where the coupling between the matter and the higher derivative curvature terms describes an exchange of energy and momentum between both.
Consider the equation of state for a perfect fluid T where ρ is the energy density and p, the pressure, respectively. The four-velocity, U µ , satisfies the conditions U µ U µ = −1 and U µ U µ;ν = 0. Introducing the projection operator h µν = g µν + U µ U ν , gives rise to non-geodesic motion governed by the following equation of motion for a fluid element: where the extra force, f µ , is given by Note that an intriguing feature is that a specific choice for the matter Lagrangian yields different dynamics. For instance, in [180], it was argued that a "natural choice" for the matter Lagrangian density for perfect fluids is L m = p, based on Refs. [181,182], where p is the pressure. This choice has a particularly interesting application in the analysis of the R−matter coupling for perfect fluids, which implies in the vanishing of the extra force [169]. However, it is important to point out that despite the fact that L m = p does indeed reproduce the perfect fluid equation of state, it is not unique [183]. Other choices include, for instance, L m = −ρ [182,184], where ρ is the energy density, or L m = −na, were n is the particle number density, and a is the physical free energy defined as a = ρ/n − Ts, with T being the temperature and s the entropy per particle (see Ref. [182,183] for details).
Hence, it is clear that no immediate conclusion may be extracted regarding the additional force imposed by the non-minimum coupling of curvature to matter, given the different available choices for the Lagrangian density. One may conjecture that there is a deeper principle or symmetry that provides a unique Lagrangian density for a perfect fluid [183]. This has not been given due attention in the literature, as arbitrary gravitational field equations depending on the matter Lagrangian have not always been the object of close analysis. See Ref. [183] for more details.

The Palatini Approach
The discussion outlined in the beginning of this section hints at reconsidering the foundational aspects of gravity as a geometric phenomenon. The key point of such an approach is to make one step back and reconsider the role that basic geometric objects play in gravitational physics. The first such step is to recognize the different role that metric and affine connection have as geometrical structures, and take them as independent objects (Palatini approach [27]), which implies the presence of non-metricity (this is actually the simplest extension towards more general geometric theories of gravitation), i.e., the failure of the independent connection to be metric, Q λµν = ∇ λ g µν = 0. Note that whether the underlying structure of spacetime is Riemannian or not is a fundamental question of gravity as a geometric phenomenon which has received little attention so far, and its associated phenomenology has been scarcely studied. Indeed, a glance at the modern description of ordered structures with defects in condensed matter physics (like Bravais crystals) teaches us that the presence of point defects induces a continuous description in terms of a geometry with independent metric and affine structures [185][186][187], which could provide important tools for gravitational physics [188,189]. In the case of GR, one finds that the Palatini approach gives exactly the same result as the standard approach of fixing the Levi-Civita connection a priori and making variations with respect to the metric g µν (metric approach); this result being attached to the particular functional form of the Einstein-Hilbert Lagrangian of GR. Nonetheless, for modified gravity it turns out that it avoids many of the problems of the metric formulation, and generically yields a set of second-order field equations and absence of ghost-like instabilities.
To fix ideas let us consider an f (R) extension of GR, defined by the standard action where κ 2 = 8πG, g is the determinant of the space-time metric g µν , which is independent of the connection Γ λ µν . The Ricci tensor, , is an object entirely constructed out of the independent connection, while the matter sector, S m = where ψ m denotes the matter fields) only couples to the metric g µν . We point out that the independence between metric and connection in the Palatini approach also allows for the presence of torsion (the antisymmetric part of the connection). However, for simplicity, in this section, we will assume that torsion vanishes, Γ [µν] = 0, and, in addition, that the Ricci tensor is purely symmetric, R [µν] = 0 (for a derivation of the field equations in full generality under the presence of torsion see [190]). Under these conditions, variation of the action (79) yields Since we are working on the Palatini approach, the terms on the variations with respect to δg µν and δR µν (Γ) must vanish independently. This brings two systems of equations of the form where f R ≡ d f /dR and the matter energy-momentum tensor T µν is defined by Equation (64). Note that in the GR case, f (R) = R, Equation (82) simply expresses the standard metric-connection compatibility condition, while Equation (81) become the Einstein equations. In the f (R) case, tracing with g µν in (81) yields R f R − 2 f = κ 2 T, which is an algebraic relation whose solution, R = R(T), generalizes the GR relation R = −κ 2 T. The field equations (81) can be solved by introducing a new rank-two tensor h µν = f R g µν in such a way that the connection equation (82) is rewritten as , the independent connection Γ λ µν becomes the Levi-Civita connection of h µν . This tensor is connected to the physical metric g µν via the conformal transformations given above. With this result, the remaining equations (81) can be written as [191] where R µ ν (h) is the Ricci tensor computed with the auxiliary metric h µν . This is a set of second-order Einstein-like equations with all the objects in the right-hand-side depending only on the matter sources (because f (R) ≡ f (R(T)) is a function of the trace of the energy-momentum tensor), which can thus be read as a modified matter source, namely, one can write G µν = τ µν , where τ µν accounts for an effective energy-momentum tensor. This property, which results from the special way gravity and matter couple in the Palatini approach, is in sharp contrast with the generic fourth order equations and ghosts that one finds in the (more conventional) metric formulation of these theories [see e.g., [192,193] for a comparison between both approaches]. Note that in the absence of matter, T µ ν = 0, the field equations boil down to those of GR with a cosmological constant term and thus these theories contain no extra propagating degrees of freedom (in particular, the theory is free of ghosts). As a consequence, the dynamics of these theories can be only excited in regions of large curvatures/short distances or for strong fields. Solving the equations (83) just requires to specify the matter source, T µ ν , and the object f (R(T)), and to transform the solution h µν back to the physical metric g µν using the conformal transformations above.

Quadratic Cosmology from Gravity-matter Coupling
This framework can be extended to include a gravity-matter coupling term g(R)L m (ψ m , g µν ) replacing the matter sector in the action (79), which yields a new set of equations [194] where now Φ ≡ f (R) + 2κ 2 g R L m , h µν = Φg µν and h µν = Φ −1 g µν . Again, these equations are manifestly second-order due to the same considerations of the previous case.
In the following, we consider a FRW cosmology, in the presence of a perfect fluid. Assuming an equation of state P = ωρ (with ω =constant), and inserting these inputs into the field equations (84) one gets two sets of equations [194] where Φ and a subindex means a derivative. Both the f (R) and the matter-gravity coupling frameworks can successfully reproduce the background evolution of loop quantum cosmology, whose chief accomplishment is the resolution of the Big Bang singularity applying non-perturbative canonical quantization techniques [195,196]. The singularity is replaced by a non-singular cosmic bounce at an energy of Planck's density order ρ c ≡ c 5 /(hG 2 ) via a quadratic cosmology model of the form with ± 1, where = +1 corresponds to the solution arising in braneworld scenarios [197] and = −1 in loop quantum cosmology [198]. In f (R) gravity, quadratic cosmology can be numerically fit with a specific model [199] (which is consistent with the fact that space-time singularities can also be avoided in spherically symmetric charged black hole scenarios of Palatini f (R) gravity [200,201]), while in f (R) theories with gravity-matter coupling different choices of Φ ≡ Φ(ρ) yields different models able to analytically reproduce this cosmic evolution [194]. This result can further be extended to isotropic and anisotropic cosmologies in Palatini gravity including quadratic R µν R µν corrections [202].

Bouncing Cosmologies in Born-Infeld-like Gravity Models
A somewhat different proposal for modified gravity models was introduced by Deser and Gibbons, coined Born-Infeld gravity [203,204], where the action is given by where λ is a small parameter, which must be suitably tuned to recover GR at low energies/curvatures. This means that, from an experimental point of view, the action (88) represents a viable theory of gravity, while encompassing new effects at high energies. Born-Infeld gravity has raised great interest due to its many astrophysical and cosmological applications, including the existence of singularity-free cosmologies.
To investigate the latter point, let us rewrite, by convenience, the first term in Born-Infeld gravity as √ −g detM, where a hat denotes a matrix andM = ĝ −1q , with q µν ≡ g µν + λ −2 R µν . Therefore we can rewrite (88) as Writing in this way, it is clear that Born-Infeld gravity (89) contains two of the invariant polynomials of the matrixM and, therefore, a natural extension of this model would be [205] where e n are the five invariant polynomials of the matrixM (see [205] for details) and β n is a set of free dimensionless parameters. In this formulation, the original Born-Infeld gravity comes from the contributions e 0 (M) = 1 and e 4 (M) = detM. When testing the non-singular character of solutions in the Palatini approach it is important to establish its robustness, namely, that the singularity resolution is not attached to a particular choice of the functions defining the model, or to some ad hoc procedure to obtain a desired result, like a fine tuning of the parameters. This program has been very successful when applied to spherically symmetric black holes replacing the point-like singularity by a wormhole structure giving continuity to all null and time-like geodesics in a number of gravitational and matter scenarios [206][207][208]. Similarly, to study bouncing solutions, let us focus on the simple extension of Born-Infeld gravity given by where we have redefined the constants to match with GR results at low energies, 1, introduced a new matrixΩ ≡M 1/2 , for notational simplicity, and added the matter contribution. As usual, in the Palatini formalism, we treat metric and connection as independent entities and obtain the field equations by variation with respect to both of them. After a bit of (lengthy) algebra one gets the field equations for the auxiliary metric t µν as [209] R µ ν (t) = with the following definitions: are functions of the auxiliary scalar (Brans-Dicke-like) field A, the auxiliary metric t µν (associated with a Levi-Civita connection, ∇ λ [ √ −tt µν ] = 0) is related to the physical metric g µν as t µν = φ(detΩ) 1/2 g µαΩ α ν , and L G ≡ (φ(detΩ) 1/2 − V(φ) −λ)/( κ 2 ) is the Lagrangian density. Note that the equation resulting from variation of Equation (91), namely, detΩ = dV/dφ, establishes an algebraic relation φ = φ(detΩ) for every choice of f (A). Moreover, due to the fact that the equations for the metric g µν can be written as this implies thatΩ =Ω(T µ ν ) and, consequently, φ = φ(T µ ν ). Therefore, all terms appearing on the right-hand-side of the field equations (92) are just functions of the matter, like in the case of f (R) theories. Thus, finding a solution of (92) for t µν immediately yields a solution for g µν .
Let us assume a perfect fluid energy-momentum tensor with N species of non-interacting fluids, where ∆ = To obtain explicit solutions one needs to further specify the family of Born-Infeld gravity extensions. Let us consider the power-law family f (detΩ) = (detΩ) n where the case n = 1/2 corresponds to the original Born-Infeld gravity. Expanding in series of 1, the corresponding Lagrangian yields the result which, via a redefinition of constants, is just GR plus a cosmological constant term. Therefore, this kind of theories is in agreement with observational data in the low energy/curvature regime, but provides new physics regarding the early universe. In this sense, at high energies two kinds of solutions are found (see [209] for details): (i) those with ρ > 0 and ω > 0 have a vanishing Hubble parameter for some finite energy density ρ B characterized by and represent unstable solutions with a minimum value at past infinity, where a(−∞) = a B ; (ii) those with ρ < 0 represent genuine non-singular bouncing solutions (as long as ω > −1) with This is in agreement with the fact that the avoidance of spacetime singularities in black holes within Born-Infeld gravity is achieved when < 0.
As the short discussion above shows, high-energy modifications of GR within the Palatini approach ( f (R), quadratic, Born-Infeld-type, and other generalizations) are able to avoid the big bang singularity via a cosmic bounce followed by a period of de Sitter expansion in a generic and robust manner. At the same time, Palatini theories generically produce an effective cosmological constant at low energies [see Equation (96)]. The specific sign and value of this constant depends both on the gravitational action and on the specific choices of parameters. This suggests to consider the combination of these two properties to achieve a model in which a sufficient small cosmological constant compatible with cosmological observations arises at late times as a result of high-energy or low-energy modifications of the theory. This could provide a unified and consistent scenario accommodating the two known de Sitter phases of the cosmic expansion without invoking dark sources of energy. On the other hand, as first suggested by Hawking [210], the existence of electrically charged microscopic black holes originated by large density fluctuations during the early universe (primordial black holes) could have a substantial impact in accounting for the matter density of the hypothesized dark matter particles. The prediction of some Palatini models of the existence of singularity-free, horizonless massive objects below a certain charge/mass scale [211], could offer new avenues for the consideration of black hole remnants [212] as a reasonable source of dark matter. Obtaining precise estimates on the abundance of such objects, as resulting from their production in the early universe, could allow to make detailed cosmological predictions according to this viewpoint.

Hybrid Metric-Palatini Gravity
The hybrid metric-Palatini gravitational theory is a novel approach to modified theories of gravity that consists of adding to the Einstein-Hilbert Lagrangian an f (R) term constructed a la Palatini [29]. Using the respective dynamically equivalent scalar-tensor representation, it was shown that the theory can pass the Solar System observational constraints even if the scalar field is very light. This implies the existence of a long-range scalar field, which is able to modify the cosmological and galactic dynamics, but leaves the Solar System unaffected. The absence of instabilities in perturbations was also verified and the theory provides explicit models which are consistent with local tests and lead to the late-time cosmic acceleration.

General Formalism
The hybrid metric-Palatini gravity theory proposed in [29], is given by the following action can be recast into a scalar-tensor representation [29,213,214] given by the action where S m is the matter action, κ 2 = 8πG/c 3 , and V(φ) is the scalar field potential. Note that the gravitational theory given by Equation (98) is similar to a Brans-Dicke scalar-tensor action with parameter w = −3/2, but differs in the coupling to curvature. The variation of this action with respect to the metric tensor provides the field equations where T µν is the matter energy-momentum tensor. The scalar field φ is governed by the second-order evolution equation which is an effective Klein-Gordon equation [29,213].

Unifying Local Tests and the Late-time Cosmic Acceleration
In the weak field limit and far from the sources, the scalar field behaves as φ(r) ≈ φ 0 + (2Gφ 0 M/3r)e −m ϕ r ; the effective mass is defined as where φ 0 is the amplitude of the background value. The metric perturbations yield where the effective Newton constant G eff and the post-Newtonian parameter γ are defined as As it is clear from the above expressions, the coupling of the scalar field to the local system depends on φ 0 . If φ 0 1, then G eff ≈ G and γ ≈ 1 regardless of the value of m 2 ϕ . This contrasts with the result obtained in the metric version of f (R) theories, and, as long as φ 0 is sufficiently small, allows to pass the Solar System tests, even if the scalar field is very light.
In the modified cosmological dynamics, consider the spatially flat Friedman-Robertson-Walker (FRW) metric ds 2 = −dt 2 + a 2 (t)dx 2 , where a(t) is the scale factor. Thus, the modified Friedmann equations take the form respectively. The scalar field equation (100) becomes As a first approach, consider a model that arises by demanding that matter and curvature satisfy the same relation as in GR. Taking the trace equation automatically implies R = −κ 2 T + 2V 0 [29,213]. As T → 0 with the cosmic expansion, this model naturally evolves into a de Sitter phase, which requires V 0 ∼ Λ for consistency with observations. If V 1 is positive, the de Sitter regime represents the minimum of the potential. The effective mass for local experiments, m 2 ϕ = 2(V 0 − 2V 1 φ)/3, is then positive and small as long as φ < V 0 /V 1 . For sufficiently large V 1 one can make the field amplitude small enough to be in agreement with Solar System tests. It is interesting that the exact de Sitter solution is compatible with dynamics of the scalar field in this model.
Relative to the galactic dynamics, a generalized virial theorem, in the hybrid metric-Palatini gravity, was extensively analyzed [215]. More specifically, taking into account the relativistic collisionless Boltzmann equation, it was shown that the supplementary geometric terms in the gravitational field equations provide an effective contribution to the gravitational potential energy. The total virial mass is proportional to the effective mass associated with the new terms generated by the effective scalar field, and the baryonic mass. This shows that the geometric origin in the generalized virial theorem may account for the well-known virial mass discrepancy in clusters of galaxies. In addition to this, astrophysical applications of the model where explored, and it was shown that the model predicts that the effective mass associated to the scalar field, and its effects, extend beyond the virial radius of the clusters of galaxies. In the context of the galaxy cluster velocity dispersion profiles predicted by the hybrid metric-Palatini model, the generalized virial theorem can be an efficient tool in observationally testing the viability of this class of generalized gravity models. Thus, hybrid metric-Palatini gravity provides an effective alternative to the dark matter paradigm of present day cosmology and astrophysics.
In a monistic view of physics, one would expect Nature to make somehow a choice between the two distinct possibilities offered by metric and Palatini formalisms. We have shown, however, that a theory consistent with observations and combining elements of these two standards is possible. Hence gravity admits a diffuse formulation where mixed features of both formalisms allow one to successfully address large classes of phenomena.

Future Outlook: More General Hybrid Metric-Palatini Theories
The "hybrid" theory space is a priori large. In addition to the metric and its Levi-Civita connection, one has also an additional independent connection as a building block to construct curvature invariants from [216]. Thus one can consider various new terms such as Though an exhaustive analysis of such hybrid theories has not been performed, there is some evidence that the so called hybrid class of theories we are exclusively focusing our attention upon here is a unique class of viable higher order hybrid gravity theories. In the more restricted framework of purely metric theories, it is well known that the f (R) class of theories is exceptional by avoiding the otherwise generic Ostrogradski instabilities by allowing for a separation of the additional degrees of freedom in a harmless scalar degree of freedom: as we have already seen, such a separation is possible also for the hybrid theories. Furthermore, it turns out that this feature is a similar exception in the larger space of metric-affine theories, since a generic theory there is inhabited by ghosts, superluminalities or other unphysical degrees of freedom.

Slow-roll Inflation in f (R) Gravity
Cosmological inflation has been a great success over the last decades as it solves in an elegant way the problems raised in the Big Bang model: the so-called horizon and flatness problems, basically initial conditions problems. Inflation consists of initial super accelerated phase, during which the above-mentioned problems are solved by requiring a sufficient duration of such a phase (see Ref. [217][218][219] and references therein). After the end of inflation, the universe reheats and the original hot state of the Big Bang model is recovered. The original idea was proposed by Alan Guth in 1981, which was later improved by Andrei Linde in 1982. Over the years inflation has become a predictable theory, since it provides the seeds, coming from quantum fluctuations, that form the basis of the formation of large scale structures as well as the anisotropies observed in the CMB. Here, we review the essentials of the inflationary paradigm and we focus on the inflationary models within the context of extended theories of gravity, particularly the so-called f (R) gravity. Furthermore, the usual observables predicted by inflation are obtained and compared with the latest data released by the Plank mission [220].

Inflationary Parameters
Slow-roll inflation has been widely studied in the literature, particularly within the framework of scalar fields, since by the appropriate scalar potential, inflation can be easily reproduced and its exit achieved at the corresponding number of e-folds (see Ref. [221] and references therein). This is probably the simplest model for slow-roll inflation, whose general Lagrangian is given by a single minimally coupled scalar field, Then, by assuming a flat FLRW metric, the FLRW equations are obtained 3 Whereas the scalar field equation is given bÿ Inflation occurs basically when the scalar field behaves approximately as an effective cosmological constant, and the friction term in (110) dominates Hφ φ while V φ 2 . Inflation ends when the kinetic term of the scalar field increases and dominates over the potential, oscillating about the minimum of the potential and reheating the universe. A useful way of describing slow-roll inflation is to define the so-called slow-roll parameters: where the primes denote derivatives with respect to the scalar field φ. The magnitude of the quantities in (111) have to be small enough during inflation in order to attain the sufficient number of e-folds, or in other words 1 and η < 1. When inflation is ending, ≥ 1 is required. Thus, for the above model the spectral index n s of the curvature perturbations, the tensor-to-scalar ratio r of the density perturbations, and the running of the spectral index α s can be written in terms of the slow-roll parameters as follows The values of the above parameters can be compared then with the values provided by the latest Planck results (see [220]). Note that the aim here is to provide the above expressions (112) in terms of a particular underlying model, capable of reproducing slow-roll inflation, in order to compare with the observational data. For the simple model (108), the slow-roll parameters are already expressed in terms of the scalar potential in (111), while the reconstruction of the scalar field Lagrangian (108) given a particular inflationary model H(N), is provided by the FLRW equations (109), which yield [222] where ω(φ) is the kinetic term and we have redefined the scalar field as φ = N, where N is the number of e-folds N = log(a/a i ). Hence, by assuming a particular Hubble parameter, the corresponding Lagrangian is easily obtained. Note also that the slow-roll parameters (111) can be expressed in terms of the Hubble parameter by using (113), leading to [223][224][225] respectively. By assuming a particular Hubble parameter, the slow-roll conditions are straightforward calculated, while the corresponding Lagrangian for the scalar field is also obtained through the expressions (113). In the next section, we review previous works where instead of considering a scalar field, inflation is driven by modifications of the Hilbert-Einstein action.

Inflation in f (R) Gravity
Over the last years, f (R) gravity has been widely analyzed, since it can reproduce late-time acceleration as well as inflation with no need of an extra field and/or cosmological constant (for a review see [9]). Particularly, the reconstruction of the corresponding f (R) action that leads to the appropriate cosmological evolution has drawn much attention (see Ref. [226,227]). Here we review some aspects of f (R) gravities and specifically its roll within the inflationary paradigm.
Let us start with the general action for f (R) gravity, given in the action (62). The field equations are obtained by varying the above action with respect to the metric tensor, yielding Equation (63). Then, by assuming a spatially flat FLRW metric, the modified FLRW equations are obtained for f (R) gravity: As shown in Ref. [228,229], the corresponding slow-roll parameters (117) can be expressed in terms of the f (R) action as well. Nevertheless, here we are interested in reviewing slow-roll inflation in f (R) gravity by using the scalar-tensor equivalence of f (R) gravity and expressing the action (62) in the Einstein frame. Firstly, note that the action (62) in the Jordan frame, can be rewritten by action (67), so that by applying the conformal transformatioñ the action (67) is transformed as follows where Hence, the Lagrangian is now written in the form of (108), such that slow-roll inflation occurs when the conditions described above are satisfied. In addition, the slow-roll parameters (111) can be expressed in terms of the original f (R) action by using the correspondence expressions (66) and (121), leading to Then, given a particular form of the f (R) action, the corresponding slow-roll parameters can be calculated in terms of the Hubble parameter, which is obtained eventually by solving the FLRW equations. In addition, the corresponding spectral index and scalar-to-tensor ratio are obtained and for a particular number of e-folds, (most of inflationary models require N e ∼ 55 − 65), both magnitudes can be compared with the constraints from the Planck results [220].
Let us now consider the so-called Starobinsky inflation [230], which has been a very successful inflationary model, taking into account the observational constraints of the Planck satellite. The model has also been analysed from the quantum point of view, where it was shown to provide a correct behaviour and a successful inflationary period [231]. The model corresponds to a particular form of the f (R) function consisting of the Einstein-Hilbert term and a quadratic correction, given by The parameter m 2 in the Starobinsky model is a free parameter which is constrained by the observations, and it is in fact related to the mass of the new, dynamical scalar degree of freedom present in the model, usually denoted as the "scalaron". Although the inflationary dynamics of the model can be analysed either in the Jordan or in the Einstein frame, here we will focus our discussion on the Einstein-frame action, applying the procedure described earlier. Notice that according to the previous discussion, the requirement of the positivity of the second derivative of f , f RR > 0, translates into the requirement that the squared mass m 2 remains positive, i.e., m 2 > 0.
In order to re-express the model in the Einstein frame, we note that for the Starobinsky model the Einstein-frame scalar is given byφ Given the last relation, and using Equation (121), it is straightforward to derive the Einstein frame potential asṼ Therefore, the Einstein-frame action will be described by the action (120), with the potential defined through Equation (125).
During inflation, the background spacetime is described by a FRW spacetime, which after the conformal transformation (119) is expressed as follows whereã(t) denotes the scale factor in the Einstein frame. Similarly to Equation (109), the Friedmann equation in this spacetime is described bỹ while the equation of motion for the scalar φ takes the form whereṼ (φ) ≡ dṼ(φ)/dφ. The first and second term on the right hand side of the Friedmann equation (127) describe the kinetic and potential energy of the inflaton. Recall that slow-roll inflation occurs in the regime where κφ 1, i.e., as long as the fieldφ rolls down the flat part of the potential (125) sufficiently slowly, and the kinetic energy of the inflaton is much smaller than its potential energy,φ 2 Ṽ (φ). Then the Friedmann and the scalar field equation, Equations (127) and (128), respectively, become approximatelyH 2 κ 2 6Ṽ (φ), and 3Hφ −Ṽ (φ). The Hubble functionH is approximately constant during inflation, and its time variation is described by the slow-roll parameter defined as In the slow-roll approximation it is required that is small, 1. We can approximate it as The last step in the above relation assumed that κφ 1. Inflation ends as long as 1, which corresponds to the inflaton's kinetic energy becoming important.
The amount of inflation is measured by the number of e-foldings N which is defined as and for agreement with observations one usually requires that N 55 − 65. Now, under the slow-roll approximation, the above relation takes the following approximate form In the last step, we used the fact that in the slow-roll approximation, κφ 1, and we performed the integration with respect toφ. Notice that the number of e-foldings is related to the slow-roll parameters as 3 4 The scalar and tensor fluctuations produced during inflation have a characteristic amplitude and scale dependence, which in the Einstein frame, under the slow-roll approximation, are given by the standard relations for GR minimally coupled to a scalar field [217][218][219], while the tensor to scalar ratio follows from the above relations as r ≡ P GW /P R (48/π)N −2 16 . Notice that in the final approximate equalities of the above relations we restricted ourselves to the Starobinsky potential for κφ 1, under the slow-roll approximation. The spectral indices corresponding to the scalar and tensor fluctuations, respectively, under the slow-roll regime are given by the relations (112).
Note that Planck provides a value for the spectral index n s = 0.968 ± 0.006 and an upper bound for the scalar-to-tensor ratio r < 0.07 coming from a joint analysis of the BICEP2/Keck Array and Planck data [220]. In order to illustrate the power of Starobinsky inflation, let us assume a number of e-folds N e = 65, which leads to the following values of the inflationary observables: Hence, Starobinsky inflation satisfies quite well the constraints coming from Planck data.
Let us consider now a slight modification of Starobinsky inflation to show the uniqueness of the above model, where λ is a positive parameter that accounts for the deviation from Starobinsky inflation. In this case, and following the same steps as above, the scalar potential in the Einstein frame yields, Here A = λ 2κ 2 (6m) 2+λ and B = (6m) 2+λ 2+λ . This actually does not lead to an effective cosmological constant when κφ 1 since at that limit, the scalar potential becomes Since λ is positive, the potential is an exponential potential at large values of the scalar field which only mimics an effective cosmological constant in the case of λ = 0 when Starobinsky inflation is recovered. Otherwise, for κφ 1, the scalar potential is exponentially suppressed, and inflation does not occur. In the case of λ < 0, the scalar field will roll fast, invalidating the slow-roll conditions. Hence, deviations from λ = 0 will not reproduce slow-roll inflation, what may lead to other kinds of inflationary models as chaotic inflation which present deviations from the values of the spectral index and the scalar-to-tensor ratio as pointed by Planck Collaboration [220]).
Other extensions of f (R) gravity can also generate slow-roll inflation. In this sense, the inflationary paradigm has been recently explored within the so-called f (R, L φ ) gravities, where [232]). Such kind of theories extends f (R) gravities by including a nonminimally coupling scalar field. As shown in Ref. [232], this model is equivalent to multifield inflationary scenarios, where isocurvature modes may become important, while adiabatic perturbations can be analytically calculated (see [233]). However, following [232] the adiabatic perturbations can be obtained and the isocurvature modes may be ignored under certain conditions on the action.

Topological Defects
The Standard Model of Particle Physics describes three of the fundamental interactions between particles -the electromagnetic, weak and strong forces -in a unified quantum field theory. According to this model, the universe underwent, in its early stages, a series of symmetry breaking phase transitions, some of which may have led to the formation of topological defects.

Spontaneous Symmetry Breaking and Topological Defects
To illustrate the process of topological defect formation as a result of spontaneous symmetry breaking, let us consider the simplest model that admits topological defect solutions: the Goldstone model with a single scalar field, φ. In this model, the Lagrangian density is given by where λ is a coupling constant. Although this Lagrangian is invariant under transformations of the form φ → −φ -thus having a Z 2 symmetry -the Goldstone potential has two degenerate minima at φ = ±η. Hence, once φ acquires a vacuum expectation value (VEV), the Z 2 symmetry is spontaneously broken, leading to the formation of a kink separating regions where the scalar field has different VEVs (here, |0 is the ground state of the model).
In a 1 + 1 dimensional Minkowski spacetime, the equation of motion of the scalar field, admits the following solution with γ = (1 − v 2 ) −1/2 , that corresponds to a kink moving with velocity v in the positive x-direction.
Since φ interpolates between φ = −η at x = −∞ and φ = η at x = ∞, there is a region (the domain wall) centred at x = x 0 + vt which has non-vanishing energy. This configuration has a conserved topological charge, Q = dxj 0 = φ(+∞) − φ(−∞) = 2η (where j µ = µν φ ,ν is the topological current), and, for this reason, it is stable. The constant vacuum state has Q = 0, and, as a result, this configuration is unable to relax into the vacuum while conserving its topological charge. In N + 1 dimensions, the solution in Equation (143) remains valid and represents a (N − 1)-dimensional planar surface. These surfaces also separate two domains with different VEVs and, for this reason, they are commonly denominated Domain Walls. Domain walls arise due to the breaking of a discrete symmetry, whenever the vacuum manifold is disconnected. However, different types of defects may be formed depending on the (non-trivial) topology of the vacuum manifold (or the type of symmetry being broken). Line-like defects, or Cosmic strings, are formed if the vacuum is not simply connected or, equivalently, if it contains unshrinkable loops. This type of vacuum manifold results, in general, from the breaking of an axial symmetry. Moreover, if the vacuum manifold contains unshrinkable surfaces, the field might develop non-trivial configurations corresponding to point-like defects, known as Monopoles. The spontaneous symmetry breaking of more complex symmetry groups may lead to the formation of textures, delocalized topological defects which are unstable to collapse.

A VOS Model for Topological Defect Networks of Arbitrary Dimensionality
The production of topological defects is predicted in several grand unified scenarios. In most instances, these topological defects may survive throughout cosmological history until late cosmological times, and leave a variety of imprints on different observational probes. A unified phenomenological framework for the description of the cosmological evolution of topological defects of arbitrary dimensionality was developed in [234][235][236]. This model-which is a generalization of the Velocity-dependent One-Scale (VOS) model for cosmic strings [237]-provides a statistical description of the large-scale dynamics of these networks. This generalized VOS model characterizes quantitatively the macroscopic dynamics of the p-dimensional topological defect network in (N + 1)-dimensional Friedmann-Robertson-Walker backgrounds by describing the evolution of the characteristic length of the network, L, and of its root-mean-squared velocity (RMS),v.
The characteristic length of a statistically homogeneous defect network may be defined as whereρ is the average defect energy density, and σ p is the defect mass per unit p-dimensional area. For infinitely thin and featureless topological defects, the following equations describing the evolution of the characteristic length and velocity of the network (L andv, respectively) may be obtained directly from the generalized Nambu-Goto action [236] dv where D = N − p, and −1 d = (p + 1)H + −1 f is a damping length scale that includes, not only the effects of cosmological expansion, but also the effect of the frictional forces caused by the scattering of particles off the defects (which is encoded in the friction length scale, f ). Here, k andc are phenomenological parameters that account, respectively, for the effect of defect curvature and the energy loss associated to defect reconnection on the macroscopic dynamics of the network. The main advantage of this phenomenological model lies precisely on these two parameters. They may either be calibrated using numerical simulations-as was done for cosmic strings [238] and domain walls [239]-to accurately describe the large scale dynamics of standard networks, or be used in a phenomenological description of a large variety of alternative defect scenarios.
This simplest version of the VOS model has been extended in order to allow for the description of non-standard defect scenarios, with additional energy-momentum decay channels (e.g., direct energy loss through gravitational radiation [238]) or additional degrees of freedom (e.g., the existence of conserved currents [240]). Other extensions allow for the description of more complex defect scenarios-such as multi-tension cosmic (super-)string networks with junctions [241], semi-local cosmic string networks [242], or hybrid defect networks [243]-by including terms to account for the energy-momentum transfer between different network components.

Observational Signatures of Topological Defects
The observational consequences of topological defect networks can be very diverse, depending both on the type of defects formed and on the evolution of the universe after they are generated. For instance, the average domain wall density grows more rapidly than the background density in the radiation and matter eras and, consequently, domain walls must necessarily be very light in order not to completely dominate the energy density of the universe at late times. Although the possibility of a significant contribution to the dark energy budget has been ruled out both dynamically [36] and observationally [37], it is possible for light domain walls to leave behind interesting astrophysical and cosmological signatures.
Other topological defects, such as cosmic strings in a linear scaling regime, have an average density that is roughly proportional to the background density. Hence, standard cosmic string networks do not tend to dominate the energy density of the universe and they naturally generate a roughly scale-invariant spectrum of density perturbations on cosmological scales. Although recent cosmological data, in particular the cosmic microwave background anisotropies, rules out cosmic strings as the main source of large-scale cosmological perturbations (see, however, [244]), it remains plausible that they may play an important role on small scales. In fact, it has been shown that the power spectrum of a string-induced density perturbation component is expected to dominate over a primordial inflationary component on small enough scales [245]. In this context, cosmic strings have been suggested to have significant impact on the formation of ultracompact minihalos [39], globular clusters [40], super-massive black holes [41] and to provide a significant contribution to the reionization history of the Universe [42]. The dynamics of cosmic defects, in particular of domain walls, has also been associated to possible variations of the fundamental couplings of nature (see, e.g., [38]).
In this section, we shall focus on the contribution of both cosmic strings and domain walls to two of the most significant observational probes: the temperature and polarization anisotropies of the cosmic microwave background and the stochastic gravitational wave background.

CMB Anisotropies
The presence of a topological defect network at late cosmological times would necessarily leave imprints on the Cosmic Microwave Background (CMB). Although current CMB observations seem to be consistent with the inflationary paradigm (in which the fluctuations are seeded in the very early universe) they also allow for a subdominant topological defect contribution [246]. Topological defects source metric perturbations actively throughout the cosmological history. For this reason, their CMB signatures are expected to be fundamentally different from those due to primordial fluctuations. In particular, cosmic defects are expected to generate a significant vector component that would not be present in inflation-seeded scenarios, because vector modes decay rapidly in the absence of a source [247]. The B-mode polarization signal originated by topological defects has then contributions from both tensor and vector modes, and may, for this reason, produce an observationally relevant signal in this channel, despite providing only subdominant contributions to the temperature and E-mode power spectra. The B-mode polarization channel thus offers a relevant observational window for the detection of topological defects. With several CMB experiments planned for the near future-e.g., LITEBIRD, PRISM and CoRE+ (see [248] for a recent overview)-there is the prospect either for the detection of topological defects or for the tightening of current observational constraints in this scenario.
Cosmic string and domain wall networks produce temperature anisotropies characterized by significantly different angular power spectra. Domain walls only become cosmologically relevant at late times, when the characteristic length scale of the network is large. For this reason, these networks contribute mostly to the TT angular power spectrum at large scales (or small multipole modes ). Since low-observational data is highly affected by cosmic variance and has thus large associated uncertainties, current data allows for a significant contribution of standard domain wall networks to the temperature power spectrum. In Ref. [249], the authors have derived a conservative constraint on the domain wall mass per unit area of σ < 3.52 × 10 −5 kg m −2 allowed by current observational data (which corresponds to a constraint on the wall-forming symmetry breaking scale of η < 0.92 MeV). Cosmic strings, on the other hand, generate a TT angular power spectrum that has a broad peak at ∼ 300, and therefore contribute dominantly to the angular power spectrum at higher multipoles (which have small observational uncertainties). For this reason, observational data allow for a significantly larger contribution of domain walls to the temperature power spectrum on large angular scales than it does of cosmic strings. Current CMB observational data constrain the cosmic string tension to be Gµ < 1.3 (2.4) × 10 −7 for Nambu-Goto (Abelian-Higgs) strings [37] (which is significantly smaller than the corresponding constraint on domain wall mass per unit area GσL 0 = 5.6 × 10 −6 ). Theses differences in the shape and amplitude of the temperature power spectra generated by domain wall and cosmic string networks are clearly illustrated in the left panel of Figure 5, where the TT power spectra generated by cosmic string and domain wall networks with the maximum tensions allowed by current observational data are plotted. We chose GσL 0 = 5.6 × 10 −6 for the domain wall networks, and the weakest constraint on cosmic string tension obtained using Planck data [37], Gµ = 2.4 × 10 −7 (which corresponds to the constraint on Abelian-Higgs strings). Figure 5. The total temperature (left panel) and B-mode (right panel) angular power spectra generated by cosmic string and domain networks with the maximum tension allowed by current observational data: GσL 0 = 5.6 × 10 −6 and Gµ = 2.4 × 10 −7 for domain walls and cosmic string, respectively. The spectra were obtained using, respectively, the CMBACT code [250,251], and the extension introduced in [249]. These plots have been adapted from [249].
Although cosmic defect contributions to the temperature power spectrum are necessarily subdominant when compared to the primordial inflationary contribution, they may produce a significant B-mode polarization signal. As for primary temperature anisotropies, cosmic strings and domain walls contribute to the BB angular power spectrum at different scales: domain walls contribute dominantly at low-, while strings have a significant contribution at small scales. Moreover, the B-mode signal generated by domain wall networks has, in general, smaller amplitude than that of cosmic strings. In the right panel of Figure 5, we plot the total BB power spectrum generated by domain wall networks and cosmic string networks that have the maximum fractional contribution to the TT power spectrum allowed by current observational data. This figure shows that, since the observational constraints on Gµ are more stringent than those on GσL 0 , the domain wall contribution may dominate over that of cosmic strings for low multipole modes (particularly if their energy density is close to the maximum allowed by current data). Moreover, the fact that domain walls and cosmic strings have spectra with different shapes and that these defects contribute mostly at different scales should, in principle, allow one to distinguish between their contributions if a signal is detected.

The Stochastic Gravitational Wave Background
Cosmic string interactions often result in the formation of loops that detach from the long string network (for instance, in self-intersections or when kinky strings collide). After formation, these loops oscillate relativistically under the effect of their tension, and lose their energy in the form of gravitational waves. These loops, thus, have a finite life span and, as a consequence, the cosmic string network loses energy through this channel. The dynamical effects caused by this energy loss are taken into account in the VOS model through the inclusion of the last (phenomenological) term in Equation (146).
The VOS model may be used to estimate the rate of loop production (per comoving volume) as a result of cosmic string interactions dn c dt =c αv L 4 , where α is a parameter that characterizes the size of loops, l c = αL(t c ), at the moment of creation, t c . This rate of loop production is strongly dependent on the large-scale properties of the network (through its dependence onv and L) and, therefore, its accurate characterisation requires an accurate description of cosmic string network dynamics. The production of cosmic string loops is expected to occur copiously throughout the evolution of cosmic string networks. At any given time, there are several loops emitting gravitational waves in different directions. The superimposition of these emissions generates a stochastic gravitational wave background (SGWB) that spans a wide range of frequencies [252][253][254][255]. The SGWB is often quantified using the spectral density of gravitational waves per logarithmic frequency ( f ) interval, in units of critical density, where ρ GW is the gravitational radiation energy density. This quantity is strongly dependent on the rate of loop production as a function of time, and, consequently, on the large-scale cosmic string parameters. This spectrum has a characteristic shape: it has a flat portion in the high frequency range, generated by the emissions of loops created during the radiation era, and a pronounced peak at low frequencies that is produced during the matter era. Note however that the amplitude of the spectrum and the characteristics of the peak (its slope, broadness, location and height) vary greatly: they are not only dependent on the cosmic string tension and the large scale properties of the network, but also on the small-scale properties of the loops that generate it. For instance, the size of loops and their emission spectrum dramatically affects the overall shape of the spectrum. The dependence of the SGWB spectrum on these parameters is illustrated in Figure 6, where SGWB spectra generated by cosmic string networks are plotted, for a wide range of values of these macroscopic and microscopic parameters. There is the prospect of probing different parts of the SGWB spectra using current and upcoming gravitational wave detectors in the near future-most notably (Advanced) LIGO [256], (Advanced) VIRGO [257], eLISA [258] and pulsar timing arrays [259,260]. Understanding the shape of the spectrum generated by cosmic strings and its dependence on macroscopic and microscopic parameters is, therefore, crucial. Here, the loop size parameter is set to α = 10 −7 , and n s = 1. 4) For an increasing number of harmonic modes of emission, n s . For larger values of n s , the shape of the spectrum saturates and is identical to that of the spectrum with n s = 10 4 . Here, α = 10 −8 , Gµ = 10 −7 , and we have only considered loops with cusps. These plots were adapted from [255].
As in the case of cosmic strings, as domain wall networks evolve and interact, a fraction of their energy is released in the form of gravitational waves. In fact, in [261,262], it has been shown using field theory simulations that domain wall networks produce a spectrum of gravitational waves that peaks at the scale corresponding to the size of the Hubble radius at the time of their decay. Studies of the gravitational wave spectrum generated by domain wall networks are, however, preliminary and limited to the case of metastable networks.

Cosmological Tests with Galaxy Clusters at CMB Frequencies
The Cosmic Microwave Background (CMB) has revolutionized the way we perceive the Universe. The information encoded in the temperature and polarization of the CMB provides one of the strongest evidences in favour of the hot Big-Bang theory and has enabled ways to constrain cosmological models with unprecedented accuracy [44]. The CMB also provides a unique way to study the formation and evolution of cosmic structure. Physical effects acting on the CMB after decoupling originate secondary CMB anisotropies that encode additional information about the growth of cosmological structure and the dynamics of the Universe [45]. The Sunyaev-Zeldovich (SZ) effect [46,47] is the scatter of CMB photons by plasma in galaxy clusters and filamentary structures. It is among the strongest and most studied sources of secondary anisotropies. Whereas the primordial CMB signal has been widely used to constrain cosmological parameters, the SZ effect is only now giving the first solid steps into adulthood as a probe of Cosmology and structure formation. Largely contributing to this state of affairs are ground-based and satellite observations such as those by SPT [263] and Planck [264] and proposals for future experiments such as in [265].
In this section we review some of the key features of the SZ effect and its application to test cosmological models and the physics of galaxy clusters, in light of Planck cluster observations. We end with an example of the application of cluster data to the study alternative cosmological scenarios involving primordial non-Gaussianities.

Galaxy Clusters
Galaxy clusters are the largest gravitationally bound objects in the Universe. They are formed by gravitational accretion of infalling material at the intersection of filamentary structures. Their typical masses range from 10 14 − 10 15 h −1 M and contain hundreds to several thousand galaxies within radii of the order of 1 − 1.5 h −1 Mpc. The intra-cluster medium (ICM) is filled with hot ionized gas, typically at temperatures 1-15 keV, which produces strong X-ray emission and microwave spectral distortions via the SZ effect, see Figure 7. Because clusters are among the latest bound structures forming in the Universe their number density is highly sensitive to the growth rate and power spectrum of density perturbations as well as to cosmological background parameters like the Universe's matter density and the Hubble parameter. The central quantity that goes into the computation of the expected abundance of clusters for a given cosmological scenario is the mass function, n(M, z). It gives the comoving number density of virialized halos forming at the redshift z with mass M. As introduced by [266], this can be expressed as: where f (σ, z) is the so called multiplicity function,ρ is the mean background density and σ(M, z) is the variance of the linear density field at redshift z, smoothed on a scale R, enclosing the mass M, Here, D(z) is the growth factor of linear perturbations, P(k) is the linear density power spectrum extrapolated to redshift zero, and W(k, M) is the Fourier transform of the smoothing kernel in real space. The multiplicity function can be computed analytically or by assessing the number of halos forming in numerical N-body simulations. The following formulae give the multiplicity function for the popular Press-Schecter, f PS , [268] and Jenkins, f J , [266] mass functions: Here δ c is the extrapolated linear overdensity of a spherical perturbation at the time of collapse. The Press-Schecter mass function was the first being derived using analytical arguments, but it overpredicts/underpredicts the abundance of low/high mass halos in simulations. The Jenkins mass function was established using simulations. Several other simulation based mass functions have been proposed to improve the functional form of the multiplicity function, see e.g., [269] for a summary with other multiplicity function expressions. Most of these studies have been focusing on improving systematic uncertainties related to the halo mass definition and quantifying deviations from the universality of f for different cosmologies. This includes extensive mass function studies using primordial non-Gaussian N-body simulations, see e.g., [270], and f (R) modified gravity simulations, see e.g., [271].
The redshift distribution of observed galaxy cluster abundances can be compared against model predictions using the expression: Here d 2 N/dzdΩ is the number of clusters per unit redshift per unit solid angle, d 2 V/dzdΩ is the volume element of the underlying model, and f survey (M, z) is the cluster survey selection function. Setting f survey = 1 and assuming that its possible to observe all cluster masses, Equation (152) gives the full redshift distribution of halos in the Universe, see e.g., [272] for an application to dark energy cluster counts. The f survey function depends on survey parameters such as sky coverage and noise distribution of the observed quantity S ν , which is usually a cluster flux or band luminosity. In a best case scenario, where the sky noise is uniform and the survey is complete above a given flux limit, S lim , Equation (152) can be computed as an integral over the observed flux, (dn/dM) (dM/dS ν ) dS ν , with S ν ≥ S ν,lim . Here dM/dS ν is the derivative of a scaling relation between M and S ν . Galaxy cluster number counts are therefore a powerful probe of cosmology provided the relation between total mass and the observed survey flux are known. Poor knowledge about the later may imply large uncertainties in the determination of cosmological parameters, as well as leaving a number of unanswered questions regarding the underlying cluster formation mechanism.
If gravity is the dominating force driving the formation and evolution of collapsed structures, virialized clusters should have their thermodynamic properties related to the total mass (through the viral theorem). The simplest model describing cluster scaling relations was presented by [273]. This assumes, among other things, that gravitational collapse is scale-free, i.e., clusters of different masses are self-similar replicas of each-other. If non-gravitational effects, such as radiative gas physics and energy feedback, influence the thermodynamic state of clusters one should expect modified scalings. These are usually parameterized as power laws of the form [269], where Y, X are cluster properties, A is a normalization parameter giving the amplitude of Y at X = X 0 , E(z) = H(z)/H 0 , α is the power law index of the (independent) property X, and β ss is a parameter giving the self-similar evolution of the normalization. In general, the amplitude parameter A may itself be a function of redshift. In this case, A evolves in a non self-similar way, usually parameterized by a power-law of redshift, A(z) = A 0 (1 + z) β . The parameters α and β ss can be computed analytically within an extended self-similar modelling that also assumes that the baryon gas component in clusters is in hydrostatic equilibrium [269]. For example, under these conditions, the scaling between the SZ integrated flux and mass, Y − M, is given by Equation (153) with the indexes α = 5/3, and β ss = 2/3. Hydrodynamic N-body simulations are the appropriate method to study galaxy cluster scaling relations as they allow us to model the baryon physics throughout the full non-linear evolution of the ICM gas and determine the parameters α, β ss , β and A 0 . We refer the reader to references [274][275][276][277] for examples of simulations studies showing the impact of various gas physical models, dark energy and primordial non-Gaussianities on several galaxy cluster scalings. Numerical N-body simulations have also been intensively used to study the internal structure of galaxy clusters. These studies demonstrate that galaxy clusters profiles can be used as a sensitive probe of cosmology, cluster physics and structure formation mechanisms. One of the most well known results from simulations, in the context of hierarchical CDM models, is that the radial distribution of dark (collissionless) matter in clusters follows a simple expression known as the Navarro-Frenk-White (NFW) density profile [278] where x ≡ r/r s , r s is a scale radius, and ρ s is the characteristic density at r = r s . Subsequent numerical simulations allowed one to improve this result and to study the impact on cluster profiles of alternative cosmological models such as modified gravity, see [269] for a review. Gas profiles have also been extensively investigated with hydrodynamic N-body simulations. For example, [279] found that the radial pressure profile, P(r), in SZ observations of resolved clusters follows a universal profile, P(r) P 500 where x = r/r s , r s is the scale radius, and P 0 , α, β, γ are fitting parameters to derive from observations. P 500 is the integrated pressure at r 500 (the radius at which the cluster density exceeds 500 times the universe mean background density). Equation (155) is known as the generalized NFW (GNFW) profile of the SZ signal in clusters.

The SZ Effect: Temperature and Polarization Signal
The Sunyaev-Zel'dovich effect is the scattering of CMB photons by electrons in hot reservoirs of ionized gas in the Universe. A particularly effective scattering medium is the ICM where the gas is at temperatures of the order of 10 8 Kelvin. When the much colder and isotropic CMB radiation propagates through these plasma regions, a fraction of the original CMB photons scatter off the moving electrons and suffer frequency Doppler shifts towards higher energies (inverse Compton scattering) causing a spectral distortion in the CMB spectrum. In the process the total number of photons is conserved. In fact, the SZ effect can be regarded as a composition of two effects: the thermal and kinetic SZ effects. The thermal effect is due to the thermal (random) motion of the electrons in the gas cloud. If the gas has also a bulk peculiar velocity in the CMB radiation frame, there is a kinematic Doppler effect in addition to the thermal effect. The spectral distortions caused by the thermal, ∆I th , and kinetic, ∆I k , SZ effects are given by simple formulae [47]: where I 0 = 2k 3 B T 3 /h 2 c 2 , T is the CMB temperature, y is the line-of-sight Compton SZ parameter, τ is the optical depth, v r is the line-of-sight velocity of the center of mass of the gas cloud (positive/negative when the cloud is receding/approaching), T e is the electron density, n e is the electron number density, and g(x), h(x) are functions giving the frequency dependence of the effects. Figure 8 shows maps of the thermal and kinetic SZ distortions imprinted on the CMB by cosmic structures in hydrodynamic N-body simulations. The mapped quantities are the y parameter and the temperature fluctuations of the kinetic SZ effect, (∆T/T) k = ∆I k /(I 0 h(x)) = −v r /c τ, see e.g., [280]. For typical cluster velocities and optical depths, the thermal SZ effect is the dominant effect. Therefore the total SZ signal integrated over the sky angular size of a cluster is well approximate by where d A is the angular diameter distance from the observer to the cluster. This is known as the volume integrated SZ Y-flux or Y-luminosity and is a measure of the total thermal energy density of the cluster. The scattering of CMB photons by hot gas in clusters also generates new (secondary) polarization signatures in the CMB that carry information about the intra-cluster medium. The dominant CMB secondary polarization effect in clusters is due to the existence of an intrinsic quadrupole component from the CMB temperature anisotropy around the scattering media. This effect is known as the CMB quadrupole induced polarization. In addition, the bulk motion of the gas cloud may also induce CMB polarization proportional to the square of the sky transverse velocity of the cloud. This effect is often referred to as induced kinetic polarization. Finally, double scattering events inside the cluster gas cloud also originate polarization, known as double scattering induced polarization. The bottom panels in Figure 7 are detailed simulations of the linear polarization degree and polarization angle of the CMB quadrupole and kinetic induced polarizations in a cluster. These simulations are an essential tool to assess the degree of contamination of the primordial CMB signal by secondary polarization effects and allow to prepare observational strategies to use CMB cluster polarization as an additional probe of cosmology and structure formation, see [267] for a detailed description of cluster polarization simulations and their use as a new observing tool.

Cosmological Tests Using SZ Cluster Surveys
SZ galaxy cluster surveys are recognized as having a great potential to test a wide range of cosmological and astrophysical scenarios. SZ galaxy cluster counts, profiles, scaling relations and SZ angular power spectra are most promising probes to confront model predictions with observations. In this section we will briefly review the current state of cosmological constraints from observations with the Planck satellite.

SZ Cluster Counts
The Planck SZ cluster survey is the largest and deepest all-sky cluster survey yet produced at CMB frequencies. Its latest galaxy cluster catalogue (PSZ2, [281]) contains 1653 detections, of which 1203 are confirmed clusters with identified counterparts in other observing bands (X-rays, optical, near infra-red and other SZ experiments). More than a thousand of these clusters possess redshift information. From this catalogue, the authors in [264] selected a sample of 439 clusters with high signal to noise (q = S/N > 6) to infer constraints on cosmological parameters using SZ galaxy cluster number counts.
For a given cosmological model, the redshift distribution of clusters abundances detected by Planck above a given noise level, q > q cat can be computed from Equation (152) as: see [264]. Here M is the mass contained within a cluster radius, R 500 , where the cluster density exceeds by 500 times the mean background density of the Universe,χ(M, z, l, b, q cat ) = ∞ q cat P[q|q m (M, z, l, b)] dq is the Planck selection function, (i.e.,χ = f survey in Equation (152)) and P[q|q m (M, z, l, b)] is the distribution of q given the mean signal-to-noise value,q m (M, z, l, b), predicted by the model for a cluster of mass M and redshift z located at Galactic coordinates (l, b). This expression reflects the fact that the survey noise level is not uniform on the sky and that clusters are detected above a given noise level, q cat . The distribution P[q|q m ] incorporates noise fluctuations and intrinsic scatter in the actual cluster SZ Y signal, (158), computed within R 500 , around the mean value,Ȳ(M, z), predicted from a Y − M scaling relation calibrated using Planck cluster data.
This method was used in [264] to impose constraints on cosmological and cluster parameters, such as the variance of the density fluctuations smoothed on a scale of R = 8 h −1 Mpc, σ 8 , see Equation (150), the matter density parameter, Ω m0 , the Hubble constant, H 0 , and the slope of the Y − M scaling relation, α. Assuming a ΛCDM concordance (base) model their analysis show that Planck SZ cluster constraints on the σ 8 − Ω m0 plane are in tension with those obtained from Planck primary CMB anisotropies, unless priors on the hydrostatic mass bias parameter, 1 − b, (which relates the observed baryon mass with the total mass of the cluster) and the effect of gas physics on the growth of structure are changed to values that are in conflict with constraints from other cosmological probes and the predictions from numerical simulations, [264,282]. A way to try to alleviate the tension is to combine cluster and primary CMB data to explore constraints on extensions to the base flat CDM model. This was carried out in [264] by allowing an extra varying parameter (to the base model), namely, the neutrino mass scale, the Thomson optical depth to reionization, the dark energy equation of state, and the Universe curvature parameter. The analysis seems to favor non-minimal neutrino masses (see Figures 11 and 14 in [264]), but this "solution" also implies a decrease on the Hubble constant to values that further deviate from most constraints by astrophysical determinations of H 0 .
The present Planck analysis also puts in evidence the importance of improving the precision of the cluster mass bias parameter, which in fact may depend on cluster mass and redshift. In the future, large lensing surveys, such as Euclid, WFIRST and LSST, are expected to determine mass biases with one percent level accuracy. With this level of precision on 1 − b, the present cluster combined analysis of Planck data could by itself provide a stringent test of the ΛCDM concordance model, see Figure 13 in [264].

SZ Power Spectra
An important statistic to derive from sky maps of the SZ effect is the angular power spectrum of the thermal and kinetic SZ temperature fluctuations. The theoretical computation of these spectra is straightforward from simulations, such as those in Figure 8, but the computational resources involved to run simulations limit their use if one is required to compute the SZ power spectra for many model realizations. Analytical methods are much faster and allow to explore many model realizations, but they are much less capable of modelling the cluster gas physics and the contribution from the diffuse gas in filamentary structures.
Analytically, the angular power spectrum of the thermal SZ signal can be computed as a sum of two terms C total = C Poisson + C clustering , where C Poisson is the power spectrum resulting from a Poisson distribution of objects (clusters) and C clustering is a "2-halo" or clustering term giving the power arising from object (cluster two-point) correlations. These terms are given by (see [45] and references therein): where y is the 2D Fourier transform of the projected 3D Compton y-parameter radial profile of individual cluster, P(k, z) is the 3D matter power spectrum at redshift z, and b(M, z) is the linear bias factor that relates the matter power spectrum, P(k, z), to the cluster correlation power spectrum P cluster (k, M 1 , M 2 , z) = b(M 1 , z)b(M 2 , z)D 2 (z)P(k, z = 0). The mass and redshift range of the integrations are chosen so as to cover all relevant SZ source contributions up to high enough redshift.
Observationally, the SZ all sky anisotropies can be separated from the overall CMB signal by applying component separation techniques that make use of the spectral signature and scale of the SZ effect. The first all-sky maps of the thermal SZ effect were produced from Planck observations [283,284]. These studies allowed to investigate the distribution of the y signal and characterize its higher order statistics, including the power spectrum and bi-spectrum. The observed thermal SZ power spectrum confirmed to be very sensitive to cosmological parameters, especially σ 8 and Ω m0 , but these are strongly degenerate with the mass bias parameter b (they found that σ 8 (Ω m0 /0.28 These results are consistent with the constraints on σ 8 (Ω m0 /0.28) 3/8 from cluster number counts but are again in "soft tension" with the constraints obtained from primary CMB data assuming a ΛCDM concordance base model, see Figure  17 in [284].

Probing New Physics with the SZ Clusters
The present best CMB spectral observations by COBE/FIRAS are unable to determine if the CMB spectrum preserved its planckian shape at all redshifts and if the CMB temperature-redshift scaling is the one predicted by the standard Big-Bang model, T = T 0 (1 + z), where T 0 is the CMB temperature observed at z = 0. However it is possible to consider alternative models where the CMB has a planckian spectrum at z = 0, but a different temperature-redshift scaling. Observational tests of this scaling require measurements of the CMB temperature at high z, i.e., throughout the Universe's history. The SZ effect has long been recognized as a most promising method to measure the CMB temperature at the redshift and location of clusters, allowing to test the CMB temperature-redshift relation at high redshifts, see e.g., [285] and references therein. Recent observations by Planck allowed to impose constraints on a simple deviation of the standard scaling, T = T 0 (1 + z) 1−β , using CMB primary data [44] and SZ effect cluster data, see e.g., [286]. The emerging view is that Planck data impose tight constraints on the T(z) that are consistent with the standard CMB temperature scaling predicted by the standard Big-Bang model (i.e., β = 0).

Probing Primordial Non-Gaussianities with Galaxy Clusters
The formation of the structures observed at large scales in the Universes is one of the most challenging and complex problems still open in modern cosmology. The widespread understanding is that these structures are the result of tiny primordial density fluctuations that might have existed at early times and that eventually grew under gravitational instability collapse. The origin of these tiny primordial density density perturbations is a cause of great debate, yet the inflationary paradigm remains the consensus framework for their generation.
The statistical properties of the primordial density perturbations predicted by most standard single-field-slow roll inflationary models are nearly Gaussian in nature, and thus the level of non-Gaussianity is not strong enough to be observationally detected [55][56][57]. Although recent CMB anisotropies and LSS data seem to support such prediction, more complex models may lead to the generation of larger a non-Gaussianity amplitude, which substantially increases the chances of detection with future high precision observational experiments. For example, under certain conditions, some multi-field inflation models where slow-roll holds, can lead to larger departure from Gaussianity. However, non-Gaussianities may arise in inflationary models in which any of the conditions leading to the slow-roll dynamics fail [58]. The spectrum of models capable of generating larger non-Gaussianities amplitudes include the curvaton scenario [59][60][61], the ekpyrotic inflacionary scenario [62,63], vector field populated inflation [64][65][66] and multi-field inflation [67][68][69][70]. The current efforts that are being done to detect and constrain primordial non-Gaussianities (hereafter PNG) using a variety of cosmological observables is of the tremendous importance. A positive detection of non-Gaussianity signal would allow to asses the different proposed mechanisms for the generation of primordial density perturbations and in particular, it would make possible to discriminate between different inflation models and rule out some of them. Moreover, our understanding of the key physical processes at the early Universe would benefit greatly from such outcome.
There is a wide range of observables to which we can resort to, in order to search for primordial non-Gaussianities. The traditional way to proceed is to measure three-point statistics of the CMB temperature anisotropies maps. Using this observable, there have been in the past reports of a positive detection of non-Gaussianity at a significant level [287]. However, although some contradicting analysis exist (see [288,289] and more recently [290,291]), the interest in testing deviations from Gaussian initial conditions has gained momentum ever since. Nevertheless, the possible existence of primordial non-Gaussianities could also be tested using other cosmological probes such as the characterization of the statistical properties of the structures at large scales. These probes include the the bispectrum and/trispectrum of distributions of galaxies and galaxy clusters (e.g., [292][293][294]), weak-lensing observations, (e.g., [295,296]) and CMB-LSS [297] as well as CMB-21cm line [298] cross-correlations. Another set of cosmological of probes very discussed in the literature to contain PNG, is the time evolution of the numerical abundances both of massive collapsed objects such as galaxy clusters, and its counterparts, large voids (see e.g [299,300] and [301][302][303][304]).

Parametrizing Primordial Non-Gaussianity
The cosmological information contained in the power spectrum (or equivalently the two-point correlation function) of the primordial density perturbations is not enough if one wants to look for deviations from Gaussian initial conditions. Therefore we must resort to higher order correlation functions in order to gain further insight on nature of non-Gaussianities. Thus, the bispectrum, B Φ is the lowest order statistics sensible to primordial non-Gaussianities and it is defined as where Φ is the Bardeen's potential and δ (3) D is the three dimensional Dirac delta. The bispectrum can be be written in the alternative form [305], where the dimensionless parameter f NL accounts for the amplitude and the shape function A stores the information on the shape of the bispectrum. Due to the large number of inflationary models capable of generating primordial non-Gaussianities, it is common practice to group them according to their resulting bispectrum shape. There are generically four classes of bispectrum shapes frequently studied in the literature: Local, Equilateral, Folded and Orthogonal. However, from these shapes, the most studied in the literature is the Local one, which can be written in terms of Bardeen's potential and performing a Taylor expansion around an auxiliary Gaussian random field Φ, (see eg. [57]) The corresponding bispectrum of the Local shapes is given by, To date, the strongest constraints on primordial non-Gaussianities of the Local type were obtained by the Planck collaboration [291] (see also [290] for previous constraint PNG by Planck collaboration), f local NL = 2.5 ± 5.7 from temperature alone and f local NL = 0.8 ± 5 combining temperature and polarization data. Note however, that although these constraints are consistent with Gaussian initial conditions, they were obtained assuming that the non-linear parameter f NL is scale-independent. Nevertheless, some inflationary models may lead to a scale dependent f NL (see [306] and references therein).

The Non-Gaussian Mass Function
There is no unique prescription to incorporate the effects induces by primordial non-Gaussianities in the mass function (see [300,[307][308][309][310][311][312]). The majority of the of the non-Gaussian mass functions are obtained resorting to some expansion of the probability density distribution (PDF) of the primordial density field and either using an extension of the Press-Schecter formalism [268] or the excursion set theory [313,314]. Generically, the modified mass function can be written in the form where R NG encodes the correction to the Gaussian mass function due to PNG. In (166), the Gaussian mass function can be simply the one of the PS formalism or any other motivated by N-body numerical simulation with gaussian initial conditions (references of fitting formulas), such as Equation (151). Two of the most discussed and used expression for the R NG correction in the literature are the MVJ [300] and the LV [306]. Both rely on the extension of the PS formalism to non-Gaussian initial conditions. The respective functional forms are, where δ c (z) is the collapse threshold at the redshift z, S 3 is the reduced skewness and δ * (z) = δ c (z) 1 − S 3 δ c (z) /3. Figure 9 show the dependence of the non-Gaussian corrective term to the mass function, R NG , the mass and as a function of redshift, for the Equations (167) (blue line) and (168) (red dashed line).
In both cases f NL = ±100 was used. It is clear that a higher positive/negative f NL value increases the number of predicted dark matter halos with higher mass. The effects are stronger at high redshift. The MJV prescription seems to predict a higher number of halos than LV for positive values of f NL . The opposite occurs when negative f NL values are considered.  (167)), while red dashed line corresponds to LV formula (Equation (168)). The level of non-Gaussianity of the Local type used was f NL ± 100.

Biased Cosmological Parameter Estimation with Cluster Counts
The statistical properties of the primodial density perturbations has a great deal of influence on the redshift distribution of the number of galaxy clusters. Thus, these objects can be used as cosmological probes to search for deviations from Gaussian initial conditions. However, neglecting erroneously the existence of primordial non-Gaussianities, can lead to biases in the estimation of the cosmological parameters using clusters number counts. In [315], the authors investigated how primordial non-Gaussianities influence the estimations of the effective dark energy equation of state parameter, when the information contained in the abundance of galaxy clusters is used to explore different cosmological scenarios. In their work, they computed the effective dark energy equation of state per redshift bin, assuming Gaussian initial conditions, that would allow them to recover the galaxy clusters counts in different non-Gaussian models. From their findings, there resulted a new diagnostic feature for the presence of PNG, in the form of an apparent evolution with time of the effective dark energy equation of state, characterized by the appearance of a discontinuity (more details see [315]). The same authors also estimated the magnitude of the biases that may arise in the determination of a larger set of cosmological parameters with galaxy clusters number counts, again ignoring that PNG may exist (see [316]). Their results indicate that, although the estimation of the present-day dark energy density and respective equation of state parameter do not seem to be sensible to non-Gaussian initial conditions in the primordial matter density field, the biases can be much larger in the case of the amplitude of the primordial density perturbations. Furthermore, the authors also argue that the conflicting constrains on the amplitude of the primordial perturbations obtained by the Planck collaboration using galaxy cluster number counts from the Planck Sunyaev-Zel'dovich Catalog and the primary Cosmic Microwave Background temperature anisotropies, could be alleviated, if a significant non-Gaussianity level at clusters scales exists (see [316] for more details).

The Impact of pRimordial Non-Gaussianities on Galaxy Clusters Scaling Relations
In order to use galaxy clusters as cosmological probes, it is essential to have an accurate estimate of their mass. However, this task is not easy, since it is not possible to directly estimate the mass of a galaxy cluster from observational data, and thus one must rely on indirect methods, such as the scaling properties of clusters, to obtain a reliable estimate. The effect of PNG on galaxy clusters scaling relations was investigated for the first time in [277]. The authors resorted to a series of hydrodynamical N-body numerical simulations featuring adiabatic gas physics and allowing for different levels of non-Gaussian initial conditions within the fiducial ΛCDM model. In particular they studied the T − M, S − M, Y − M and Y X − M scalings relating the total cluster mass with temperature, entropy and SZ cluster integrated pressure that reflect the thermodynamical state of the intra-cluster medium. Figure 10 shows as an example, the distribution of galaxy clusters for the Y − M, for three different levels of PNG, i.e., f NG ± 500, 0 obtained by [277], while Figure 11 shows the result of fitting Equation (153) for the Y − M scaling to the corresponding distribution of clusters, using the procedure suggested by [276,317]. The main findings of [277], indicate that PNG do have an impact on clusters scaling laws, in particular on the amplitude and redshift evolution of the scalings normalization. In the case of Y − M, the redshift evolution of the its normalization can change as much as 22% when f NL varies from −500 to 500. If f NL is a scale independent parameter and has a small value, as suggested by the current Planck constraints [290,291], then the impact of PNG on clusters scaling laws is negligible when compared with the effect of additional gas physics and other cosmological effects, such as dark energy. On the other hand, if f NL is instead a scale dependent parameter, as suggested by some inflationary models, then PNG could have larger positive/negative amplitude at clusters scales. In this case, the effects of PNG must be taken into account when galaxy cluster data are used to assess the constraining power of future cluster surveys or to estimate the cosmological parameters from observational cluster counts data.

Testing Gravity with Weak Lensing
Gravitational lensing is the name commonly given to the phenomenon of deflection of light by gravitational fields. A lensing system consists of three parts: the source of light, the lens (a inhomogeneous distribution of mass, or a spacetime geometry) and the observer. Gravitational lensing of a point source by intervening matter produces image displacement, or even multiple images of the source. On extended sources, tidal gravitational fields lead to differential deflection, producing shape distortion and size magnification. Gravitational lensing on cosmological scales is a direct probe of Poisson-like constraint equations and probes structure formation independently of mass-to-light properties. It is thus a promising powerful way of testing gravity on cosmological scales and it is a core probe in the Euclid space mission (http://www.euclid-ec.org).

Gravitational Lensing
Gravitational lensing [318][319][320][321] involve different types of astrophysical objects that spread in a broad range of mass, mass density and length scale. Various types of systems have been studied, both theoretically and observationally, with the aim of investigating the source, the lens, or the geometrical properties of the spacetime, from the observed images. Sources may range from stars, galaxies, to the last scattering surface, while lenses may be stars, galaxies, galaxy clusters or the large scale structure. The lensing properties are the same in all cases: the effect depends only on the projected two-dimensional mass density distribution of the lens and is independent of the luminosity or composition of the lens; gravitational lenses do not have a well-defined focal point; the effect is achromatic; there is no emission or absorption of photons; the surface brightness is conserved. The conservation of brightness implies that a change in size changes the flux, as if the source was closer to the observer, creating thus a zooming effect without loss of resolution, the so-called natural telescope.
The lens equation is a mapping between source (β) and image (θ) planes, The deflection field α(θ) is a gradient of the gravitational potential and it is a central quantity of gravitational lensing. It is determined by the lensing potential, i.e., the gravitational potential integrated along the line-of-sight, and it depends thus on the geometry of the system (or the spacetime metric in cosmological-size systems). It is useful to Taylor-expand the lens equation, a procedure valid for small sources in angular size and small deflections (implying small angular size images). To linear order, the lens equation will thus be described by a Jacobian mapping, known as the amplification matrix: When applied to a light ray bundle emitted from an extended source, the symmetrical traceless part of the matrix produces a sheared image, the trace is responsible for an isotropically convergence (or divergence) and the antisymmetrical component introduces a rotation. This decomposition defines the optical components fields in the image plane: the scalar convergence field κ( θ) and the shear field γ(θ). The convergence is given by the Laplacian of the lensing potential, and it corresponds thus, through the Poisson equation, to the surface mass density in the effective two-dimensional plane where ψ is defined. The shear is a spin-2 quantity and is given by, The rotational field is not present, since the deflection field is a gradient and the amplification matrix is thus curl-free.
Lensing effects are classified as strong or weak. Strong lensing effects occur at positions of the image plane where the convergence and shear fields take large values (> 1). In these cases multiple images may be formed, as well as strongly deformed multiple images that may merge forming giant arcs or even a full ring, the Einstein ring. The magnification of the images µ(θ) is given by the inverse of the determinant of the amplification matrix. Images are magnified or demagnified depending on µ being greater or smaller than 1. The sign of the determinant defines the parity of the images, images in regions of det A < 0 are inverted with respect to the source. The points in the image plane where det A = 0 form closed curves named critical curves. The observed giant arcs form near critical curves. Strong lensing occurs at various scales: a quasar (source) -galaxy (lens) system produces multiple images of the quasar with a typical separation of arcsec, while a background star -foreground star system produces multiple images with a typical separation of miliarcsec, which is known as microlensing. If κ 1 everywhere, multiple images do not form. An example in this intermediate regime are the arclets formed in the images of galaxy clusters, by lensing of background galaxies. The separation of the lensing regimes in an image depends mainly on the position with respect to the source-lens-observer alignment, as illustrated in Figure 12 (see also [322]). In the linear regime (weak lensing), when κ 1, the distortion on one source is not noticeable. However, the lensing effect can still be detected statistically when there is a large enough sample of sources suffering a coherent lensing effect. This is the case with images of dense populations of faint galaxies present in deep optical and near-infrared observations. The images of these background galaxies are coherently sheared by a large-scale inhomogeneous matter distribution: the large scale structure of the universe. The shear field of such images is known as cosmic shear and shows angular correlations up to degree scales for sources at z ∼ 1. Cosmic shear correlations are a direct measure of the statistical properties of the lensing potential, which depend on the statistical properties of the cosmological density field and the geometrical properties of spacetime. For this reason cosmic shear is a probe of both geometry and structure formation and it is the lensing effect we deal with in this contribution.

Cosmic Shear
Observationally, a pioneering attempt to measure a cosmological weak lensing signal was made in [323], having produced an upper limit for the signal. A decade later, [324] measured shear correlations with error bars of order 100% and verified the existence of instrumental and atmospheric anisotropical contaminations. A first significant detection, cosmic variance dominated, was made in [325] on a very small field (4 sqarcmin). Finally in 2000, the first definitive detection of the cosmological shear field signal was reported almost simultaneously by four independent teams [326][327][328][329]. A large number of surveys from various ground-based observatories was soon undertaken, confirming the detections and using two-point cosmic shear correlations to constrain Ω m and σ 8 parameters (see [330] for a compilation). Those surveys culminated on the Canada-France-Hawaii Telescope Legacy Survey (CFHTLS), a dedicated cosmic shear survey covering 170 deg 2 that successfully detected a signal on degree scales [331]. In the past 5 years the focus of cosmic shear observations moved from the detection itself to the robust estimate of residual systematics. Currently, the state-of-the-art cosmic shear survey is the Canada-France-Hawaii Telescope Lensing Survey (CFHTLenS) [332], a reanalysis of CFHTLS where new methods have been applied to measure the galaxy distortions [333], verify the source redshift distributions [334] and remove intrinsic alignment correlations [335]. These robust data were used in various cosmological applications, including testing GR on cosmological scales, both alone and in combination with other probes, both for specific gravity theories and parametrized deviations to GR [336][337][338].
Cosmic shear studies deal with the light from background galaxies that propagates in the perturbed FRW spacetime [48,339,340]. In the Newtonian gauge, considering scalar perturbations and a flat universe, this spacetime is described by the metric where Ψ and Φ are the Newtonian gauge potentials. The lens equation is the solution of the equation of movement for the comoving transverse separation between the perturbed null geodesic and a reference unperturbed one. The deflection at a given point in the image may be seen as the integral of multiple local deflections over cosmological distances. The local deflection of a null geodesic at a given redshift depends on the traveling time of a light ray, and thus on the sum of the metric perturbations, (Ψ + Φ) . The resulting lens equation has a simple form in the Born approximation, i.e., when the line-of-sight integral is performed along the unperturbed geodesic, which is a good approximation for small deflections. In that case the lens equation is where χ is the comoving radial coordinate. From here we can read off the local deflections as α = ψ, i , defining the lensing potential as Note that the dependence on the sum of the two potentials, which is a relativistic feature, is responsible for the well-known factor of 2 discrepancy between General Relativity and Newtonian physics predictions for the deflection of light. The convergence and shear fields may then be computed from the second-order derivatives of the lensing potential, according to Equations (171) and (172), with the following caveat. The convergence and shear fields at an image location θ have contributions from all light rays that may reach that point from a population of galaxies and not necessarily from only a single source at radial coordinate w. The converge then writes where the integration is made over the sources, which are distributed over redshift as p(χ) dχ = p(z) dz.
Cosmic shear can only be detected statistically and the most used observables are two-point correlation functions and power spectra. Considering the angular correlation function of the convergence field in the Limber approximation (i.e., neglecting contributions from pairs at large radial separation) and assuming no source clustering (i.e., assuming two independent source distributions), the convergence power spectrum may be written as The lensing power spectrum on an angular multipole scale has contributions from the power spectra of the Laplacian of the lensing potential at various redshifts z(χ) on a redshift-dependent scale k = l/χ. The functions q include the dependences on the sources distances and distributions. Notice also that the statistical properties of the convergence and shear fields are related and in particular their power spectra are identical: P κ ( ) = P γ ( ).

Testing Deviations from General Relativity
In order to test a cosmological model with cosmic shear data, we need to be able to compute the evolution of the lensing potential. In General Relativity, Einstein equations at perturbation level provide a set of evolution equations for the metric and a set of constraint equations that relate the metric perturbations (the potentials Ψ, Φ) to the perturbations of the cosmological fluid (matter density contrast δ, pressure perturbation δp, fluid velocity v, and anisotropic pressure σ).
A combination of the 00 and 0i Einstein equations relates Φ with the comoving matter density perturbation: it is a generalized Poisson equation [341]. In Fourier space, with k 2 being the spatial Laplacian, it writes where ∆ is the comoving density perturbation ∆ = δ + 3H(1 + w)v/k 2 and w is the equation-of-state parameter w = p/ρ. The ij Einstein equation gives a second constraint, relating the difference Φ − Ψ with the stress σ: it is the anisotropy equation [341], The lensing convergence is determined by the Laplacian of the lensing potential, and thus by, Inserting this in Equation (177), the dependence of the convergence power spectrum on the matter power spectrum becomes explicit. The latter needs then to be worked out from the evolution equations, possibly with the help of Boltzmann codes and N-body simulations. Now, if we consider gravity theories that are modifications of General Relativity, the above equations are modified. Similar to the PPN formalism used in solar systems tests, it is useful to encapsulate the modifications in a set of parameters that can be computed in the various theories and can be easily constrained by experiments [342]. Various simplified versions of this so-called Parameterised Post-Friedmannian formalism have been proposed [343][344][345][346], where the modifications are described with two independent parameters. These are the gravitational screening Q = G eff /G, which parametrizes a change in the local gravitational force on cosmological scales, and the gravitational slip η = Ψ/Φ − 1, which parametrizes an anisotropy responsible for a difference between the two potentials. The 2 parameters may be functions of scale and redshift and their values in GR are Q = 1 and η = 0. This simple relations are valid for a large number of theories of gravity. Inserting these parameters in the constraint equations Equations (178) and (179), we can combine them to get Inserting this in Equation (177), the convergence power spectrum may be written as, Here it was assumed that ∆ ∼ δ and the Ω m factor comes from the matter density in Equation (180). The parameter Σ(a) = Q(a)(1 + η(a)/2) is the combination of PPF parameters cosmic shear is sensitive to [346] . The amplitude of the power spectrum depends thus on a degenerate combination of: modified gravity via constraint equations (Σ); distances to the sources (functions g, which depend both on the knowledge of the source redshifts and on the cosmological background); matter density (Ω m h 2 ); amplitude of the matter power spectrum (σ 8 and structure formation in modified gravity); and external systematic effects.
We applied cosmic shear data from the COSMOS survey to constrain the Σ parameter [347]. The COSMOS cosmic shear survey covers an area of 1.64 deg 2 where cosmic shear angular correlations were measured on a narrow scale range going up to 20 arcmin [348]. It is however deep, with a limiting magnitude i = 26.7 and a high density of 76 galaxies/arcmin 2 , allowing one to cross-correlate populations at different redshifts and perform a tomographic analysis [349]. We assume Σ is scale-independent in the narrow range of scales probed. Following [344], we also assume that the departure from GR originates an effective smooth dark energy Λ with constant w = −1, which implies the background evolution is identical to ΛCDM and the PPF parameters deviate from GR proportionally to ρ DE (a)/ρ m (a). This implies that Q(a) = 1 + Q a a 3 , η(a) = η a a 3 .
Assuming the PPF parameters remain small, the evolution of Σ also scales as Σ(a) ≈ 1 + Σ a a 3 . The last item needed for the evaluation of Equation (182) is the evolution of the matter overdensities, i.e., the growth function. We consider the system of the evolution equations and constraint equations for the case of sub-horizon scales, assuming time derivatives of the perturbations are negligible compared to the spatial derivatives, and ∆ ≈ δ. The evolution equation of the cold dark matter perturbation greatly simplifies, containing only a Hubble drag term and a term on k 2 Ψ [350]. In this regime, the growth parameter approximation [351]  = Ω m (a) γ − 1 .
The growth parameter γ has the value γ = 6/11 in GR and in PPF-parameterized modified gravity models [352].
With these assumptions, we attempt to constrain the parameters using only the cosmic shear data in a Bayesian analysis. The cosmic shear data consists of 15 two-point cross-correlation functions, between 6 galaxy redshift bins. All 15 functions have similar shapes but different amplitudes depending on their redshifts. Notice that our parametrization only affects the amplitude of the cosmic shear functions, since we work in the ansatz that Σ and γ are scale-independent. This implies that our constraint will necessarily be a degenerate combination of the 2 modified gravity parameters, the redshift and the standard cosmological parameters affecting the amplitude of the power spectrum. For this reason, we need to use either additional data to constrain the standard cosmology or to introduce strong priors. We introduce a prior in the form of importance sampling [353]. In particular, we take WMAP7 ΛCDM Monte Carlo Markov Chains [354] and evaluate a cube of models in a (Q a , η a , f z ) space for each point of the MCMC chain. In this way, the ΛCDM MCMC weights will be changed as a function of the cosmic shear likelihoods in the cubes of models. The (Q a , η a ) likelihoods are marginalized over the redshift-calibration nuisance parameter f z that accounts for the source redshift uncertainty. The shear correlation functions are computed using the transfer function of Eisenstein & Hu [355] with the modified linear growth γ and the non-linear halofit prescription [356]. Besides addressing the degeneracies with the standard cosmological parameters and the redshift uncertainties, we also minimise contamination from intrinsic alignments by not using cosmic shear auto-correlation functions in the analysis.
The result obtained in this analysis is compatible with GR. We find the following 1-σ constraint for the following combination of the 2 parameters:

Testing Gravity with Future Cosmic Shear Data
To take full advantage of future high-precision cosmic shear data, accurate modelling is needed at several levels, and the analyses become much more complex than the one described above. When studying specific modified gravity theories, structure formation needs to be worked out for those models, including non-linear scales, which calls for model-specific perturbative analytical approaches [357] or N-body simulations [50]. Numerical simulations are also required for other purposes. Lensing simulations are needed to estimate the covariance of lensing observables, including cosmic variance. Hydrodynamical simulations are increasingly needed to model various baryonic effects on lensing observables, such as supernova and AGN feedback, star formation or radiative cooling [51].
Working in a model-independent approach requires detailed modelling of the deviations in Einstein equations. Besides applying the PPF parameters, it is also interesting to parametrize the deviations at the observables level. [358] parametrizes deviations of the lensing power spectrum from the non-linear GR matter power spectrum. The latter can be computed from the linear one with a fitting formula [359] or from templates [360]. Scale and time dependencies of the parametrized deviations may also be probed in bins by performing a principal components analysis [361].
Cosmic shear probes deviations from GR through the integrated effect of Ψ and Φ on null geodesics. It has a high signal-to-noise on the scales where the deviations are maximal and is thus more promising than other probes of Φ + Ψ, such as the integrated Sachs-Wolfe effect [362]. However, to be able to break the degeneracy between the two potentials in the modified gravitational effects, cosmic shear needs to be combined with a local probe of matter overdensities, which are sensitive to Ψ alone. One possible combination of cosmic shear and redshift-space distortions is the so-called E G test [363,364], which has the additional advantage of being independent of galaxy bias. A first detection of the E G signal was made in [365]. The combination of weak lensing and galaxy clustering observables is an important asset of the forthcoming Euclid space mission to be launched in 2020, which combines a weak lensing and galaxy clustering wide survey and aims at measuring the rate of cosmic structure growth to a 2% precision [49].

Angular Distribution of Cosmological Parameters as a Measurement of Spacetime Inhomogeneities
On large scales, the matter distribution is statistically homogeneous. Locally however, matter is distributed according to a pattern of alternate overdensed regions and underdensed regions. Since averaging inhomogeneities in the matter density distribution yields a homogeneous description of the Universe, then the apparent homogeneity of the cosmological parameters could also result from the averaging of inhomogeneities in the cosmological parameters, which would reflect the inhomogeneities in the density distribution. In the context of backreaction models, angular variations in the parameters could also source a repulsive force and potentially emulate cosmic acceleration. Hence the reasoning was to look for these inhomogeneities across the sky and then to use an adequate toy model to compute the magnitude of the acceleration derived from angular variations of the parameters compared to the acceleration driven by a cosmological constant.
Recently, a method was proposed to search for inhomogeneities based on a local estimation of cosmological parameters [71]. We demonstrated the method with supernovae (SNe), since they are a sensitive probe of the late-time expansion history of the Universe by means of the dependence of the luminosity distance on both the geometry and the growth of structure. Hence the parameters to be estimated are those that affect the luminosity distance, namely, {Ω M , Ω Λ , H 0 }.
Supernovae have previously been used to constrain inhomogeneities in the cosmic expansion, namely, by measuring the hemispherical anisotropy of H 0 [367], by assessing the dependence of H 0 with the position on the cosmic web [368] and by mapping {H 0 , q 0 } in the sky [369]. Supernovae have also been used to constrain inhomogeneities in the dark energy, namely by constraining fluctuations of a dynamical dark energy [370] and by constraining radial inhomogeneity for a Lemaître-Tolman-Bondi metric with a cosmological constant [371].
We used the type Ia supernova sample compiled by the SNLS-SDSS collaborative effort Joint Light-curve Analysis (JLA) [366] , totalling N SNe = 740 SNe with redshift 0.01 ≤ z ≤ 1.30 and distributed on the sky according to Figure 13. The sample was obtained from http://supernovae.in2p3.fr/sdss_snls_jla/ReadMe.html.

Method of Local Parameter Estimation
Our method consists in dividing the SNe over a pixelated map of the sky, with pixels of equal surface area according to the HEALPix pixelation scheme [372]. The SNe subsample that falls in each pixel k is used to estimate the cosmological parameters in that pixel, thus the designation of local estimation. For each cosmological parameter being estimated, there results a map with the same pixelation as that of the SNe subsamples. Each pixel is assumed to be described by a Friedmann-Lemaître-Roberston-Walker metric so that the full sky is an inhomogeneous ensemble of disjoint, locally homogeneous regions. In each pixel, we estimated the cosmological parameters that minimise the chi-square of the fit of the theoretical distance modulus i . Since the SN sample used is a collection from different SN surveys, the SN subsampling in pixels is highly inhomogeneous. As a result, well sampled pixels return robust estimations, whereas poorly sampled pixels return weak constraints. This generates a noise bias. We hypothesised that the noise bias due to the inhomogeneity in the SN subsampling is modelled by the local estimation obtained from randomizing the dependence of redshift with position and then averaging over the different randomizations. For each randomization in z, the SN positions in the sky remain the same as in the original sample. Subtracting the noise bias off the maps previously estimated we obtained unbiased maps x unbias ik ≡ x ik −x bias ik + x fid i , wherex bias ik ≡ x bias ij k j is the average over j randomizations in z. For better visualization, we define the unbiased difference maps as ∆x unbias ik ≡ x unbias ik − x fid i = x ik −x bias ik . In Figure 14 we show the values of the difference maps at each pixel, before (left panel) and after (right panel) the noise bias subtraction. We also show the corresponding errors. We observe a negative correlation of the number of SNe in each pixel both with the fluctuations about the fiducial values, and with the corresponding errors. These results reflect the poor sampling derived from an inhomogeneous coverage of the sky, both in angular and redshift space, by the SN surveys. A poor sampling into pixels implies a poor sampling in redshift, which limits the constraint of the degeneracies among the parameters and hence compromises the parameters estimation.
We also observe that, after the bias subtraction, the fluctuations about the fiducial values are reduced, albeit the errors increase. Consequently, for the correlation of the number of SNe per pixel, the unbiased maps in comparison with the original maps show a smaller correlation with the difference values and a larger correlation with the errors, as expected. These results indicate the validity ofx bias ik as a measure of the noise bias due to the inhomogeneity of the SN sampling. Comparing the unbiased difference maps with the fiducial values, we measure fluctuations of order 5%-95% for Ω M , 1%-25% for Ω Λ and up to 5% for H 0 .

Average over the Local Parameter Estimation
When averaging over the homogeneous regions, angular fluctuations in the expansion factor (and consequently in H 0 ) induce a backreaction term in the form of an extra positive acceleration [373]. We derived the analytical extra positive acceleration for a toy model of an arbitrary number of disjoint, homogeneous regions and computed the overall deceleration parameter assuming a) no backreaction and b) backreaction for the measured angular distribution of H 0 .
In the absence of backreaction, the averaging consists in taking the mean weighted by the variance's inverse w k = 1/Var[x ik ] of parameter x i in each pixel k, For the derived parameters, we find Ω κ = 0.078 ± 0.131 andq 0 = −0.451 ± 0.159, and after the noise bias removal we find Ω κ,unbias = 0.002 ± 0.135 andq 0,unbias = −0.526 ± 0.172. After the noise bias removal, the pixel average values of both Ω κ and q 0 become closer to the corresponding values estimated using the complete sample.
In the presence of backreaction, the averaging consists in taking the mean weighted by the 3-volume V k of each pixel k,x Since all pixels have the same surface area and hence the angular directions expand the same way today, we can assume that the radial direction also expands the same way today, which implies that the 3-volume is the same today for all pixels. This amounts to identifying a volume as a pixel.
Defining v k = V k / ∑ k V k , then for N pixel disjoint regions, the average of the Hubble parameter is given by Taking the derivative, it follows that which decomposes into a linear term in the pixel average of (ä/a) and a quadratic term in differences of H 0 between pairs of pixels. The quadratic (backreaction) term generates an acceleration that is due not to regions speeding up locally, but instead to the slower regions becoming less represented in the average. Note that, in the absence of the quadratic term, the volume average reduces to the pixel average above. Then the volume average of q 0 becomes Using the fluctuations in H 0 measured in the local parameter estimation, we findq 0 = −0.452 ± 0.159 and after the noise bias removal we findq 0,unbias = −0.527 ± 0.173. The volume average of the other parameters is the same as the pixel average. After the noise bias removal, the volume average values become closer to the value for the deceleration estimated from the complete sample. The results from both averaging methods strengthen the validity ofx bias ik as a measure of the noise bias due to the inhomogeneity of the SN sampling.
The quadratic term is of order 10 −3 (or lower) times the linear term, which means that the error of the difference is below the standard deviation. Hence for the angular fluctuations in H 0 measured with this SN sample, the contribution of the quadratic term in Equation (192) is insignificant, which renders the volume averaging equivalent to the pixel averaging. Hence in the context of this toy model of an inhomogeneous spacetime, backreaction is not a viable dynamical mechanism to emulate cosmic acceleration. A subsequent study using a kinematic parametrisation of the luminosity distance expressed in terms of time derivatives of the scale factor yielded concordant results [374].

Summary and Conclusion
In this work, we have explored the dynamics and evolution of the Universe at early and late times, focusing on both dark energy and extended gravity models and their astrophysical and cosmological consequences. More specifically, we presented scalar-tensor theories, focussing on a brief review of the general formalism, the conformal picture, discussed the experimental and observational status of the theory, presented a detailed dynamical system analysis of the cosmological dynamics, and concluded with a quintessence scenario where a multi-scalar-tensor theory is responsible for the present accelerated expansion of the Universe. Furthermore, we considered Horndeski theories and the self-tuning properties, and explored realistic cosmologies that can be constructed, by studying the properties of this intriguing theory. We also reviewed f (R) modified theories of gravity and extensions, where the foundational questions and astrophysical/cosmological issues of the Palatini approach were extensively explored; we also focussed on recently proposed theories such as on curvature-matter couplings and the hybrid metric-Palatini theory. The inflationary epoch was also presented in f (R) gravity, where the so-called Starobinsky inflation is known as one of the more reliable candidates for explaining the inflationary paradigm as provided by the recent Planck/Bicep2 constraints. Moreover, we have reviewed the evolution and cosmological consequences of topological defect networks that may be formed in the early universe. In particular, we focused on two of the most significant signatures of these networks: the cosmic microwave background anisotropies and the stochastic gravitational wave background.
As mentioned in the Introduction, due to the extremely important aspect of the synergy between theory and observations, the second part of this review was dedicated to observational cosmology. In particular, cosmological tests with galaxy clusters at CMB frequencies were presented. More specifically, galaxy clusters and the thermal and kinetic Sunyaev-Zel'dovich (SZ) effects were reviewed, and the cosmological tests using SZ cluster surveys were presented, focussing on SZ cluster counts, SZ power spectra and the possibility of probing new physics with the SZ clusters; a review on probing primordial non-gausianities with galaxy clusters was also presented, with an emphasis on the parametrization of primordial non-Gaussianities, the non-Gaussian mass function, biased cosmological parameter estimation with cluster counts, the impact of primordial non-Gaussianities on galaxy clusters scaling relations. Furthermore, gravitational lensing was also explored, with a focus on cosmological scales, with the possibility of testing deviations from General Relativity, in particular, with future cosmic shear data. In conclusion, the angular distribution of cosmological parameters as a measurement of spacetime inhomogeneities was presented, using a method of averaging over the local parameter estimation. Thus, cosmological tests were discussed, and their relevance in constraining our cosmological description of the Universe.
(EC) by FEDER funding through the program Programa Operacional de Factores de Competitividade -COMPETE. The authors also acknowledge the COST Action CA15117, supported by COST (European Cooperation in Science and Technology).