Cosmology for quadratic gravity in generalized Weyl geometry

A class of vector-tensor theories arises naturally in the framework of quadratic gravity in spacetimes with linear vector distortion. Requiring the absence of ghosts for the vector field imposes an interesting condition on the allowed connections with vector distortion: the resulting one-parameter family of connections generalises the usual Weyl geometry with polar torsion. The cosmology of this class of theories is studied, focusing on isotropic solutions wherein the vector field is dominated by the temporal component. De Sitter attractors are found and inhomogeneous perturbations around such backgrounds are analysed. In particular, further constraints on the models are imposed by excluding pathologies in the scalar, vector and tensor fluctuations. Various exact background solutions are presented, describing a constant and an evolving dark energy, a bounce and a self-tuning de Sitter phase. However, the latter two scenarios are not viable under a closer scrutiny.


I. INTRODUCTION
General Relativity (GR) is the most widely accepted theory of gravity. It has been subjected to intense experimental confrontation, but yet there is no direct evidence that would signal a need for modifications of the theory [1]. Nevertheless, research of different gravity theories is flourishing. This is partially because GR is not found to be theoretically satisfactory: it is not a quantum theory, and it predicts singularities, for example. Also, coupled with the Standard Model of elementary particles, it predicts a catastrophically large cosmological constant, if one does not resort to technically highly unnatural fine-tuning that is unstable against quantum corrections [2,3]. Furthermore, if we accept that GR is valid at all scales, we then need to invoke new ingredients in the cosmological matter sector: the notorious dark energy and the dark matter that has escaped all our attempts at direct or indirect detection [4][5][6][7].
A defining property of GR is the unity of the metric and the affine structures of spacetime. In GR, the latter is tied to be fully determined as the Christoffel symbols of the metric tensor which is thus the only independent field. Most modifications of GR have been introduced in this context, by considering then more general actions for the metric and by introducing additional fields that couple non-minimally to the metric. The cosmological constant problem, for example, might be addressed by devising an action for the metric that would degravitate [8][9][10][11], adjust suitably the underlying symmetries such as by introducing unimodularity [12,13] or feature self-tuning with or without extra fields [14][15][16]. On the other hand, the ultraviolet singularities and the obstacles to renormalisation could be smoothened by quantum effects or by higher-derivative terms for the metric in the gravitational action [17][18][19][20].
In the more general context of non-Riemannian geometries, richer possilibilities emerge as the connection may carry also torsion [21] and non-metricity [22][23][24]. As also Einstein himself argued, the connection appears to have a more fundamental ontological status than the metric field [25], which has become apparent especially in the developments of gauge theories of gravity [26,27]. In the Palatini approach, wherein one regards the spacetime connection as an independent field besides the metric, it is well known that GR is dynamically recovered for the Einstein-Hilbert action (and all the subsequent Lovelock invariants [28,29]), but more general actions become different theories when subjected either to the metric or to the Palatini variation [30][31][32][33]. A perhaps more conservative than the Palatini approach is however to consider that only some of the degrees of freedom residing in the non-metricity or the torsion sector are physically relevant. A prototype example is the Weyl geometry, wherein the only additional field is the trace of the non-metricity, the so called Weyl vector.
Since the first conception of the idea of gauge symmetry, Weyl geometry has provided an important framework for the development of fundamental theories [34]. In Weyl's unified theory of electromagnetism and gravity, the principle of relativity applied not only to the choice of reference frames, but also to the choice of local units of length. The corresponding spacetime geometry can accommodate scale-invariance, meaning that physics becomes insensitive to for example the difference of masses of various particles. Though the original formulations of Weyl's unified theory can be abandoned as not viable, the idea that physics at some fundamental level is scale invariant remains alive and appealing. Obviously, such a symmetry should be somehow broken at our low energy world, but it still could provide a key to a solution of the most fundamental problems of gravitational physics, if it is accepted that the elimination of the physical propagation of the conformal mode of GR could redeem gravity from both the ultraviolet singularities and the infrared vacuum's weight [35,36]. It was recently argued that local conformal invariance has to be an exact symmetry and further, broken in a spontaneous manner [37].
In this paper we study cosmologies in an extended Weyl geometry. It is not our aim to formulate a scale-invariant theory (for such, see e.g. [38,39]) but simply explore the cosmological implications of a fundamental gravitational vector mode. It is perhaps surprising that the leading order, quadratic curvature corrections to the gravitational action were elucidated in the context of Weyl geometry only very recently [40,41]. The result is interesting as it introduces a novel four-parameter class of potentially viable (at least, ghost-free) vector-tensor theories. We considered an extension of the Weyl geometry wherein the distortion (i.e. the non-Levi-Civita part of the connection) is linear in a vector field in the most general non-derivative way, featuring then two extra terms besides the pure Weyl-type non-metricity [42]. As one of the special cases, this geometry includes the one-parameter non-metric spaces of Ref. [43]. Here we will systematically derive the quadratic curvature theory in the geometry with linear vector distortion (for studies of other quadratic non-metric theories, see e.g. [44,45]), and apply the results to study various new cosmological scenarios. Vector field cosmologies have been already studied extensively in the literature, with interest in e.g. dark energy [46][47][48], inflation [49][50][51], anisotropies [52][53][54][55], dual fields [56,57], Gauss-Bonnet couplings [58,59], screening [60,61], a cosmological constant [62][63][64] and stability [65,66]. Our framework however suggests a fundamental geometric origin for the possible existence of a cosmic vector, furthermore predicting a rather well-specified vector-tensor action.
We begin in the next Section II by deriving this action from the most general quadratic-in-curvature theory, on which we then impose restrictions by various consistency requirements (we also take a brief look at possible higher order curvature invariants). Remarkably, it turns out that we are then restricted to the class of generalised Weyl geometry, whose special status was already recognised from another aspect [42]. In Section III we then study the existence of de Sitter solutions, and find out that there are 0-4 such solutions depending on the theory parameters. We go further and analyse the propagation of scalar, vector and tensor fluctuations for generic models in de Sitter background. In Section IV we then study in more detail some specific analytical solutions that are found for some given parameter combinations. We present both interesting late time (dark energy, cosmological constant) and early time (self-tuned de Sitter with first order phase transition, bouncing solution) scenarios. Finally, we conclude in Section VI and complete some derivations with details in the Appendix.

