Dark energy with a shift-symmetric scalar field: obstacles, loophole hunting and dead ends

We discuss the possibility of a scalar field being the fundamental description of dark energy without the need of self-tuning any parameter. Thereby, we focus on shift-symmetric scalar-tensor theories satisfying also that the propagation speed of gravitational waves is equal to the speed of light. Analysing the stability of scalar linear perturbations, we discuss the conditions that seems to be necessary to describe (super) accelerated cosmic expansion without introducing instabilities. It has been established, however, that this stability can be ruined when taking into account the interaction with tensor perturbations (essentially gravitational waves). Indeed, although we shall point out that the standard proof of absence of dark energy stable braiding models due to this interaction has a possible way-out, we find general arguments suggesting that there are no dark energy stable solutions that can exploit this loophole. Thus, we discuss future research directions for finding viable fundamental descriptions of dark energy. We also provide a dictionary between the covariant version of the scalar field theory and that of the E ff ective Field Theory approach, explicitly computing the parameters in the latter formalism in terms of the functions appearing in the covariant version, and its derivatives. To the best of our knowledge, this is the first time these expressions are explicitly obtained up-to arbitrary order in perturbation theory.


Introduction
Scalar field theories could be crucial in providing an underlying theoretical framework for addressing the mysteries of dark energy (DE).These theories introduce a scalar field that can be used to describe dynamical DE, in contrast to the cosmological constant case (see, for instance, references [1][2][3]).To select the appropriate scalar field theory for cosmology, it is reasonable to demand that the corresponding field equations to be of second order.Otherwise, the new degrees of freedom (DOF) may destabilise the dynamics of the theory by means of the Ostrogradski ghost [4].In this context, a natural framework is that provided by Horndeski theory [5] (see also [6] for a review), although some theories with higher-order derivatives in the action can also avoid the Ostrogradski instability through degeneracy conditions in the Lagrangian (see, for instance, [7,8]).Rather recently, the speed of propagation of gravitational waves (GWs) has also become an important consideration when formulating theories beyond general relativity.After the detection of the event GW170817 [9] and its electromagnetic counterpart [10], it has been concluded that the propagation speed of GWs is very close to that of light.Selecting the Horndeski models that trivially propagate GWs at the speed of light reduces the original parameter-space of the theory to the so-called viable subclass of Horndeski [6].This includes the well-known k-essence scenario [11][12][13][14][15], the more general Kinetic Gravity Braiding (KGB) theory [16], and a possible non-minimal coupling to gravity via the scalar field (see, for instances, reference [17]).Scalar fields minimally coupled to gravity and with a canonical kinetic term have been already used to describe the current accelerated expansion of the Universe.However, in that simple case one has to rely on the particular form of a potential term, re-introducing the fine-tuning problem usually associated with the cosmological constant.This issue can be avoided by focusing on models that are invariant under constant shifts in the scalar field, for which potential-like terms do not appear in the Lagrangian.This symmetry-based argument is an elegant way to make the evolution of the system to depend only on the rate of change of the scalar field but not on the scalar field itself (see discussion in reference [16]).However, the simplest minimally coupled scalar field with a canonical kinetic term in the action is only able to describe a stiff fluid when imposing shift-symmetry (see, for example, references [2,18]); leading to a non-interesting phenomenology for the (late-time) DE sector.On the other hand, in the Horndeski subfamily previously mentioned this argument prevents from having a non-minimal coupling to gravity via the scalar field itself, although a derivative coupling to matter is still allowed in principle.In either case, evading fifth-force constraints coming from astrophysical scales is naturally easy under this symmetry consideration [16].The Horndeski models that satisfy this symmetry condition and predict a propagation speed for GWs equal to the speed of light are known as shift-symmetric KGB theories [16].
Shift-symmetric KGB theories are well-known for producing self-tuning de-Sitter future solutions [16] (see also references [19][20][21][22][23][24][25]), but they can also describe different future phenomenology [26,27].Contrary to the k-essence scenario, in the KGB theory second order derivatives of the metric are present in the scalar field equation of motion and vice versa, second order derivatives of the scalar field appear in the metric field equation [16].As a result, these theories can describe a phantom regime that is stable in the sense of no ghost-or gradientlike instabilities for linear scalar perturbation [16].In fact, the KGB theory can also produce a smooth phantom crossing in the DE sector [16], a feature not possible with a single scalar field à la k-essence [28,29].With regard to phantom energy, it should be highlighted that the violation of the Null Energy Condition (NEC) in the DE sector is not only still allowed by the observational data (see discussion in reference [30]), but it has been indeed proven that a phantom regime (with phantom crossing) is a necessary prerequisite to ease both the H 0 and S 8 cosmological tensions simultaneously by taking into account new physics that is relevant only at late cosmic times [31,32].Therefore, the motivation, simplicity and apparent stability properties of the shift-symmetric KGB scalar field models makes this framework an interesting proposal for modelling DE.Nevertheless, the absence of instabilities in linear scalar perturbations is a necessary but not a sufficient condition for the stability of the classical theory.Indeed, it has been shown that the interaction mediated by the braiding term between tensor perturbations (essentially GWs) and DE fluctuations may induce a ghost-like and/or gradient-like instability in the scalar sector [33].(Se also references [34,35] for a discussion on the decay of GWs into scalar fluctuations when Lorentz invariance is spontaneously broken.)Consequently, it was concluded that the braiding term in the KGB theory should not produce any sizeable effect in order to trivially escape from the GWs-induced instabilities [33].According to these results, the available parameter-space of the shift-symmetric KGB models seems to shrink to that of the shift-symmetric k-essence theory only [33].
In the present work, we first re-analyse the inter-relation between phantom behaviour and classical and semi-classical instabilities of the scalar field in the framework of shiftsymmetric KGB models.By implementing a novel characterization of the already-known stability conditions, we shed some light on previous results in the literature for k-essence models [36] (see also, for instance, references [11,12,[37][38][39]) and also obtain new compact results for the braiding term.In the second place, we reflect about the instability induced by the interaction of the braiding term with GWs [33].We propose a possible way to circumvent this obstacle when departing from the Galileon terms for the KGB functions.Unfortunately, we shall argue that the apparent loophole seems to lack practical appli-cability for constructing viable dark energy models with a nonvanishing braiding term.In order to perform this study, we have derived the detailed form of the mass parameters of the KGB theory in the Effective Field Theory (EFT) approach.Although the expressions for these parameters were already computed at leading order, to the best of our knowledge, this is the first time that the complete family of parameters up-to arbitrary order is presented.
This work is organised as follows: we review the general properties of shift-symmetric KGB theories in the first part of section 2. The stability of linear scalar perturbations in the wellknown k-essence subclass of the more general KGB framework is reviewed in section 2.1.In sections 2.2 and 2.3 we discuss the viability of linearised scalar perturbation for the remaining part of the KGB theory.Section 3 is devoted to the interaction of the braiding term with GWs.We review the arguments presented in reference [33] in section 3.1.We identify a potential mechanism to ease the constraints found in reference [33] in section 3.2.In section 3.3 we discuss that the possibility for a viable KGB model to implement such a mechanism is an elusive ambition in practice.Therefore, in section 4, we reflect about the open paths that can still be further explored to find fundamental viable descriptions for dark energy.Finally, auxiliary results for discussing the stability of linear perturbations are summarised in Appendix A. A self-contained guide to the EFT approach to DE and useful calculations are presented in Appendix B.
It should be noted that one of the arguments most commonly used to go beyond the standard model of cosmology is the existence of a fine-tuning problem regarding the value of the cosmological constant.However, having potential terms in a scalartensor theory could reintroduce the fine-tuning problem now regarding the parameters that fix the minimum of the potential, if the accelerated expansion is reached in the potential domination regime.One possible way of trying to avoid this problem is the consideration of a shift-symmetric scalar field, which guarantee that the cosmic acceleration is driven by the kinetic energy.Those are models that are invariant under constant shifts in the scalar field, for which potential-like terms do not appear in the Lagrangian.This symmetry-based argument is an elegant way to make the evolution of the system to depend only on the rate of change of the scalar field but not on the scalar field itself (see discussion in reference [16]).Therefore, in this work we will focus on the shift-symmetric sector of the KGB theory [16] (see also [19][20][21][22][23][24][25][26][27]).This is when the above action is invariant under constant shifts in the scalar field given by being c some constant.In practice, this implies that the functions K and G do not depend on the scalar field but only on its kinetic term.Hence, we will focus on the particular case of K = K(X) and G = G(X).(Nevertheless, this assumption will be relaxed when discussing the EFT approach in Appendix B.) In addition, we will also consider the spatially flat cosmological background described by the homogeneous and isotropic Friedmann-Lemaître-Robertson-Walker (FLRW) space-time.That corresponds to the line element where a(t) stands for the scale factor and dx 2 3 are the spatial three-dimensional flat sections.Taking matter and radiation as external sources to the action (1), the field equations of the theory for this background read [16] where a dot represents derivation w.r.t. the cosmic time t and G X stands for dG/dX.In addition, is the shift-current related to the global shift-symmetry of the theory [16].In view of these equations, the energy density and pressure of the dark fluid can be defined as and respectively.Then, the equation of state parameter and the partial energy density for the scalar field take the form w ϕ p ϕ /ρ ϕ and Ω ϕ ρ ϕ /(3H 2 ), respectively.Contrary to k-essence [11], the energy density, ρ ϕ , now depends on the Hubble parameter through J while the pressure, p ϕ , contains φ [16].Moreover, the pressure will also show an explicit dependence on the Hubble function (and the external sources) when φ is removed using the field equations.
The field equation of the scalar field can be obtained varying the action (1) w.r.t. the scalar field.This yields where has been defined for the compactness of the notation.This equation is completely equivalent to the covariant conservation of the ϕ-fluid in the hydrodynamic approach, that is where ρ ϕ and p ϕ correspond to the scalar field energy density and pressure defined in expressions ( 9) and (10), respectively.It should be noted, however, that the scalar field equation of motion (either in the form of equation (11) or, equivalently, in the hydrodynamic approach depicted in equation ( 13)) is not independent from the system of equations given in expressions (4) to (7).
Alternatively to the scalar field equation (11), one can obtain the same information from the conservation equation related to the shift-symmetry of the theory; that is As noted in the reference [16], this conservation trivially leads to a first integral of motion for the system given by where Q 0 is a constant.This feature of the shift-symmetric KGB theories can be used to significantly simplify the discussion on the evolution of these cosmological models [16].This is because J = 0 defines a surface in the configuration space towards which all trajectories evolve if the scale factor is growing large.Therefore, the locus J = 0 can be used to simplify the field equations, at least asymptotically.However, this simplification should be taken with due caution: the vanishing of the current J may not always imply that the term φJ should be dropped from equations ( 4) and (5); see discussion in references [26,27].In general, the future phenomenology of this shift-symmetric scalar field theory has been shown to encompass a great variety of options.From self-tuning de-Sitter-like future solutions [16] (see also references [19][20][21][22][23][24][25]) to DE-driven cosmological singularities [26,27].The stability of scalar perturbations around a FLRW background have already been addressed at linear order for this scalar field theories [16] (see also references [24,41,42]).Indeed, the quadratic action describing scalar perturbations reads [16,24] where ζ represents the curvature perturbation as commonly defined (see, for example, reference [16]), Q S is the normalization factor of the kinetic term, and c s denotes the speed of propagation (hereon referred to as the "speed of sound") of the scalar perturbation.From the above action, the absence of ghost and gradient instabilities in the scalar sector implies [16,41] respectively, where the dimensionless functions α K and α B were first introduced in reference [41].They are defined as Note that condition (17) trivially reduces to D > 0, since the denominator of Q S is always positive.Violation of this condition leads to scalar perturbations whose associated kinetic energy takes negative values [43]; i.e. the perturbations become ghost-like.The presence of a ghost DOF may be detrimental (at the classical level) if the ghost interacts with a positiveenergy mode, as this may lead to runaway solutions where the total energy is conserved but the relative energies associated to the ghost and no-ghost sectors diverge (see, for instance, reference [44]).At the quantum level, the corresponding ϕparticles would have negative energies.In this case the vacuum becomes unstable due to the spontaneous production of ghost-particles together with normal-particles with arbitrarily high energies and momenta [44] (see also reference [43]).Consequently, backgrounds with ghosts are generally considered pathological at both classical and quantum levels.On the other hand, violation of condition (18) introduces the so-called gradient instability.This is when the leading order spatial derivatives have the wrong sign w.r.t. the time derivatives in the perturbed action (16) (see also references [43,45]).In Fourier space, the frequencies of the oscillations become imaginary at high momenta, resulting in perturbations that grow exponentially fast.Thus, it precludes a stable classical model (at least, as long as a perturbative treatment is still valid).
It is interesting to note that and, therefore, the ghost-free condition D > 0 (see equation ( 17)) implies in our notation.Moreover, replacing the time derivative of the Hubble rate that appears in the scalar field equation ( 11) by means of the Raychaudhuri equation ( 5) leads to Hence, the ghost parameter D is just the factor ahead of the scalar field acceleration.
In the following we will analyse the stability conditions described in equations ( 17) and (18) for three different scenarios: (i) the k-essence model (G ≡ 0), (ii) when K ≡ 0, and (iii) the most general KGB scenario.These are sections 2.1, 2.2 and 2.3, respectively.In doing so we will review some well-known facts and also face new discussions about the viability of these shift-symmetric theories from the point of view of linear cosmological perturbations.

K-essence
Let us first re-analyse the stability for the well-known kessence case [11][12][13][14][15].This corresponds to setting the braiding function G to zero in action (1), although a non-zero but constant braiding function would also reduce to the k-essence scenario.From this point of view, one may consider the k-essence theory as a special sub-case of the more general KGB action (1).However, it is worth noting that the two scalar field theories have very different properties.This is because second order derivatives of the metric and scalar field are no longer mixed when G is constant.Thus, some of the most interesting features of the KGB set-up (such as stable phantom crossing) are not present in the k-essence scenario.Still, k-essence represents the most general scalar theory whose action contains up-to first order derivatives of the scalar field and has been extensively explored in a wide variety of cosmological scenarios, from inflation [11,37] to DE [12][13][14][15].
Unlike the scalar field with a canonical kinetic term, for a k-essence field a phantom behaviour is not necessarily due to the presence of a ghost [39] (see also [36]), since equation (17) reduces to whereas the violation of the NEC by the scalar field depends on the slope of the k-essence function K, being cf. with expressions (9) and (10) when G X is trivial.The slope of the function K should be negative for a phantom field, but this do not necessarily force the energy density ρ ϕ to be a decreasing function of the kinetic term.So, the sign on the r.h.s. of equation ( 24) is in principle independent of that of equation (25); see reference [39].However, the relation between instabilities and phantom behaviour appears when considering the squared speed of sound of the perturbations, i.e. c 2 s .Taking into account the Friedmann equation ( 5), we can recast the gradientfree condition (18) in term of the background quantities as where Ω ϕ ρ ϕ /3H 2 is the partial energy density for the scalar field and we have also used the fact that α B ≡ 0 for k-essence.This new reformulation of the expression for c 2 s is extremely useful since we can address the stability of linear perturbations in terms of the background quantities for which we have a strong physical intuition.Since we are focusing on the case Ω ϕ > 0, this inequality shows an inter-relation between the ghost-and gradient-free conditions with the possible phantom character of the ϕ-fluid.In particular, if the NEC is satisfied for the scalar field (i.e. if 1 + w ϕ > 0), there is not a gradient instability if and only if the ghost-free condition is also satisfied.As a result, there are k-essence models with a completely stable scalar spectrum.On the other hand, the gradient-free and ghost-free conditions are anti-correlated when the NEC is violated for the ϕ-fluid.That is to say, c 2 s is positive for a phantom fluid only if a ghost-like instability is present and vice-versa, scalar perturbations are ghost-free only if a gradient instability is triggered.Please note that the same conclusion was reached in reference [36], albeit through a different reasoning (see also references [39,46]).
As is well-known, it is not possible to describe an effective phantom fluid via a scalar field à la k-essence in a way that both ghost and gradient instabilities are absent [36,39,43,46].However, let us reflect about the nature of these instabilities (see, for example, reference [43] for a review on the topic).On the one hand, the presence of gradient instabilities is related to a wrong sign for the spatial derivatives in the perturbed action (16).This makes the frequencies of the fluctuations imaginary at high momenta, resulting in perturbations that grow arbitrary fast.Thus, it precludes a stable classical model (at least, as long as a perturbative treatment is still valid).For a ghost mode, on the other hand, the associated kinetic energy is negative.If the ghost DOF interacts with a positive energy mode (at the classical level) there may be runaway solutions, where the total energy is conserved while individual energies diverge; we referred the interested reader to the discussion in, for example, reference [44].Thereby, the common lore states that ghost instabilities are catastrophic (already) at classical level.Still, the classical background is stable against high momenta perturbations [43].Nevertheless, it is interesting to note that there are examples in the literature in which the presence of a ghost DOF interacting with a positive energy DOF do not lead to runaway solutions and, therefore, to the destabilization of the classical motion of the system [47].(See also "islands of stability" [48,49] and meta-stability of ghosts [50,51].)Upon canonical quantisation of theories with ghosts, the energy conservation does not forbid pair creation from vacuum of ghosts-particles together with normal-particles: the vacuum becomes quantum mechanically unstable (see, for instance, references [43,44]).So, the presence of a ghost is also considered as pathological at the semi-classical level.However, the possibility of safely living with ghosts at quantum level has already been explored with positive conclusions [52][53][54][55] (see also references [56,57] for applications to quadratic gravity).Therefore, even if it is clear that the presence of any of these instabilities (ghost and/or gradient) is a red light that has to be seriously taken into account, it would be interesting to consider the possibility that the scalar field (understood as being of gravitational nature) should be quantized in a different way to matter fields (for example, with the space-time itself [58,59]) or not been quantized at all, when dealing with the ghost instability.For a review on the topic see, for example, reference [60].
In this context, a gradient instability is hardly tameable in the classical theory, whereas a ghost-like instability could be further discussed at both classical and quantum levels.So, if one considers that the H 0 tension has to be solved by means of a phantom fluid, then, our best hope in the k-essence case would be selecting the function K in such a way that only the ghost-free condition is violated.Otherwise, a negative c 2 s will inevitably jeopardise the perturbative approximation due to fast growing solutions.On the other hand, it is also worth to mention that the scalar field in (shift-symmetric or not) k-essence theory cannot cross the phantom divide [28] (see also [29]).The equation of state parameter for the scalar field, namely w ϕ , is always greater or less than −1 for a single scalar field à la kessence.Nevertheless, a phantom-crossing in the DE sector was analytically shown to be a prerequisite for solving both the H 0 and S 8 tensions simultaneously [31,32].(Similar conclusions were also reached by addressing the resolution of the H 0 tension only [61,62].)Therefore, if both tensions were considered on the same footing1 , a scalar field theory beyond k-essence should be explored.See also references [61,62] for previous works on the H 0 tension.

Braiding model
We now analyse the case of K ≡ 0 in action (1).This is a straightforward limiting scenario of the complete KGB framework where the contribution of the k-essence part is negligible.However, contrary to pure k-essence discussed in the previous section, second order derivatives of the metric and scalar field are still mixed in the corresponding field equations due to the function G; see equations ( 5) and (11).As a result, this is a scenario simpler than the complete KGB theory that still provides us with an interesting candidate for addressing the stability properties of the kinetic braiding.
The expressions for the energy density and pressure of the scalar field in equations ( 9) and (10), respectively, simplify to where the expression for J in equation ( 8) has been also considered.Then, the equation of state parameter for the scalar field reads The scalar field acceleration, φ, can be removed from the previous expression by means of the Raychaudhuri and scalar field equations, these are expressions ( 5) and ( 11), respectively.Then, after simple manipulations we can re-write the scalar field equation of state as where Ω r ρ r /3H 2 is the partial energy density for radiation and the contribution of matter has been removed through the The lower oblique-lined region corresponds to the presence of a the ghost instability.Scalar perturbations are both ghost-and gradient-unstable in the central vertical-lined region.Stability is ensured in the w eff > 1/3 zone, since w ϕ is clearly positive there in virtue of the relation w eff = Ω r /3 + w ϕ Ω ϕ .In the dotted region, which corresponds the horizontal band 0 ≤ w eff ≤ 1/3, perturbations may be healthy if and only if w ϕ is non-negative, otherwise they will be both ghost-and gradient-unstable.Nevertheless, the sign of w ϕ in that region of the (Ω ϕ , w eff )-plane cannot be inferred without solving the whole evolution of the model and, therefore, the stability conditions there cannot be further addressed without specifying the function G.
constraint Ω m = 1 − Ω r − Ω ϕ coming form the Friedmann equation (4).Please note that the above expression allow us to reexpress the ghost-free condition (17) in terms of background quantities as As emphasised in the previous section, this new approach to the evaluation of the stability conditions is extremely useful since we have a strong physical intuition for the values of background quantities.The functional form of the equation of state for the scalar field, that is to say w ϕ = w ϕ (a, ȧ, φ, φ), will depend on the choice for the braiding function G.However, if we assume a positive Ω ϕ , it follows from the above inequality that a ghost instability is always present if w ϕ is negative.A significant implication of this result is that any model that aims to describe a DE component (not necessary phantom-like) will inevitably have this instability, since w ϕ should be less than −1/3.A similar derivation can also be applied to the expression for the speed of sound squared in equation (18).Straightforward differentiation over the definition for α B in equation (20) leads to where w eff p tot /ρ tot represents the effective equation of state parameter of the total fluid on the r.h.s. of equations ( 4) and (5).Taking this result into account, and also substituting Ḣ by means of equation ( 5), the gradient-free condition ( 18) can be re-formulated in terms of background densities as The numerator in the above expression can be positive or negative depending on the value of Ω ϕ and w eff .In Figure 1, we represent the fulfilment of the stability conditions in the (Ω ϕ , w eff )plane.As there can be checked, it is possible to have ghost-free and gradient-free scalar perturbations during radiation and matter dominated epoch; that would correspond to w eff ≈ 1/3 and w eff ≈ 0 for Ω ϕ ≪ 1, respectively.A necessary and sufficient condition for this to happen is that of w ϕ being positive; see inequality (33).On the contrary, perturbations will always become unstable at some point if w ϕ is negative; for example, if the braiding scalar field acts as DE.Hence, an accelerating universe cannot be safely modelled when K ≡ 0 in the KGB action (1).(Recall that acceleration demands the violation of the Strong Energy Condition and, therefore, w eff < −1/3 at some moment in the recent past.According to Figure 1, this would inevitably produce ghost and gradient unstable scalar perturbations.)Interestingly, c 2 s could be positive for a super-accelerated regime (w eff < −1), although this would only be possible at the price of having a ghost mode.This feature may be interesting for modelling phantom DE since, as we have discussed in the previous section, the positivity of the speed of sound squared can be argued to be more fundamental from a classical point of view than the absence of a ghost.However, from Figure 1 it is clearly not possible to have a braiding model (with K ≡ 0) that connects this gradient-stable supper-accelerated (late-time) regime with a realistic early-universe description in a way that c 2 s remains always positive.Therefore, modelling DE with only the braiding function G of the KGB theory seems rather unwisely.

The complete KGB theory
In the previous sections we followed a novel approach to the evaluation of the stability of linear scalar perturbations by rewriting the ghost-free (17) and gradient-free (18) conditions in terms of background quantities, like the partial energy densities Ω i and the equation of state parameters w ϕ and w eff .Our main purpose with this approach is to have a better physical intuition about the requirements for a given model to produce stable scalar perturbations.From this analysis, we found that instabilities are always triggered when 1 + w ϕ < 0 (violation of the NEC in the DE sector) in the case of k-essence.Moreover, modelling kinetic braiding only with the function G surprisingly makes things even worse, since now instabilities appear as w ϕ becomes negative.Thus, no stable DE component (for which w ϕ should be less than -1/3) can be produced in the case of K ≡ 0. In this section, we will apply the very same analysis now for the full KGB theory.
When both K and G functions are non-trivial, the ghost-free condition ( 17) cannot be re-expressed solely in terms of the previous background quantities but where are the partial energy densities for the contributions of the shiftcurrent, J, and that of the braiding function, G, to the different terms in the expression for scalar field energy density ρ ϕ in equation ( 9), respectively.Please note that Ω B J coincide with the definition of α B [41], see equation (20).It is also important to highlight that the physical constraint of Ω ϕ ∈ [0, 1] does not imply, in general, any bound on the value of these auxiliary partial energy densities.Their value can be arbitrary large and/or negative.
Similarly, the gradient-free condition ( 18) can be translated to a relevant quantity for the discussion in the next section.Please note that when G or K are trivial, the above expressions reduce to those obtained in sections 2.1 and 2.2, respectively.Contrary to the limiting cases analysed in the previous sections, conditions (36) and ( 39) can be satisfied simultaneously for a wide variety of models including those with phantom behaviour.For instance, considering η to be positive (since this will be crucial in the incoming sections) we show all the possible cases where both inequalities are satisfied in Appendix A. Nevertheless, since Ω B J and Ω J are not physically meaningful quantities as Ω ϕ is, it is difficult to extract the most general conditions the functions K and G should meet in order to produce a stable model.Anyway, there are examples of such models in the literature.For example, in reference [16] the authors analyse the case of K and G being proportional to X.They also show there that the resulting cosmological model has completely stable linear scalar perturbations even though the scalar field is violating the NEC.(Note, however, that their model does not feature super-acceleration since 1 + w eff > 0 always.)More examples of stable KGB models can be found in references [24,[68][69][70].Hence, it is possible to construct reliable cosmological models with the scalar field coming from shift-symmetric KGB theories when both k-essence and braiding functions are non-trivial.

Summary
Could we describe a stable DE component with a scalar field coming from shift-symmetric KGB theories?There is no problem to describe it with a k-essence term if it satisfies the NEC.However, it is well-known that k-essence suffers from instabilities if the NEC is violated [36,39,43,46].Moreover, the scalar field energy density should be always phantom or non-phantom since k-essence cannot produce a phantom-crossing [28] (see also [29]).Hence, if we interpret current data as favouring the presence of a slight violation of the NEC in the DE sector at present times, then, the scalar field should have always violated the NEC.In this scenario, a k-essence proposal will be plagued with a gradient or a ghost instability.With an appropriate selection for the function K, however, the gradient-free condition can be satisfied for (even) NEC-violating models.Nevertheless, the resulting scalar field will feature a ghost DOF; yet this may not always be catastrophic.As we have discussed, it would be interesting to explore whether it is possible to have a ghost mode in cosmology without destabilising the classical nor the semiclassical regimes of the theory.While this may be feasible, addressing both the H 0 and S 8 tensions (in case the latter is still present in future surveys) by means of late-time physics only continues to pose an insurmountable obstacle in the k-essence scenario, since no phantom-crossing is possible (see discussion in references [31,32]).
The situation is significantly worsen when only the braiding function G is present in the scalar field's action (1).Scalar perturbations are always ghostly if the scalar field acts as a DE component, i.e. if w ϕ < −1/3 (see condition (33)).This is aggravated even further by the fact that a gradient instability is also present in any realistic model that interpolates between an initial radiation dominated epoch and a late-time accelerated expansion regime; see Figure 1.Hence, the term G2ϕ alone in the action (1) is not a viable option for modelling DE.
Finally, when considering the sum of k-essence and the braiding term the situation is different.Following the same strategy we have used to analyse the preceding marginal examples, we have recovered the well-known result that the conditions for the avoidance of ghost and gradient instabilities can be simultaneously fulfilled even when the ϕ-fluid violates the NEC (see Appendix A).Although the lack of a physical intuition for the auxiliary quantities Ω J and Ω B J did not allow us to restrict the possible phenomenology of the theory as stringently as for the previous cases, note that linearly stable KGB models have been already reported in the literature (see for instance [24,[68][69][70] and references therein).This fact points to the viability (from the point of view of the stability of linear scalar perturbations) of modelling a realistic DE component with the scalar field from the (complete) KGB theory, at least in principle.

Beyond scalar perturbations
In section 2 we have discussed the absence of ghost and gradient instabilities in the linearised scalar perturbations around a flat FLRW background.In that approach, interactions between scalar, vector and tensor perturbations were not present, since they decoupled at the linear level.However, these interactions can appear when going beyond the linear regime and may be important.In fact, the interaction between tensor perturbations (essentially GWs) and DE fluctuations may induce ghost and/or gradient instabilities in the scalar sector, as it was pointed out in reference [33].In that case, the stability conditions discussed in the previous sections do not guarantee the viability of the theory 2 , since the braiding term in the action can destabilize the theory through the GWs-DE interaction.As a concrete example, for cubic Galileon (which corresponds to G(X) ∝ X) this effect is detrimental when |α B | ≥ 10 −2 ; that is to say, if the effects of the braiding are sizeable [33].
In general, the result of reference [33] raises serious doubts on the viability of any KGB theory (different from k-essence) aiming to describe the currently observed abundance of DE in the Universe.Nevertheless, as the authors of that reference commented (see also reference [71]), it is also important to analyse the fate of the instabilities once they are originated.The fate of the theory after the instability is reached is, however, uncertain since it depends on an unknown UV-completion.
In the following, we review the arguments and notation presented in reference [33]; that is section 3.1.In section 3.2 we obtain the general conditions a KGB model should meet in order to not be destabilised when interacting with GWs, even when the contribution from the braiding is non-negligible.The latter possibility turns to be attainable only if the braiding function G is not linear in the kinetic term X (i.e. if it is different from that of cubic Galileon).In section 3.3 we explore whether these conditions could be satisfied by some reasonable KGB model.

Interaction with gravitational waves: a review
In reference [33] the authors use the EFT approach to discuss the interaction between DE perturbations and GWs.(See Appendix B for a review on the EFT approach to DE.)For this effect to be sizeable one needs a cubic coupling ππγ in the perturbed action, where γ stands for tensor perturbations (essentially GWs) and π denotes scalar field fluctuations [33] (not to be confused with the scalar field ϕ itself).In the covariant version of the theory, the interaction ππγ is related to the Cubic Horndeski operator; i.e. the braiding term G2ϕ in action (1).Nevertheless, this interaction vertex is also present in beyond Horndeski theories [33][34][35].
The interaction between DE perturbations and GWs is conveniently addressed in the Newtonian gauge [33] where Φ and Ψ are the Newtonian potentials, and γ i j is transverse and traceless.At leading order, the Lagrangian density 2 See also references [34,35] for a discussion on the decay of GWs into scalar fluctuations when Lorentz invariance is spontaneously broken.describing this interaction reads 3 [33] where the expansion of the Universe has been neglected and indices are raised/lowered with the Minkowski metric.In addition, c 2 s is given by expression (18) and being D the ghost factor (17).Moreover, the parameter η is defined through the relation [33] 4 m3 where m3 1 and m3 2 are functions on the background time, t, that appear in the EFT action (B.15).This parameter is useful for measuring deviations with respect to cubic Galileon (i.e.G(X) ∝ X) for which η is zero and, therefore, 4 m3 2 = − m3 1 .(Further details can be found in Appendix B.) The covariant version of this parameter can be obtained from the expression for the mass parameters m3 1 and m3 2 in terms of the function G and its derivatives (see Appendix B).Taking the expressions (B.43) and (B.44) into account, it is straightforward to obtain that the covariant version of η is nothing but the quantity we already have defined in equation (40).That is η = 2XG XX /G X .
The field equation for the scalar field perturbation, π, follows from varying the action corresponding to the Lagrangian density (45).This reads [33] π 3 This Lagrangian density follows from restoring time diffeomorphism in the EFT action (B.15) for the KGB theory using the Stückelberg procedure.That is, by performing in the action (B.15) the transformation t → t + π(t, x i ), where π denotes the scalar field fluctuations.The relevant transformations in this case are [72] (see also [33]) where the former transformation is an exact relation to all order in π, whereas the latter has been expanded up-to second order in perturbations.After restoring full gauge invariance, the Newtonian gauge (41) can be implemented.The Lagrangian density (45) is finally obtained after solving (at linear order) for the Newtonian potentials in terms of π, and canonically normalizing the scalar and tensor fluctuations as [33] π → √ DHπ and respectively, being D the ghost function (17).
where ∇ 2 ≡ η i j ∂ i ∂ j .Here we will restrict our discussion to the case where π propagates subluminally, that is when c 2 s < 1 [33].(The luminal and superluminal cases are discussed in reference [33].)In this scenario, the lightcone for π is narrower than the one for the GWs, since c 2 GWs = 1 [33,35].Consequently, the DE fluctuations are not sensitive to the source of the GWs, provided we are far enough from the source of emission [33,35].To address the DE-GWs interaction, therefore, we can consider the following classical GW solution with linear polarization "+" travelling in the z direction [33,35] being ϵ + i j = diag(1, −1, 0) and h + 0 the dimensionless strain amplitude.Then, it follows from the definition (46) that being a parameter.The field equation ( 50) simplifies significantly when written in the null coordinates [33] u t − z and v t + z. ( Since there is no intersection between the region where the source of emission of the GWs is active and the past lightcone of π, there is a translational invariance along the latter null coordinate [33].This suggests the änstaz π = π(u) for the solutions of the field equation (50).By defining equation ( 50) reduces to [33] φ where prime denotes derivation w.r.t. the null coordinate u.In addition, the relations (52) have been used in the second equality.From this result it follows that φ is always decelerating along the null direction u.
The stability of a general solution to the equation ( 56), namely π, can be addressed by analysing the tensor structure of the kinetic part related to the quadratic action for small fluctuations around that solution.That is, by considering π = π + δπ and expanding the Lagrangian density (45) around the background solution π coming from equation (56).This leads to [33] (see also reference [73]) where the values of the kinetic matrix Z µν depend on the spacetime coordinates through the background solution π.These are [33] η cos 2 (ωu) , (58b) η cos 2 (ωu) , (58c) where φ ′′ is given by equation ( 56).(Recall that we are considering here only the case of subluminal propagation for π, i.e. c 2 s < 1. See reference [33] for a discussion when c 2 s ≥ 1.)It should be highlighted that for a general time-dependent kinetic term like (57) there is no clear definition of stability [73].However, in the limit where the time and length scales considered are much shorter than the rate of variation of π, it is perfectly acceptable to analyse the stability of the system as if the kinetic matrix Z µν were constant [73] (see also reference [33]).Within this local approximation to the stability of the theory, the absence of a gradient instability will be determined by the correct relative signs between the terms involving time and spatial derivatives in the Lagrangian density (57).It is convenient, however, to move to Fourier space and directly analyse the dispersion relation for a mode with four-momentum k µ = (ω, k i ) since the kinetic matrix Z µν is not in diagonal form.For the Lagrangian density (57), the dispersion relation reads [73] Moreover, multiplying this relation by Z 00 , it is possible to reexpressed it as [73] when Z 00 is not trivial.We recall that indices have been raised/lowered with the Minkowski metric.The absence of a gradient instability demands all frequencies, ω, to be real for any wave vector k i .Otherwise, exponentially growing modes will be present.It is straightforward to notice from the above expression that this stability is guaranteed if the spatial matrix Z 0i Z 0 j − Z 00 Z i j is positive defined [73].Since this matrix is diagonal for the Z µν introduced in equations ( 58), in our case of interest the positive-definedness condition reduces to for all spatial directions.In addition to the previous condition, the ghost-like instability is avoided if [73] Z 00 > 0. ( It is interesting to note that a theory with Z 00 < 0 can be also made stable provided it features superluminal excitations and that one can boost to a frame with Z 00 > 0 [74].However, we will not discuss this case here.Also note that the conditions (61) and ( 62) reduce to the standard ones when the matrix Z µν is in diagonal form [43].When η is always trivial (that is when considering the case of cubic Galileon only, i.e.G ∝ X), the ghost-free condition (62) leads to which follows from demanding the factor ahead of the trigonometric term in equation (58a) to be less than one.Moreover, the expression (Z 03 ) 2 − Z 00 Z 33 is always positive when η ≡ 0. Therefore, there is not a gradient instability in the z direction.
The background solution π, on the other hand, does not affect the entries Z 11 and Z 22 .For | β | > 1, these entries oscillate between positive and negative values in a way that condition (61) could never be fulfilled.Avoiding this oscillations in the (x, y)-plane demands For a non-relativistic speed of sound, that is c s ≪ 1, condition ( 64) is clearly more restrictive than the ghost-free condition (63).Therefore, by demanding | β | < 1 the stability conditions ( 61) and ( 62) are simultaneously fulfilled in the non-relativistic regime.
The complete stability region (including relativistic c s ) is shown in green in figure 2. As there can be seen, | β | < 1 is a necessary condition for the absence of ghost and gradient instabilities when η is trivial 4 .Owing to the definition of β in equation ( 53), this condition results in a tight constraint over m3 1 [33].Please note that there is a direct relation between m3 1 and α B , which essentially measures the amount of braiding in the theory.This is as follows from comparing equations ( 20) and (B.43).The observational bounds for α B coming from the stability condition of | β | < 1 are discussed in reference [33].To that aim, the authors there have considered a GWs-background as that generated by the event GW150914 [75] with a chirp mass of M c = 28M ⊙ and f = 30 Hz at a distance of 1 Mpc.They found that the condition of | β | < 1 directly translates to |α B | ≤ 10 −2 for that event [33], see definition of β in equation (53).Thus, they conclude that the braiding part should not have any sizeable effect in order to avoid the destabilising interaction with GWs when η is zero [33].Similar conclusions apply also to the case of c s = 1 [33].For the superluminal case, c s > 1, it has been argued that the same results should hold as well, but no explicit calculations were provided due to technical difficulties [33].
On the contrary, if η 0, then, the oscillations in the spatial entries Z 11 and Z 22 can be avoided even for | β | > 1 provided that η is sufficiently large and positive (see appendix A of reference [33]).Hence, this may provide a mechanism for having a stable theory with | β | > 1.If that is the case, then, the strong observational bounds for α B discussed in reference [33] shall be revised.The larger the value of η, however, the harder avoiding a ghost instability will be since this quantity enters linearly in the factor ahead of the trigonometric term in expression (58a).Therefore, one should carefully explore the region of viability for the optimal effect when η > 0. We address this discussion in the next sections.

Loophole hunting
Condition (64) places a strong restriction to the maximum amount of braiding a theory may have in order to avoid the destabilising interaction with GWs described in reference [33].As a result from that constraint, the braiding part should not have any sizeable effects when data from GWs events are taken into account [33].However, this restrictions holds only for the special case of η ≡ 0. In a more general situation, it would be possible to ease the constraint in equation (64) provided that η is positive enough to ensure that all spatial direction in the kinetic matrix Z µν are negative.This is obviously not the case for G ∝ X for which η is trivial [33], but may be possible for models beyond a linear function G.
Inspired by the fact that the ghost-free condition ( 63) is less restrictive than the constraint in equation ( 64) for nonrelativistic sound velocities when η is trivial, we explore here the possibility of a ghost-free and gradient-free theory with | β | > 1 when η 0. If that possibility comes to be true, then, the previously mentioned observational bound of |α B | ≤ 10 −2 [33] should be revised.(Recall the definition of β in equation (53).)This may indicate that KGB theory could have sizeable effects at cosmological scales after all.
In general, the ghost-free condition (62) implies which follows from demanding the factor ahead of the trigonometric term in equation (58a) to be less than one.The gradientfree condition ( 61) is automatically satisfied in the z direction since the combination is always positive for the values of interest for η.Whereas in the xy-plane the absence of a gradient instability implies Z 00 Z 22 < 0, (69) as follows from comparing the gradient-free condition (61) with the expressions for the Z µν found in equations (58).Contrary to the case of η ≡ 0, the oscillations in the spatial entries Z 11 and Z 22 can be avoided even for | β | > 1 provided that η is sufficiently large and positive.Taking into account their explicit expressions in equations (58b) and (58c), respectively, we found that for these terms are always negative.A model satisfying the above condition will feature positive entries Z 11 and Z 22 .Therefore, it will be free from gradient instabilities in the xy-plane provided that Z 00 is positive.That is to say, if the ghost-free condition (66) is also satisfied.In other words, condition (66) and the new constraint found in equation ( 70) ensure the total stability of the system.We should emphasise the fact that condition (70) does not depend on the specific properties of the GWs (local properties) but only on the value of c 2 s , which depends on the KGB functions.It is also interesting to note that if scalar perturbations propagate close to the speed of light, then, condition ( 70) is easily satisfied if η is just slightly positive.However, this regime would strongly constraint β to vanish, see equation (66).For a model with a non-negligible amount of braiding, this suggests that avoiding GWs-induced instabilities becomes more difficult if DE fluctuations propagate with a relativistic speed of sound.Conversely, the ghost-free condition (66) s is small.Nevertheless, η should diverge strongly than 1/c 2 s in order to avoid the gradient instability in this regime.This behaviour is summarised in figure 2. There we show different slices at constant η of the stability region where conditions ( 66) and ( 70) are simultaneously fulfilled.The green region corresponds to the case of η ≡ 0 (compare with figure 2 in reference [33]).As discussed before, the absolute value of β is always confined there to be less than one.Alternatively, this region becomes significantly bigger for small speed of sound if η is positive.In fact, stability at | β | ≈ 10 is possible for η of the order of twenty.Please note that this suppose an enlargement of the stability region by a factor of ten w.r.t. the green zone.Moreover, this region can be amplified even more at low c s for larger values of η.In the high c s regime, however, the stability zone narrows compared to the case of η trivial.
In view of these results, we postulate that the would-be KGB candidate that could ease (to some extent) the restrictions from the interaction with GWs described in reference [33] must meet the following necessary (but not sufficient) conditions.It should produce a non-relativistic speed of sound, c s , and a large value of η during the expansion history where interactions with GWs are expected.The definition of both quantities in terms of the KGB functions and its derivatives can be found in equations  66) and ( 70) are simultaneously satisfied.Within these regions, the interaction between DE fluctuations and GWs described in reference [33] does not induced a ghost or gradient instability.The green zone represent the stability at η = 0.This is the corresponding region for cubic Galileon, i.e.G ∝ X (compare with figure 2 in reference [33]).The darker shaded areas represent the stability at different non-zero values of η.
(18) and ( 40), respectively.However, if for the selected model DE perturbations propagate at relativistic speeds, then, smaller values of η are preferred.In any case, the contribution of the braiding part in that regime should be negligible.
It should be also mentioned that, within the spirit of the discussion in past sections, one may consider the possibility to relax the ghost-free condition (66) and be concerned only with ensuring a gradient-free spectrum.Of course, this will drastically reduce the restrictions on the theory.Although the possibility of a ghost-mode in cosmology sets an interesting discussion, here there is a technical difficulty in pursuing that direction.Given the structure of Z 00 in equation (58a), if the condition (66) is violated then the sign of Z 00 is not always negative (ghost mode excited) but fluctuates.Moreover, this fluctuation will jeopardise the relations ( 68) and ( 69) since the functional dependence of Z 00 , Z 11 and Z 22 on cos(ωu) is different.Thus, if the ghostfree condition is violated, then, gradient instabilities in the xyplane will also (periodically) appear within each fluctuation of the GWs.

No natural dark energy
A natural KGB candidate to test our hypothesis is that when η is a constant.From its expression in terms of G and its derivatives (40), this implies that being c G a coupling constant.That is a power-law prescription for the braiding function G. Consequently, the case of η constant provide not only a natural choice to test our proposal but a physically relevant model.For this model, the resulting speed of sound squared should be constrained as in order to fulfil condition (70).This results in a lower limit for c 2 s that should be satisfied (at least) during the cosmological evolution where interactions with GWs are expected.Please note that this is in agreement with the results presented in reference [76], where the possibility of a vanishing speed of sound was shown to be disfavoured.It should be also highlighted that, even though having a constant η is a very restrictive consideration, the results from exploring this case should also hold (at least asymptotically) for any other model for which η converges to a constant value.In other words, the physical conclusions from exploring a model like (71) should also apply asymptotically to any KGB theory for which the function G converges, at some point, to a power-law function.
In order to explore whether this mechanism can be conducted by a reasonable KGB model, we select a power-law function for the k-essence part as well.That is, we also consider where α is a constant and c K equals to +1 or −1.For the KGB model given by the functions ( 71) and ( 73), we have numerically explored the parameter-space spanned by {(η > 0) × α × c G }.As our main purpose, at this stage, is to examine the stability of the KGB theory, not to confront it with observational data, we have considered the initial conditions H 0 = 73.2km s −1 Mpc −1 , Ω ϕ0 ≈ 0.68, Ω m0 ≈ 0.32 and Ω r0 ≈ 10 −5 at face value for the numerical integrations (in line with local direct measurements as those reported in reference [64]).
Our numerical analysis has shown two main tendencies for these models.First, we found that the ghost-free and gradientfree conditions discussed in section 2 (where scalar perturbations decouple from GWs at linear order) are not fulfilled for all the parameter-space, but only in some regions.Second, for the models that satisfy the ghost-free and gradient-free at linear order, we have found that the condition ( 70) is always violated at some point during the (past) evolution of the system.We summarise these findings with the two proxy KGB candidates shown in figures 3 to 6.In the former, we consider η = 0.4, α = 0.6 and c G ranging form 9 • 10 −11 to 9.3 • 10 −11 .As can be appreciate in figure 3, this simple KGB model has a compelling behaviour on the background level.Depending on the value of the coupling constant c G , the fractional contribution of the scalar field to the total energy of the universe picks at matter-radiation equality, rendering this model as a potentially interesting Early Dark Energy (EDE) candidate (see a review on EDE in, for instance, reference [83]).Moreover, the resulting evolution of the Hubble rate at low redshift can be, in principle, made compatible with the current data available from direct measurements 5 .At the level of linear scalar perturbations, the ghost-free and gradient-free conditions are always satisfied, even though the model at hands shows a slight phantom behaviour at present time; see figure 4. Nevertheless, as can be seen in figure 5, the condition ( 70) is never fulfilled for the domain of the numerical integration (that is from z i = 10 6 to z f = −0.999).Since this bound is always violated, the signs of Z 11 and Z 22 will alternate between positive and negative values if | β | > 1.In that case, conditions ( 68) and ( 69) would not always be satisfied during the period of a single oscillation of the GWs.As a result, the interaction between DE fluctuations and GWs will induce (at least) a gradient instability for certain (periodic) values of the argument ω(t − z) if the amount of braiding is non-negligible.Hence, this model cannot provide a working example for the mechanism we discussed in section 3.2.
A second example is shown in figure 6.This corresponds to η = 2, α = 2 and c K = −1 in expressions ( 71) and (73).In this case, the coupling constant c G takes values from 4.5 • 10 −4 to 4.8 • 10 −4 .The background evolution of this proxy model is still compatible with low redshift measurements, but shows some tensions with data above z ≈ 1 for some values of c G ; see left panel in figure 6.The right panel shows the fulfilment of the condition (70).This condition is always satisfied in the asymptotic past of the model (within the limits of the integration domain).Unfortunately, it is consistently violated when the scalar field becomes dominant, if not even before that moment.As a result, interaction between DE fluctuations and GWs will induce (at least) a gradient instability if the amount of braiding is non-negligible.
The models represented in figures 3 to 6 illustrate the general behaviour we found for the functions ( 71) and (73).Even though there exist regions in the parameter-space for which D > 0 and c 2 s ≥ 0 always, the condition ( 70) is, in the best case scenario, fulfilled only for a finite time during the (past) evolution of the model.Moreover, this condition is systematically violated (at the latest) when the scalar field becomes dominant.Hence, a GWs-induced instability, like that described in reference [33] if the amount of braiding is non-negligible, seems to be inevitable for the models at hands.However, it is important to emphasize that the numerical screening of the (η, α)-plane we have performed is not an analytic proof of the non-viability of the KGB example given by the functions ( 71) and (73).Moreover, a different choice for the K function may also affect the results.Nevertheless, given the naturalness of power-law functions in physics and the lack of a positive conclusion in our numerical analysis, we suspect that finding a KGB model that could conduct the mechanism we discussed in section 3.2 for avoiding a GWs-induced instability [33] would be rather difficult, if not impossible.

What is left?
The braiding term G2ϕ in the scalar field action (1) has been shown to interact with GWs in such a way that it produces a ghost and/or a gradient instability in the scalar perturbations if the amount of braiding present at cosmological scales is nonnegligible [33].Motivated by the fact that the stability analysis in reference [33] was mainly performed for the special case of cubic Galileon, i.e.G ∝ X, in this paper we have proposed a theoretical mechanism to avoid these GW-induced instabilities for KGB models beyond cubic Galileon.This mechanism  71) and ( 73), respectively.In addition, the coupling constant c K equals minus unity, whereas c G ranges from 9 • 10 −11 to 9.3 • 10 −11 .Left panel represents the evolution of the partial energy densities for radiation, matter and the scalar field.Right panels show the evolution of the Hubble rate at low redshift compared to direct measurements collected from references [77][78][79][80][81] (see also Table II in reference [82]).71) and (73), respectively.In addition, the coupling constant c K equals minus unity, whereas c G ranges from 9 • 10 −11 to 9.3 • 10 −11 .Left panel represents the speed of sound squared, c 2 s , as defined in equation (18).Right figure shows the evolution of the ghost parameter D defined in equation (17).As can be seen, linear order perturbations (where scalar perturbations are decoupled form GWs) are always stable even though the scalar field displays a phantom behaviour at present epoch.relies on having a large and positive η (quantity that measure deviations w.r.t. a linear function G) for expanding the stability region to | β | ≫ 1 at low values of the speed of sound, as shown in figure 2. It should be mentioned, however, that the parameter β is not only proportional to the amount of braiding in the theory (through m3 1 or, equivalently, α B ) but to the strain amplitude of the GWs as well, see equation (53).If observational data of GWs events with higher and higher strain amplitudes are considered, then the resulting m3 1 (or α B ) will still be constrained to be negligible, even though our mechanism could extend the stability region to larger values of β.In that scenario, the conclusion in reference [33] about the non-viability of the braiding at cosmological scales will be only marginally affected by the discussion we have proposed here.Nevertheless, before finally concluding that the KGB theory should not have any sizeable effects, the fate of these instabilities, once originated, should be addressed.This lays beyond the scope of the present work.
We have also explored whether some reasonable KGB candidate could conducted the theoretical mechanism proposed here for evading GWs-induced instabilities.Unfortunately, using numerical simulations with a power-law prescription for the functions K and G we have not been able to find a working example for this effect.As we commented in section 3.3, the power-law choice for G provides a physically relevant case for our mechanism to work.Therefore, we consider the lack of a working example with power-law functions as a solid hint towards the non-viability of a cosmological braiding like the one present in the shift-symmetric KGB theories.
Conversely, the k-essence part in the KGB action (1) trivially evades the constraints coming from references [33][34][35].Moreover, this scalar field theory can feature ghost-free and gradient-free scalar perturbations if the NEC is satisfied for the DE sector.As we have discussed, however, the ghost-free and gradient-free conditions for k-essence are anti-correlated when the NEC is violated.Therefore, if we were to describe a DE component that violates the NEC (a possibility that is observationally viable [30]) with a single scalar field within the kessence theory, then we will inevitably face a ghost or a gradient  71) and (73).Taking c K = −1, the curves represent the ratio (1 − c 2 s )/2c 2 s from equation ( 70) for different values of the coupling constant c G .The horizontal line depicts the upper-bound for this ratio found in expression (70).Clearly the model never satisfies this bound for the domain of the numerical integration and, therefore, the interaction with GWs always introduce (at least) a gradient instability in the scalar sector.instability.Whether the instability will be ghost-or gradientlike depends, ultimately, on the function K.It is also important to note that k-essence cannot produce a phantom crossing [28,29].Hence, the scalar field will always be phantom or nonphantom, with the stability issues this may entail in the former case.This points towards the impossibility of solving simultaneously the H 0 and S 8 tensions with late-time modifications to ΛCDM à la k-essence should the S 8 tension be still present in future observations, as this resolution would demand a phantom crossing [31,32].See also previous works on the H 0 tension in the k-essence theory in, for instance, references [61,62].With this potential limitation in mind, an interesting extension of the k-essence scenario is that when more than one scalar field is present [84,85].In that case, the effective DE fluid can cross the phantom divide [86] even though each scalar field, namely ϕ i , is restricted to the phantom or non-phantom regime in the same way as described in section 2.1.
It should be noted that the discussion presented in this work will be only partially affected if we allow the shift-symmetry to be (weakly) broken at some regimes.This would introduce an explicit dependence on the scalar field ϕ in the functions K and G.The destabilising interaction between GWs and the braiding term will still occur even for G = G(ϕ, X) [33].However, it would be interesting to explore whether the mechanism proposed in section 3.2 will have a better performance in that case.If not, then, the braiding term should not have any able effect and the surviving theory will be effectively that of kessence.For this non-shift-symmetric k-essence theory, no additional terms will appear in the formulas for the scalar field energy density and pressure, or in the equations for c 2 s and D; see the corresponding expressions in, for instance, reference [16].Thus, our discussion about the linear stability of the theory shall remain the same as for the shift-symmetric k-essence case.In addition to the previous considerations, the k-essence scalar field cannot produce a phantom crossing even if K = K(ϕ, X); see reference [28,29].Therefore, the conclusions presented in section 2.1 shall remain unchanged.
Along this line of thought, another possible extension of the k-essence theory is that when a non-minimal coupling to gravity is allowed for the scalar field.Assuming that the external matter field are minimally coupled to the metric g µν , the gravitational part of the total action would read where the non-minimal coupling function f can only depend on the scalar field, but not on its kinetic term X.This is due to strong restrictions on the speed of propagation of GWs; see reference [17].Please note that this scalar field theory would also trivially evade the constraints from [33].Nevertheless, the nonminimal coupling may introduce a fifth-force at local scales.Therefore, some mechanism should be invoked to screen this force on astrophysical scales; see, for instance, references [87][88][89] and references therein.This limitation may also be present in extensions of k-essence to modified theories of gravity like those discussed in, e.g., references [90][91][92][93][94][95].Furthermore, within the spirit of analysing fundamental symmetries in order to address DE, we should also briefly mention the possibility of considering a scalar field with an action that is invariant just under transverse diffeomorphisms [82,[96][97][98].
In this case, one can describe phenomenology of interest for the dark sector even with a shift-symmetric field with a canonical kinetic term [82,98].On the other hand, one can consider more general theories invariant only under transverse diffeomorphisms [99,100], being the most prominent example that of Unimodular Gravity [101] (see also [102][103][104][105][106][107]).Finally, as it is well-known, more general alternative theories of gravity goes beyond the consideration of a scalar field to describe the dark sector [108].
Finally, we shall recall again that we have discussed the viability of the KGB theory solely from the point of view of the stability of scalar perturbations.In doing so, we have focus on theoretical considerations only.Should any of the models studied had fulfilled our stability criteria, then the analysis of background expansion by means of the evolution of the partial energy densities and the Hubble rate (like shown in figures 3 and 6) would not be sufficient to claim for the viability of the DE-candidate.In such case, a thorough analysis on the observational impact using more advanced techniques and specialised software like CLASS should be carried on; see, for instance, references [109,110].We leave for future research the stability of more involving KGB models and their compatibility with observational data.II in reference [82]).Right panel represents the ratio (1 − c 2 s )/2c 2 s from equation ( 70).The horizontal line there depicts the upper-bound this ratio should satisfy to avoid a GWs-induced gradient instability; see inequality (70).

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A. Stability of linear scalar perturbation
In this appendix we discuss the possibility of satisfying simultaneously the ghost and gradient-free conditions for linear scalar perturbations given by the inequalities (36) and (39), respectively.In order to make our notation more compact, it will be useful to introduce the two following reference scales where Ω ⋆ is well-defined only for Ω B J 2. The case of Ω B J = 2 can be directly addressed using the expressions (36) and (39).We consider here the parameter η to be positive since that is the case proposed in section 3.2 for avoiding a GWs-induced gradient instability.Also note that the quantity Ω r + 3 − 3Ω ϕ is always positive 6 .Hence, Ω ⋆ is positive if Ω B J ∈ (0, 2), whereas it is negative for Ω B J < 0 and Ω B J > 2. Combining the ghost and the gradient-free conditions lead to the following scenarios in which both conditions are satisfied for η non-negative: Scenario 1. (Ω B J < 0, Ω J > max{(1 + w ϕ )Ω ϕ , Ω ⋆ , Ω}): Both the numerator and the denominator in equation (36) are positive.The reference scale Ω ⋆ is negative.Phantom behaviour is allowed.
The numerator and the denominator in equation (36) are negative.The reference scale Ω ⋆ is positive and, therefore, so it is Ω J .Hence, from the upper bound for Ω J it follows that the scalar field cannot be phantom.

Scenario 4. (Ω
The numerator and the denominator in equation ( 36) are positive.Note that Ω > (1 + w ϕ )Ω ϕ if the scalar field is phantom-like and, therefore, Ω J should be greater that Ω.In the non-phantom regime, the lower bound for Ω J is given by (1 The numerator and the denominator in equation ( 36) are positive.The reference scale Ω ⋆ is negative.Phantom behaviour is allowed.
Even though the physical intuition behind each of these scenarios is not so transparent as for the marginal models studied in sections 2.1 and 2.2, their very existence is a solid proof that it is possible to have a ghost-free and gradient-free scalar perturbations at linear level (recall we consider here only the case of η non-negative) in the shift-symmetric KGB theory even when w ϕ < −1.

Appendix B. Quick guide to the EFT of dark energy
The framework of EFT was first applied to DE in reference [36] (see also [111,112] for applications to inflationary models) and further developed in references [113][114][115][116] (see also reference [45] for a review).This approach considers the most general form for the gravitational action (including cosmological perturbations up-to an arbitrary order) built-up only on symmetry arguments.Originally, this framework was based on two assumptions [113]: i) the validity of the weak equivalence principle and the existence of a Jordan metric g µν universally coupled to matter fields; ii) the existence of a unitary gauge compatible with the residual symmetries of unbroken spatial diffeomorphisms.
The last of these assumptions advocates for the presence of a scalar field ϕ in the DE sector.This is so because the presence of a homogeneous scalar field ϕ(t) in a FLRW background defines a preferred time slicing of space-time.These are the hypersurfaces with ϕ constant and, therefore, δϕ ≡ 0 on them.So, for this choice (unitary gauge) only metric degrees of freedom are explicitly displayed in the action.The scalar field perturbation can be reintroduced explicitly in the theory via the Stüeckelberg trick.That is, by performing infinitesimal time diffeomorphism t → t + π(t, x) being π the scalar field fluctuations.However, in this approach π does not represent the original scalar field ϕ in the DE sector, but the perturbations encoding the scalar degree of freedom in the theory.
In order to construct the most general expression for the action satisfying the previous ansätze, note that hypersurfaces with constant ϕ can be defined as those orthogonal to the unit four-vector provided that ∂ µ ϕ is time-like, i.e.X > 0, where ϵ sgn( φ) = ±1 is to ensure the four-vector is future-oriented.(Here we are not going to consider oscillating solutions where the background field velocity, φ, crossed zero since this case may be problematic in standard perturbation theory [41].)For a homogeneous scalar field, that is ϕ = ϕ(t), this four-vector reduces to This slicing induces the spatial metric where n µ n µ = −1 and n µ h µν = 0. Note that invariance under time-translations is spontaneously broken (in the sense discussed in reference [116]) since the scalar field signals out a preferred time.The terms allowed in the EFT action are, therefore, those invariant under the residual symmetries of the unbroken spatial diffeomorphisms, such as the contravariant time-time component of the Jordan metric g 00 .The extrinsic curvature of constant time hypersurfaces is also allowed to appear.This is defined as [117] and represents the spatial projection of the covariant derivative of n ν .(Note that the extrinsic curvature is sometimes defined with a different sign convention, see, for instance, references [118,119].)In addition, its trace reads It will be useful for latter calculations to note that up-to boundary terms and for any function l on time.The Ricci scalar R, any curvature invariants, and contractions of tensors with g µν , n µ and the covariant derivative ∇ µ are also invariant under spatial diffeomorphisms.Therefore, the most general EFT action can be expressed as [111] where time is also allowed to appear explicitly.For addressing the cosmological perturbation, however, it is more convenient to re-write the above action into the part contributing to the background level and that enclosing the perturbations around a flat FLRW metric (at any order).This is [113] (sea also [111,115]) where f , Λ and c are functions of the cosmic time t.Moreover, the latter two background quantities read being p DE and ρ DE the pressure and energy densities of the DE fluid.In addition, M 2 * is the Planck mass.Hereon we will take M 2 * = 1 since we have adopted the geometric unit system in this paper.Conversely, S (2)  DE contains all terms that start at quadratic order in perturbations and, therefore, they do not affect the background.(Note, however, that the background quantities do affect all perturbation levels.)This part can be expressed as follows [113]  i and µ i are mass parameters whereas λ i and γ i are dimensionless quantities.Generally, these parameters are allowed to depend on the time coordinate t.In addition, δg 00 g 00 + 1 is the perturbation of the metric, δK µ ν is the perturbation of the extrinsic curvature, and δK the perturbation of its trace.Likewise, δR µν is the perturbation of the Ricci tensor and δR the perturbation of its trace, whereas C µνρσ represents the Weyl tensor.The action (B.8) describes virtually all DE models encompassing a single scalar degree of freedom.
For the case of the KGB theories given by action (1), the background quantities in (B.8) read [113] f = 1, (B.12) with ρ ϕ and p ϕ given by equations ( 9) and (10), respectively.Whereas the action (B.11) for the perturbations reduces to [113] That is, only the mass parameters M 4 n and m3 n−1 (for n ≥ 2) are present for the KGB theory (see also, for instance, reference [45]).Recall that these parameter are, in general, functions on the time coordinate t.Moreover, note that M 4  2 is connected with the kineticity term α K , whereas m3 1 is equal to the braiding term α B up-to sign conventions; see the discussion in Appendix B.3 (see also table 2 in reference [41]).It is also important to highlight that the operator introducing the kinetic mixing with gravity at leading order is m3 1 .This can be understood in an intuitive way by noting that the perturbation δg 00 introduces a term proportional to π at linear order after the Stückelberg procedure, see equation (42).The extrinsic curvature, on the other hand, contains time derivative of the spatial metric, h i j , when time diffeomorphism is fully restored [117] (see also, for instance, references [113,120,121]).Schematically, the product δg 00 δK, therefore, introduces a term proportional to πḣ (at second order) into the Lagrangian.In the Newtonian gauge (41), this contribution is proportional to π Ψ; see reference [113].For the same reason, the operators m3 n−1 (for n ≥ 2) will also be responsible for introducing the kinetic mixing between scalar perturbations and the GWs at higher orders.
In the following part of the appendix, we compute the expression of these mass parameters in terms of the KGB functions K and G, and their derivatives.These calculations will be done without invoking a shift-symmetry and, therefore, the resulting expression are valid for the most general KGB scenario.

Appendix B.1. Mass parameters for k-essence
As a warm-up, we consider first the k-essence theory.In the most general scenario, the k-essence function depends on both the scalar field and its kinetic term; i.e.K = K(ϕ, X). (We remind the reader that we will not invoke a shift-symmetry here; thus, the results obtained in this appendix are valid for the most general case.)When perturbing this function, one may naively consider that perturbations around the background values of both ϕ and X should be taken into account.However, in the unitary gauge (defined by δϕ ≡ 0) only the perturbations of the kinetic term remains.Let the perturbed kinetic term be defined as where X is the background value and δX the perturbation.In the unitary gauge, the perturbed kinetic term X is related to the perturbed metric through X = − 1 2 φ2 g 00 = −Xg 00 , (B.17 where we have assumed that ϕ and its kinetic term are functions of time only at the FLRW background.Note that the scalar field does not appear explicitly in the unitary gauge: it is "eaten" in the metric degrees of freedom [45,113].However, in order to obtain the expressions for the mass parameters in the EFT language in terms of K(ϕ, X) and its derivatives, it will be useful to maintain the usual covariant notation for the function K and its arguments instead of K(t, g 00 ) when there is no risk for confusion.Please note that this is an abuse of notation since neither ϕ nor X can appear explicitly in this gauge.They are simply functions on the time coordinate.In this notation, the k-essence function can be expanded around its background value as K(ϕ, X) = K(ϕ, X) + where subscript "δX = 0" denotes evaluation of the corresponding quantity on the background level.For the sake of the compactness of the notation, we will avoid writing explicitly the arguments of the function K and its derivatives appearing in the r.h.s. of the expansion, where it should be kept in mind that these expressions are always evaluated on the background.We will also omit the subscript "δX = 0" in the following and express the coefficient in the expansion simply as ∂ n K/∂X n where there is no risk for confusion.Accordingly, the above expansion can be rewritten as where we have substituted δX = X − X in the first order in the expansion, and used equations (B.17 for n ≥ 2. Therefore, the k-essence function K contributes only to the mass parameters M 4 n at perturbation level.So, there is no kinetic mixing between scalar and tensor perturbations.

Appendix B.2. Mass parameters for G2ϕ
From the definition of the d'Alembertian operator and the slicing vector n µ (see equation (B.1)), it is straightforward to check that To compute the EFT parameters for the braiding part of the KGB action (1), it will be useful to consider the following expansion in polynomials of δX [113] G(ϕ, X) = ϵ where X = X + δX is the perturbed kinetic term and , (B.27) 7 Please, note that the authors in reference [113] used a different convention for the definition of the kinetic term.This is X g µν ∇ µ ϕ ∇ ν ϕ [113], in contrast with our definition for X − 1 2 g µν ∇ µ ϕ ∇ ν ϕ, where the convention used for the signature of the metric tensor is in both cases the same.Also note that it is always possible to redefine the scalar field in a way that ϕ(t) = √ 2t and, therefore, the kinetic term would simply read X = 1 when evaluated on the background [113].However, we will not consider such redefinition here in order to explicitly keep the kinetic term, X, in the final expressions for the EFT parameters.
are the coefficient in the expansion evaluated on the background.Recall that we are not imposing a shift-symmetry in this section.Consequently, the above coefficients depend on the background values of both the scalar field and its kinetic term.That is l m = l m (ϕ, X).Let us emphasise again that this is an abuse of notation, since neither ϕ nor its kinetic term can explicitly appear in the unitary gauge.As mentioned in the previous section, background quantities as ϕ and X are functions of time only, whereas scalar perturbations are encoded in time-time component of the metric.Hence, the perturbed braiding function, G, and the expansion coefficients, l m , should read G(t, g 00 ) and l m (t) in this gauge choice, respectively.Nevertheless, this abuse of notation will be convenient, where there is no risk of confusion, for obtaining the expressions of the EFT parameters in terms of G(ϕ, X) and its derivatives.

Figure 1 :
Figure1: Stability of scalar perturbations in the (Ω ϕ , w eff )-plane for K ≡ 0. The lower oblique-lined region corresponds to the presence of a the ghost instability.Scalar perturbations are both ghost-and gradient-unstable in the central vertical-lined region.Stability is ensured in the w eff > 1/3 zone, since w ϕ is clearly positive there in virtue of the relation w eff = Ω r /3 + w ϕ Ω ϕ .In the dotted region, which corresponds the horizontal band 0 ≤ w eff ≤ 1/3, perturbations may be healthy if and only if w ϕ is non-negative, otherwise they will be both ghost-and gradient-unstable.Nevertheless, the sign of w ϕ in that region of the (Ω ϕ , w eff )-plane cannot be inferred without solving the whole evolution of the model and, therefore, the stability conditions there cannot be further addressed without specifying the function G.

Figure 2 :
Figure 2: Slices at constant η of the stability region where conditions (66) and (70) are simultaneously satisfied.Within these regions, the interaction between DE fluctuations and GWs described in reference[33] does not induced a ghost or gradient instability.The green zone represent the stability at η = 0.This is the corresponding region for cubic Galileon, i.e.G ∝ X (compare with figure2in reference[33]).The darker shaded areas represent the stability at different non-zero values of η.

Figure 3 :
Figure 3: Numerical evolution for η = 0.4 and α = 0.6 in expressions (71) and(73), respectively.In addition, the coupling constant c K equals minus unity, whereas c G ranges from 9 • 10 −11 to 9.3 • 10 −11 .Left panel represents the evolution of the partial energy densities for radiation, matter and the scalar field.Right panels show the evolution of the Hubble rate at low redshift compared to direct measurements collected from references[77][78][79][80][81] (see also Table II in reference[82]).

Figure 4 :
Figure 4: Numerical evolution for η = 0.4 and α = 0.6 in expressions (71) and(73), respectively.In addition, the coupling constant c K equals minus unity, whereas c G ranges from 9 • 10 −11 to 9.3 • 10 −11 .Left panel represents the speed of sound squared, c 2 s , as defined in equation(18).Right figure shows the evolution of the ghost parameter D defined in equation(17).As can be seen, linear order perturbations (where scalar perturbations are decoupled form GWs) are always stable even though the scalar field displays a phantom behaviour at present epoch.

Figure 5 :
Figure 5: Numerical evolution for η = 0.4 and α = 0.6 in expressions (71) and(73).Taking c K = −1, the curves represent the ratio (1 − c 2 s )/2c 2 s from equation (70) for different values of the coupling constant c G .The horizontal line depicts the upper-bound for this ratio found in expression(70).Clearly the model never satisfies this bound for the domain of the numerical integration and, therefore, the interaction with GWs always introduce (at least) a gradient instability in the scalar sector.