II. GENERALIZING WEYL GEOMETRY: SPACETIMES WITH LINEAR VECTOR DISTORTION
In this section we will start by briefly reviewing the basic properties of Weyl geometry and how it can be generalized to include the most general connection linearly determined by a vector field as introduced in [42]. Then, we will proceed to the construction of gravitational actions based on these geometries.

A. Geometrical framework
The defining property of Weyl geometry is the breaking of the metricity condition by introducing a vector field as follows∇ This relation is preserved by the transformation g µν → e 2Λ(x) g µν when simultaneously we transform A µ → A µ −∂ µ Λ(x) and, thus, Weyl geometry is a natural arena to formulate conformally invariant theories, although we will not pursue this here. The above expression can be easily solved for the connection (assuming vanishing torsion) to obtain We see that in Weyl geometry, the connection acquires a distortion linearly depending on A µ .
The class of geometries introduced in [42] extends the Weyl geometry by allowing for the most general connection with a distortion tensor linearly determined by a vector field. The form of the desired connection is thuŝ where Γ α βγ are the usual Christoffel symbols of the metric and b i are arbitrary coefficients. This connection gives rise to non-metricity, but also contains the trace vector part of the torsion. Interestingly, although the conformal invariance of the metric compatibility condition is lost by the general connection (3), the extra torsion component allows to recuperate it for more general connections. To see it explicitly, we can compute the covariant derivative of the metric, which is given by∇ Thus, we recover the Weyl condition (1) provided we impose 2b 1 − b 2 − b 3 = 0. In the absence of torsion b 3 = 0, we exactly recover the Weyl condition, but the torsion b 3 term allows to maintain the gauge invariance of (4) for more general connections. It is not difficult to see that We can now introduce the Riemann tensor of our connection as usual It is important to keep in mind that, since we have torsion and non-metricity, this Riemann tensor does not have the usual symmetries of Levi-Civita connections, but only the antisymmetry in the first two indices inherited from its definition as a commutator: R µνρ α = −R νµρ α . In particular, this means that we can construct 3 independent traces, namely: the usual Ricci tensor R µν ≡ R µαν α , the co-Ricci tensor P µ α ≡ g νρ R µνρ α and the homothetic tensor 1 Q µν ≡ R µνα α . For the connection given in (3), these tensors can be expressed as where R µν and ∇ are the Ricci tensor and covariant derivative of the Levi-Civita connection of the spacetime metric, F µν = ∂ µ A ν − ∂ ν A µ is the strength tensor of the vector field and D is the spacetime dimension. Notice that the Ricci tensor R µν is not symmetric, not even in the torsion free case with b 3 = 0, since non-metricity can also induce an antisymmetric part for the Ricci tensor. It is also convenient to keep in mind that the homothetic tensor is always antisymmetric and for symmetric connections it is proportional to the antisymmetric part of the Ricci tensor. Finally, the Ricci scalar is unambiguously defined as R = g µν R µν = −P µ µ and, for our connection (4), it is given by where R is the Ricci scalar of the Levi-Civita connection.
To end this section we will give some important geometrical objects. A defining property of geometries with nonmetricity is that the length of vectors is not preserved under parallel transport. If we have a vector v α , its length v 2 = g αβ v α v β under a parallel displacement dx µ is Obviously, for a pure metric geometry with A µ = 0, the length is conserved. However, it is remarkable that the presence of torsion (b 3 ) also allows to have a wider class of geometries that preserve the length of vectors given by This actually allows to avoid one of the main problems of the pure Weyl geometries, where the length of a vector changes as it is parallel transported and, consequently, the properties of a physical object may depend on its history 2 . In our more general geometrical set-up, this can be avoided even for non-trivial connections with non-metricity. One could relax the condition of invariance of the length of vectors only when they are transported along closed loops. In that case, the variation is determined by the homothetic or segmental curvature tensor Q µν introduced above. Again, we have a special family of geometries with 2b 1 − (D + 1)b 2 + (D − 1)b 3 = 0 in which the length of a vector remains invariant when transported around a closed path. Of course, these geometries contain the aforementioned case with b 3 = b 2 = b 1 , but more general cases are possible in which the length of a vector may vary under a parallel transport while remaining the same when the trajectory closes.
Besides the variation of the length of vectors, the presence of non-metricity will also affect the properties of the geodesics. Since only the symmetric part of the connection enters the geodesic equation, the torsion term determined by b 3 will not enter here. A crucial feature of the geodesics is their projective similarity, i.e., two families of geodesics related by a change of affine parameter will describe the same class of paths. More precisely, it can be shown that two connections differing by δ α (µ ξ ν) give rise to the same geodesics up to a redefinition of the affine parameter [67]. The projective invariant object determining the class of paths is the Thomas projective parameter, which for our connection is given by: where Π α µν is the piece corresponding to the Levi-Civita part of the connection. In the above expression, only the symmetric (torsion-free) part of the connection should be considered. As expected, this expression does not depend on b 2 , since this term precisely corresponds to a projective transformation and, consequently, the geodesic trajectories will not be affected by it.

B. Gravitational actions
Theories based on the Ricci scalar given in (9) lead to interesting phenomenologies, including the Starobinsky inflationary model and its so-called α-attractor generalisation [42,68,69]. This is possible because in f (R) type of theories with (9) the vector field is dynamically constrained to be the gradient of a scalar and, thus, similarly to usual f (R) theories, they are equivalent to a scalar-tensor theory.
In order to have a fully dynamical vector field (and not imposed to be a pure gradient) we need to consider more general actions. The natural step is then to include quadratic curvature terms. For Levi-Civita connections, the requirement of having second order equations of motions (so that Ostrogradski instabilities are avoided) leads to the well-known Lovelock invariants. The quadratic term of such invariants is the Gauss-Bonnet term, which in 4 dimensions is a total derivative and, thus, it does not contribute to the gravitational field equations. Considering the corresponding Gauss-Bonnet term with our vector connection will also give a total derivative in 4 dimensions. However, as found in [40], more general quadratic curvature terms can give rise to a new interesting class of vectortensor theories in the context of Weyl geometry, see also [41]. Following the same procedure, we will write down the most general action quadratic in curvature invariants in the more general linear vector geometry defined by the connection (4). Such terms are given by where c i and b i are dimensionless constants and µ has dimension [mass] 4−D . In order not to have Ostrogradski instabilities associated to higher order equations of motion for the metric we will restrict the parameters in order to recover the Gauss-Bonnet term for the Levi-Civita part of the connection (i.e, when A µ = 0) so that Since Q µν identically vanishes for A µ = 0, the coefficients c 7 , c 8 , c 9 remain fully free. Now we can rewrite (12) as a vector-tensor theory in a Riemannian geometry. To that end, we will use the decomposition of the connection (4) in (12) and express everything in terms of A µ and the Levi-Civita connection Γ. After some straightforward algebra and a few integrations by parts as done in [40], the action can finally be expressed as where α, β, γ i , λ and ξ are dimensionless parameters that depend on b i , c i , d i and the spacetime dimension D. The combination in the brackets in the first line is nothing but the Gauss-Bonnet term for the Levi-Civita connection, as a consequence of having imposed (13). The remaining terms in the first line are the same as were obtained in [40] for the case of pure Weyl geometry. Despite the derivative interaction ξ and the non-minimal coupling β, those terms only propagate 3 degrees of freedom, very much like the simpler case of a Proca field. However, the terms in the second line were not present in the pure Weyl case and only arise for our general connection (4). These are however undesirable because they will propagate one additional degree of freedom besides the 3 polarizations corresponding to a massive vector field and this extra mode will generically suffer from the Ostrogradski instability, i.e., it will be a ghost. Thus, in order to have a stable theory we need to impose the conditions γ i = 0. The details are given in the Appendix A. Remarkably, the only solution for D ≥ 4 is b 3 = 2b 1 − b 2 , which exactly coincides with the generalised Weyl geometry discussed above that preserves the local conformal invariance of the metric compatibility condition.
In that case, we can canonically normalize the field by means of the rescaling Pl R for completeness, which will simply give a mass term for the vector field, so the final action for the vector reads with As we see, in 4 spacetime dimensions only the mass term remains while all the interactions vanish. To obtain a non-trivial theory in 4 dimensions we could consider the limit D → 4 simultaneously with α → 0 while keeping (D − 4)/α and αµ fixed. Then there remains only two free parameters, in particular the action could be reduced to (when choose the free parameters as the constants ξ and β) Throughout this work, we will focus on D = 4 and for generality consider all 4 parameters as independent. It can be useful to note that our general vector-tensor action can be rewritten in an alternative way by using that where we can recognize the typical Horndeski form for the vector field derivative self-interactions.
So far we have only considered actions up to quadratic order in curvature invariants. The same program can be straightforwardly applied to higher order curvature invariants by imposing the terms at a given order to reduce to the corresponding Lovelock invariant. Instead of discussing the general framework for this construction, we will simply comment on a class of terms that give rise to additional non-trivial derivative interactions for the vector field. For simplicity, we will consider the case D = 4. At cubic order we know that a healthy coupling is given by the Horndeski vector-tensor interaction [70] L µναβ F µν F αβ , with L µναβ = − 1 2 ǫ µνρσ ǫ αβγδ R ρσγδ the double dual Riemann tensor, so we can use it to construct new terms in our class of geometries. For that, we notice that the antisymmetric part of the 3 independent traces of the Riemann tensor given in (8) are all proportional to F µν . Thus, we can consider an interaction of the form that will give the Horndeski interaction for the non-minimal derivative coupling of the vector. Although we have written down explicitly all the possible terms, they all will contribute exactly the same interactions for the vector field. Explicitly, the above cubic interaction leads to whereF µν = 1 2 ǫ µναβ F αβ is the dual of the vector field strength tensor and µ 3 some energy scale. The first term is the advertised Horndeski vector-tensor interaction (which has been studied in, e.g., [71,72]), while the second term in the first line is a derivative self-interaction of the vector field. Although this term does not respect gauge invariance and contains derivatives of the vector field, its structure does not spoil the constraint making A 0 non-dynamical. This can be seen by noticing that the time derivative of A 0 couples toF 0iF 0 i , which is proportional to the magnetic part of F µν . Since this magnetic component only contains spatial gradients, A 0 will not acquire second derivatives in the field equations and, consequently, it will not propagate. These interactions are in fact within the class of derivative self-interactions for a massive vector field discussed in [73][74][75][76]. Moreover, it is expected that by considering higher order curvature terms within our framework, the higher order terms interactions introduced in these references will be generated.

III. DE SITTER SOLUTIONS
In this section we will show the existence of isotropic de Sitter solutions for the vector-tensor theory given by (15). In order to comply with the given symmetries, we will consider a purely temporal and homogeneous vector field A µ = A 0 (t)δ µ 0 configuration as well as a homogeneous and isotropic metric described by the Friedmann-Lemaître-Robertson-Walker (FLRW) line element where a(t) is the scale factor. The only non-trivial vector field equation for this configuration is given by with H =ȧ/a the Hubble expansion rate. As expected, this is an algebraic equation showing that A 0 is not dynamical and it is fully determined by H. This equation leads to 3 branches, namely: the trivial one with A 0 = 0 and 2 nontrivial ones with From this expression we can see that de Sitter solutions are not guaranteed for any values of the parameters. For the singular case λ = 0, the degree of the equation is reduced and only one non-trivial solution remains. This particular case will be studied separately below. In order to obtain the value of the Hubble parameter for the de Sitter solution we need to look at the corresponding Friedmann equation For the trivial branch with A 0 = 0 we obtain that H = 0, i.e., we recover the Minkowski solution as expected. On the other hand, for the non-trivial branch given by (26), the Friedmann equation gives an algebraic equation for H. Such an equation forces H to be constant and, therefore, these branches actually correspond to a de Sitter universe. However, the existence of such solutions will be subject to the existence of real solutions for the corresponding system of algebraic equations, which imposes restrictions on the parameters as already mentioned above. In fact, we can see that the equations reduce to a system of two polynomial (one quadratic and one quartic) equations so that we will have in general up to 8 different branches. They have some properties that can simplify the analysis of the solutions. First, one can easily see that there is a symmetry in the equations H is a solution for ξ, then −H will be a solution for −ξ. Thus, without loss of generality we can absorb ξ into a rescaling of H → H/ξ and focus on solutions with H > 0 (i.e., expanding solutions), keeping in mind that a corresponding contracting solution will be guaranteed to exist as well. Further, one can get rid of another parameter (up to its sign) by absorbing it into the normalization of A 0 . The most convenient one is to absorb β so that A 0 → A 0 / |β|. In order to characterize the models for which de Sitter solutions exist, we will introduce the following rescalings and a dimensionless variable x Obviously, this can only be done if both β and ξ are non-vanishing. From now on, we will assume this and treat separately below the case of vanishing β or ξ. Then, from the vector field equation we can solve for H in terms of x as This already allows us to put some conditions on the parameters to have de Sitter solutions, since we need to have H 2 > 0. Now we can plug this solution into the Friedman equation to obtain the following quartic equation for x: This is the crucial equation that will determine the existence of de Sitter solutions. From here we see that if λM 2 is negative, there will be at least a couple of real solutions. It is important to notice that, since this equation does not depend on β, we can take its real solutions, plug them into (29) and choose the sign of β in order to guarantee that H 2 is positive. In Fig. 1 we show the number of solutions in the parameter space spanned by (λ, M 2 ). For the singular case λ = 0 it is not difficult to see that there is 1 solution for M 2 > 0 and 3 for M 2 < 0.
We can see that, quite generally, the vector-tensor theories considered here give rise to de Sitter solutions. The next natural question is then if such solutions are stable. This is the subject of the next subsection.

Stability
After showing the existence of de Sitter solutions, we will now proceed to check the stability of such solutions. For that, we will study the behaviour of the inhomogeneous perturbations around the de Sitter solutions found above. Thus, the background will be given by a constant temporal component for the vector field A 0 and a constant Hubble expansion rate H. The perturbations for the metric g µν will be decomposed into irreducible representations of the SO(3) symmetry of the background in the usual manner [77] where it is understood that all the metric perturbations depend on time and space and we have the constraints On the other hand, the vector field will be analogously perturbed as with δ ij ∂ i δA j = 0. Naïvely counted, we encounter 14 dof's in this decomposition, namely: 2 × 1 in the traceless symmetric tensor (h ij ), 2 × 3 in the divergence-free vectors (B i , E i , δA i ) and 1 × 6 in the scalars (Φ, B, ψ, E, δA 0 , A s ). However, we have 4 diffeomorphism gauge symmetries that remove 2 dof's each, so we have 14 − 2 × 4 = 6 dof's. In addition, as discussed above, the temporal component of the vector is not dynamical, but an auxiliary field, so we should substract yet another dof and, thus, the number of physical propagating modes will be 6 − 1 = 5, i.e., the 2 polarizations corresponding to the massless graviton plus the 3 polarizations of the massive vector. In the following we will explicitly see this by studying the tensor, vector and scalar perturbations individually, and further establish the stability conditions for the 2 tensor, 2 vector and 1 scalar dof's in the theory. Due to the homogeneity of the background, it will be convenient to decompose the perturbations in Fourier modes with respect to the spatial coordinates. Hence, all perturbations will be expanded as where Θ represents a given perturbation. As a final remark, let us stress that we will perform the analysis without a priori fixing the gauge at the level of the action. We have also carried out the analysis in the Newtonian gauge where B = E = B i = 0 as a consistency check and found the same results.

Tensor perturbations
Let us start by studying the tensor perturbations. This is the simplest case because there are no dynamical dof's to be integrated out. After inserting the metric perturbations (31) into our vector-tensor action with the fields decomposed in Fourier modes and using the dynamical background equations from the previous section, the action quadratic in the tensor perturbations reads We see that the quadratic action for tensor perturbations is modified in the presence of a background A 0 and, as one would expect, only the non-minimal coupling to the Einstein tensor (β-term) contributes. However, there is a dependence on the remaining parameters through the background value A 2 0 , which is determined by all the parameters. From the above action, we easily conclude that tensor perturbations around the de Sitter solutions will avoid ghost-like instabilities, i.e., they have the right sign for the kinetic term, if we impose which is trivially satisfied if β ≥ 0. On the other hand, the propagation speed of the perturbations is that also must be positive to avoid gradient instabilities. These results are in agreement with those found in [78]. We can avoid both ghosts and gradient instabilities for the tensor perturbations if In any case, this non-trivial effect on the propagation speed of gravitational waves will be tightly constrained at present time since binary pulsar observations put stringent limits, at the level of 10 −2 deviations from the speed of light [79].

Vector perturbations
Let us now turn to the slightly more involved case of vector modes. As we mentioned previously, two of the vector perturbations can be integrated out and only one vector will propagate (two dof's, as a transverse 3-vector). We expand our Lagrangian to second order in the vector perturbations and immediately observe that the vector field B i does not have any kinetic terms. Therefore we can simply compute the equation of motion with respect to B i and integrate them out. This yields After plugging this expression back into the action and adding total derivatives, the dependence on E i drops. Thus, as advertised, we end up with the quadratic action for only one vector, The propagation speed is given by again in agreement with the findings in [78]. As in the case of tensor perturbations, we see that only the coupling to the Einstein tensor (β-term) modifies the quadratic action of vector perturbations, and the remaining parameters only enter through the background value A 0 . We see that the vector perturbations are never ghostly and the only instability that can appear is a gradient one. This was expected because the kinetic term for our vector field is nothing but the usual Maxwell term. In order to avoid the Laplacian instability we have to require c 2 v > 0. Interestingly, it is trivially satisfied for models with β ≥ 0.

Scalar perturbations
As usual, the scalar sector is the most involved one. We have six scalars (ψ, δA 0 , A s , E, B, Φ ), but only one propagates, corresponding to the longitudinal mode of A µ . Again, we expand the action to quadratic order in the scalar perturbations. The first thing to be noticed is the fact that the corresponding kinetic matrix (or alternatively the Hessian matrix) contains already at this stage three vanishing eigenvalues, imposing three constraint equations that makes three out of the six scalar fields not propagating. The kinetic matrix is As one can see, the quadratic action does not evolve any kinetic term for the scalar fields Φ and B and δA 0 . Thus we can simply replace them by using their equations of motion. For instance, the equation of motion for the scalar field Φ gives Similarly, the expressions for δA 0 and B are obtained using their equations of motion, which we omit here. After inserting the solutions for B, Φ and δA 0 back into the quadratic action, the resulting expression depends only on the remaining three scalar fields (ψ, A s , E). On a closer inspection, one realizes that the kinetic matrix of the three scalar fields still has a vanishing determinant, pointing to the presence of more constraints that can be used to integrate out some of the scalar fields. To be precise, the kinetic matrix contains two zero eigenvalues and one non-vanishing eigenvalue. We can diagonalize the kinetic matrix by performing the following field redefinitions, which will make the only propagating scalar field manifest: After adding total derivatives, the resulting action depends only on F 3 and reads where K s and V s are some functions of the theory parameters and the background solution. Their exact expressions are very cumbersome in the general case. However, their UV limits with k → ∞ can be expressed as For the stability of the perturbations we have to impose that both K s and V s be positive. Unlike in the previous cases, we see that for the scalar sector all the terms contribute to the perturbations and not only through the background evolution.

IV. COSMOLOGICAL SOLUTIONS
After showing the existence of de Sitter solutions and studying the corresponding perturbations around them, we will now consider a more general case in the presence of a matter component. Again, we will consider homogeneous and isotropic universes with the FLRW line element. In such a background metric and for a homogeneous vector field configuration with A µ = (A 0 (t), a −1 A(t)), the vector field equations read From these equations we see that it is a consistent Ansatz to consider purely isotropic solutions with A = 0 so that our isotropic solutions will be supported by the temporal component of the vector field. Having isotropic solutions based on the spatial part of the vector usually requires to have a set of vector fields with a global SO(3) symmetry [46,80], sometimes called triad, or rapidly oscillating fields [81,82].
Let us now consider a universe filled with the vector field plus a matter component with energy density ρ. We will still restrict to purely isotropic solutions. Then, the vector field and Friedmann equations are respectively. For A 0 = 0 we have the trivial branch that recovers standard gravity so we will assume that A 0 = 0. In such a case, we can simply integrate out A 0 by (algebraically) solving its own equation of motion, whose solution is still given by (26) and, in general, it yields A 0 = A 0 (H). Then, we can use this solution again in Friedmann equation to obtain a generalised version of it as Then, we can invert this equation to obtain a Cardassian-like model where the Hubble expansion rate is given by a non-linear function of the energy density. This is nothing but a particular example of the general result that gravity with auxiliary fields leads to a modified matter coupling. The explicit expression is where ǫ = ±1 parameterizes the two non-trivial branches. If we expand this expression for small H we obtain recovering the usual GR Friedmann equation with an effective cosmological constant and a rescaled Planck mass. Notice that the case with λ = 0 is singular, as a consequence of reducing the degree of the vector field equation (which becomes linear in A 0 in that case). This singular case will be studied in more detail separately below. Also notice that the small-H regime does not need to correspond to a low density regime. A very interesting property of the cosmological evolution in the presence of the vector field is the possibility of having a maximum and/or minimum value for the Hubble expansion rate. This is determined by the discriminant of the vector field equation, i.e., the behavior of which must be positive in order to have real solutions. If both combinations λM 2 and 9ξ 2 − 24λβ are positive, we always have real solutions, while if they are negative, real solutions do not exist. However, if they have different signs, we encounter two possibilities: • If λM 2 > 0 and 9ξ 2 − 24λβ < 0, there is an upper limit for H given by • If λM 2 < 0 and 9ξ 2 − 24λβ > 0, then we have a lower bound for the Hubble expansion rate so that H ≥ H ⋆ .
The above conditions guarantee that the solution for the vector field is real. However, we need additional conditions to guarantee the existence of physical solutions because, once the vector field has been solved for, we will obtain an equation for H(t) which must also have real solutions. The overall effect will be the presence of bounds for the energy density of the matter component. The corresponding analysis for the general theory is very cumbersome, so we will focus on some particular cases instead below.
The models with an upper bound for H could be useful to resolve the Big Bang singularities. This is expected to be a more general feature (not only for cosmology) that could help regularising other types of singularities, e.g. black hole singularities. Notice that this stems from the quadratic equation for the vector field. Pl ) we can see 6 points (marked in blue) where we encounter a divergence inḢ and, thus, the evolution has a singularity. We also see a turnover matching an expanding universe with a contracting phase at low densities. Notice that at the turnover the different branches join. Also as discussed in the main text, we see that both H and ρ are bounded to a compact region. In the right panel (β = 0, λ = 1, ξ = 6 and M 2 = 2.5M 2 Pl ) we can see the type of bounce that we can have in these theories. We see that from the bounce we can evolve either towards a low density de Sitter phase or end in a singularity. This illustrates the fact that a bounce cannot be connected with a GR regime at low densities and for small H. The GR regime in this case is beyond the singularities, in a region that is not continuously connected with the bounce.

A. Bouncing solutions
An interesting question that arises is whether it is possible to obtain viable bouncing solutions. As we have shown, the presence of the auxiliary field A 0 gives rise to a modified Friedmann equation and, therefore, it would be plausible to encounter some bouncing solutions. Such solutions are characterized by the existence of a finite (non-vanishing) energy density ρ b for which the Hubble expansion rate vanishes, i.e., H(ρ b ) = 0 for ρ b = 0. If we look at Eq. (49), we see that for a potential bouncing solution we must have so that it is only possible in models with M 2 and λ different from zero and such that they have the same sign. This relation further implies that, at the bounce, the two branches of solutions for A 0 with ǫ = ±1 coincide. Thus, a generic behaviour will be the merging of branches at the bounce (this is explicitly illustrated in Fig 2). On the other hand, from the Friedmann equation we have that the potential bounce happens for an energy density so that at the bounce we have the relation ρ b = 4λA 4 0 . Moreover, we see that λ (and therefore M 2 ) must be positive in order to have a bounce with a positive density energy. However, for this bounce to actually describe a transition from contraction to expansion, we also needḢ to be positive near the bounce. This time derivative can be expressed asḢ where we have transformed the time derivative into derivative with respect to ρ and used the continuity equatioṅ ρ = −3H(ρ + p). This expression shows the well-known fact that in GR we need to violate the null energy condition to have a bounce (since dH 2 dρ = 1/(3M 2 Pl ) > 0 in that case). In our case, the derivative near the bounce can be easily computed from the Friedmann equation to give Thus, if we want to have a bounce corresponding to a transition from a contracting to an expanding phase for a matter component satisfying the null energy condition, we need This condition implies that the effective Newton's constant in the Friedmann equation for small H in Eq. (53) is negative. Then, from that equation we see that the bouncing solution at high energy density must match a de Sitter universe at low densities. In particular, this precludes the possibility of recovering the usual GR Friedmann equation at low densities and, as a consequence, realistic bouncing scenarios will be difficult to realize.
On the other hand, as we showed above, the stability of tensor perturbations requires 1 + . Thus, we need to have ξ = 0 for the bounce to be possible and avoid tensor instabilities. It is important to stress that these are necessary conditions, but not sufficient.
Another interesting scenario corresponds to the opposite case, i.e., a turnover in the cosmological evolution with a transition from an expanding to a contracting universe. This case corresponds to havingḢ < 0 near the turnover so that we need to require to guarantee the existence of these solutions. We can see that the transition from expanding to contracting universe is possible for models with ξ = 0, unlike in the previous case where ξ needed to be non-vanishing. Interestingly, this behaviour resembles the evolution in the presence of a non-vanishing curvature for the spatial sections, but for a flat FLRW metric. Thus, it is possible to mimic the effect of a non-vanishing curvature with the vector field.

B. Singularities
In the previous section we have studied the existence of turnonvers corresponding to transitions between contracting and expanding phases. However, we can also have another class of turnovers which actually lead to cosmological singularities. These cases correspond to points where dH/dρ becomes infinity for finite H and ρ. If we consider again the expression (58), we see that at those points there is a divergence inḢ, while H and ρ remain finite. Since the divergence only appears in the derivative of the Hubble expansion rate, it is plausible that such a singularity does not lead to a singular spacetime, i.e., the spacetime could be geodesically complete. In other words, the geodesics might go smoothly through the singular point. On the other hand, another worrisome feature of having a singularity is that, even if the geodesics can go through, tidal forces might diverge so that extended objects can not go through the singularity (see e.g. [83,84] for further discussions about these points). We will not discuss in further detail the specific properties of these singularities here, but we will content ourselves with analyzing their attracting properties. The points where singularity occurs can be easily computed from Eq. (52) as those values of H for which dρ/dH = 0. Then, near the turnover the Friedmann equation will read C. Self-tuning solutions In this subsection we will investigate whether the vector-tensor theories that we are considering can lead to selftuning solutions, i.e., Minkowski solutions in the presence of an arbitrary cosmological constant. For that we need to have H = 0 for ρ = 0. The idea behind the self-tuning mechanism is to make the vector field equation trivial for A 0 in Minkowski so that, then, it can be used in the Friedmann equation to screen an arbitrary cosmological constant. In our case, making the equation of A 0 in Minkowski is only possible if M 2 = λ = 0. However, for that case all the dependence on A 0 from Friedmann equation also drops in a Minkowski solution. This means that it is not possible to realize self-tuning solutions. This is not surprising and, in fact, in the original proposal within the class of Horndeski theories considered in [15,85] the presence of non-flat FLRW was crucial. In our case, we can show that the same conclusion can be reached. Let us consider a non-flat FLRW metric with k the curvature of the spacial sections. Then, instead of imposing H = 0 as it corresponds to pure Minkowski space, we impose the Ricci flat condition H 2 + k/a 2 = 0. The vector field equation then reads This equation is trivial only if M 2 = ξ = λ = 0. Then, the Friedmann equation reduces to and we see that A 0 can dynamically compensate for the cosmological constant analogously to the original proposal. Of course, this does not come completely as a surprise since the theory with all the parameters vanishing but β resembles one of the terms identified in [15,85] to generate self-tuning solutions.

V. MODELS EXAMPLES
Now that we have studied some general properties of the theory under consideration, in this section we will explore some specific models where the general features discussed above can be illustrated by simple cases which can be treated analytically. Among the considered examples, we will study the singular cases mentioned above.
A. Case λ = ξ = 0: de Sitter self-tuning A singular case is the one with λ = ξ = 0, since in that case the non-trivial branch of the vector field equation determines the Hubble function to be H 2 = M 2 /(6β), which further implies that M 2 and β must have the same sign. The value of A 0 will then be determined by the Friedmann equation with ρ * ≡ M 2 Pl M 2 /(2β). We see that this solution features the existence of a lower or an upper bound for the energy density for positive and negative β respectively. Thus, in this case the vector field can compensate for the matter density in order to maintain a de Sitter universe. However, the presence of a bounded allowed region for the density of the matter component prevents to screen an arbitrarily large or small (depending on the sign of β) density. This observation leads to the interesting possibility of a transition between branches. If β > 0, then above the minimal density ρ * we have a de Sitter phase with A 0 compensating for the matter source. When we reach the minimal density, A 0 transits to the trivial branch with A 0 = 0, while recovering the usual Friedmann equation 3M 2 Pl H 2 = ρ. However, the transition between branches is not smooth because although it is continuous, a discontinuity in the derivatives will appear (see Fig. 3). In this case the vector field equation becomes linear so we only have one branch (besides the trivial one) with As explained above, we can further fix ξ = 1 without loss of generality (or, equivalently, rescaling A 0 → ξ 1/3 A 0 and M 2 → ξ 4/3 M 2 ). Friedmann equation then gives the following two branches: If M 2 > 0, we see again the existence of a minimal density given by ρ ⋆ = M 3 M Pl / √ 6. However, if M 2 < 0, the energy density can take any value. In the limit of small densities we then obtain On the other hand, for high densities we find We then see that the branch H + has a transition from the usual GR behaviour at high densities to a de Sitter universe at low densities. Notice that such a solution is only possible for M 2 < 0, which also guarantees that the effective cosmological constant is positive. The solutions are illustrated for both negative and positive M 2 in Fig. 4.
C. Case λ = 0: modified early cosmology The vector field equation yields in this case We see that for M 2 > 0 there is a minimal allowed value for the energy density for both branches, while for negative M 2 we have the transition from the usual GR regime at high densities to a de Sitter phase at low densities, as discussed in the main text.
Again, without loss of generality, we will set ξ = 1. From the Friedmann equation we obtain The second term in the LHS going as M 4 /H 2 of this effective Friedman equation indicates that only if M 2 = 0 we can recover the usual GR equation in the small-H regime (we disregard the possibility β = 0 as we would reduce to the case considered in the previous section). Thus we will set M 2 = 0. Then, we can rewrite the above expression as Thus, if β < 0, there is a maximum energy density given by ρ ⋆ = 3M 4 Pl /(4|β| 3 ). Also notice that in the case with β < 0, we automatically have that H 2 is positive. Interestingly, we see that when such a maximum density is reached, the Hubble rate remains finite in both branches: H 2 ± → −M 2 Pl /β 3 , although the derivative diverges |Ḣ ± | → ∞, in accordance with our general discussion on the singularities. At low energy densities ρ ≪ ρ ⋆ , we find the usual GR result in the branch H − , while the branch H + gives a de Sitter universe with H 2 + = −M 2 Pl /β 3 irrespective of the value of ρ. In fact, these behaviours are effective up to energy densities very close to ρ ⋆ (see Fig. 5). For the de Sitter branch at low densities, we can use our results on the perturbations around the de Sitter fixed points in vacuum.
When β > 0, only the branch H − exists, since H 2 + is negative in that case. We find again two regimes determined by ρ ⋆ . At low densities ρ ≪ ρ ⋆ , we recover the GR regime with the usual Friedman equation, i.e., 3M 2 Pl H 2 − ≃ ρ. However, in the high energy density regime we have This behaviour leads to an evolution of the scale factor a ∝ t 4 3(1+w) so that we have accelerated expansion for w < 1/3. In particular, for dust we have that the first slow roll parameter is ǫ = −Ḣ/H 2 = 3/4. In principle, this would be in tension with CMB observations, but we should remember that since we are not in the GR regime, scalar perturbations do not necessarily acquire a spectral index determined by ǫ. Another interesting feature of this regime is that the effective Friedmann equation does not contain any dimensionfull constant since the dependence on the Planck mass of the full equation drops in this limit. The left panel shows the case with negative β and the right panel shows the case with positive β. We see the existence of an upper energy density for β < 0. In that case both branches exist, one mimicking GR in the range of allowed energy density while in the second branch we find a de Sitter universe irrespective of the value of ρ. In the case with positive β only one branch exists, reproducing GR at low energy densities, while at high energy densities we have H 2 ∝ √ ρ.
background, with the Λ given by the constant vector field as Obtaining consistently an effective positive cosmological constant requires M 2 < 0 and λ < 0. The vector field fluctuates, and thus model could be distinguished, at least in principle, from ΛCDM by studying perturbations.
In the asymptotic future de Sitter limit, both tensor and vector perturbations acquire the standard form for their quadratic lagrangians with kinetic term K v,t = 1 propagation speed c 2 v,t = 1 as consequence of having β = 0. For the scalar sector, the kinetic and potential coefficients defined in (45,46) are given by Since the kinetic coefficient is positive definite when Λ λM > 0 (implying M 2 < 0), ghost-like instabilities are avoided. The propagation speed vanishes identically so no Laplacian instabilities will arise neither. However, this is a very singular case and means that the scalar perturbations will not actually propagate at this order and it would be necessary to include some higher order terms.

VI. DISCUSSION AND CONCLUSIONS
In this paper we have studied cosmological solutions for a class of vector-tensor theories that arise within the framework of the geometries introduced in [40,42]. Those geometries are defined by the most general connection linearly determined by a vector field (and without derivatives). We have then considered the most general action quadratic in the Riemann curvature tensor for this connection. Remarkably, we found that the absence of ghostly degrees of freedom associated to more than 3 polarizations for the vector field restricts the connection to be of the extended Weyl type. This extended version of Weyl geometry maintains the conformal invariance of the non-metric compatibility while, in addition to non-metricity, allowing the presence of the vector trace torsion. The resulting vector-tensor theory contains the usual Maxwell kinetic term with a mass and a quartic potential. Moreover, we also found a cubic derivative self-interaction for the vector field which gives the cubic Galileon in the decoupling limit and a non-minimal coupling to the Einstein tensor. Although in this work we have focused on the quadratic theory, we can argue that higher order derivative self-interactions with Horndeski vector-tensor and Galileon-like terms for the Stückelberg field in the decoupling limit can also be constructed by considering cubic and higher order curvature terms in the action.
For the vector-tensor theory arising from quadratic curvature terms, we have shown the existence of isotropic de Sitter solutions. These solutions are supported by the temporal component of the vector field, which is non-dynamical and plays the role of an auxiliary field. For the general case, we have obtained the regions in the parameter space of the theory where such solutions exist and performed a classification according to the number of de Sitter branches. We have then studied the stability of the inhomogeneous perturbations around the de Sitter critical points by computing the quadratic action for the scalar, vector and tensor perturbations. This analysis explicitly showed that the theory only propagates the 2 polarizations of the graviton plus the 3 polarizations of the massive vector field. For vector perturbations the only effect is a modification of the propagation speed so that they are never ghostly and only a potential Laplacian instability can be present. Furthermore, the effects in both the vector and the tensor sectors vanish if the coupling to the Einstein tensor is turned off. For the scalar sector there is a modification of the kinetic term and the propagation speed that depends on all the terms in the action, so the scalar mode can potentially contain both ghost and/or Laplacian instabilities around the de Sitter solutions.
After devoting a careful analysis to the de Sitter solutions we have turned to the general cosmological evolution in the presence of a matter sector. We have shown that, typically, we might expect to have upper and lower bounds for both the Hubble expansion rate and the energy density. We have also considered the possibility of having bouncing solutions and found that this is only potentially possible in a stable way if the vector Galileon interaction is present. In any case, we have shown that the effective Newton's constant in the background Friedmann equation must be negative at the bounce, although this by itself does not represent an instability. However, we ahve shown that bouncing solutions at high energies generically cannot match continuously a GR regime at low densities. Moreover, we have also studied the presence of turnovers with a transition from an expanding to a contracting universe. Another type of turnovers that we have analyzed is those happening at a finite H so that they cannot correspond to the discussed transtions. We have shown that these additional turnovers correspond to a cosmological evolution whereḢ diverges. We have also analyzed the existence of self-tuning solutions such that we can have Minkowski solutions in the presence of an arbitrary cosmological constant. As in the Horndeski class of theories, this is only possible for non-flat FLRW universes for our vector-tensor theory. Moreover, such solutions can be realized if only the coupling to the Einstein tensor is present, with all the other terms vanishing.
Finally, we have considered some specific examples to illustrate the general results discussed in the first part of the paper. We have shown that the class of vector-tensor theories studied in this work and which naturally arise in the framework of geometries with a linear vector distortion can give rise to a rich an interesting phenomenology for cosmology. The effects of the vector field can be important both in the early universe or at late times depending on the values of parameters. Therefore, they can be used to build dark energy/dark matter models or inflationary scenarios. Moreover, the natural presence of upper bounds for the energy density and/or the curvature might help evading singularities without resorting to quantum effects. In this work we have focused on purely isotropic universes, but the spatial components of the vector field and the class of anisotropic solutions that may be obtained could also deserve attention.
We can immediately see that these coefficients identically vanish for the generalized Weyl geometries with b 3 = 2b 1 −b 2 . One would expect that, given the numerous independent coefficients available, additional solutions to the equations γ i = 0 exist for any D. In the following we will show that this is not the case. Let us assume that b 3 = 2b 1 − b 2 and solve γ 3 = 0 for c 6 to obtain c 6 = 1 − 2c 3 − 2c 4 − c 5 . If we plug this solution into γ 1 and γ 2 we obtain: Now we can solve γ 1 = 0 for any of the remaining parameters, say d 2 , and insert it in the expression for γ 2 which finally reduces to Thus, we see that for D ≥ 4 there are no additional solutions besides the generalised Weyl geometry with b 3 = 2b 1 −b 2 . We can also see that in lower dimensions additional solutions exist.