Gravity in the Era of Equality: Towards solutions to the Hubble problem without fine-tuned initial conditions

Discrepant measurements of the Universe's expansion rate ($H_0$) may signal physics beyond the standard cosmological model. Here I describe two early modified gravity mechanisms that reconcile the value of $H_0$ by increasing the expansion rate in the era of matter-radiation equality. These mechanisms, based on viable Horndeski theories, require significantly less fine-tuned initial conditions than early dark energy with oscillating scalar fields. In Imperfect Dark Energy at Equality (IDEE), the initial energy density dilutes slower than radiation but faster than matter, naturally peaking around the era of equality. The minimal IDEE model, a cubic Galileon, is too constrained by the cosmic microwave background (Planck) and baryon acoustic oscillations (BAO) to relieve the $H_0$ tension. In Enhanced Early Gravity (EEG), the scalar field value modulates the cosmological strength of gravity. The minimal EEG model, an exponentially coupled cubic Galileon, gives a Planck+BAO value $H_0=68.7 \pm 1.5$ (68\% c.l.), reducing the tension with SH0ES from $4.4\sigma$ to $2.6\sigma$. Additionally, Galileon contributions to cosmic acceleration may reconcile $H_0$ via Late-Universe Phantom Expansion (LUPE). Combining LUPE, EEG and $\Lambda$ reduces the tension between Planck, BAO and SH0ES to $2.5\sigma$. I will also describe additional tests of coupled Galileons based on local gravity tests, primordial element abundances and gravitational waves. While further model building is required to fully resolve the $H_0$ problem and satisfy all available observations, these examples show the wealth of possibilities to solve cosmological tensions beyond Einstein's General Relativity.

Observational cosmology has established a simple and successful standard model of the Universe. ΛCDM is named after the dominant components: the cosmological constant (Λ) accelerates the expansion at late times and cold dark matter (CDM) drives the formation of large-scale structure (LSS). In addition, the model includes other matter species known from Earthly experiments (atoms, photons, neutrinos) and assumes the va-lidity of Einstein's general relativity (GR). This remarkably simple model successfully describes most cosmological observations in terms of a handful of parameters [1]. But despite ΛCDM's success, several datasets interpreted within the standard model are in conflict [2].
The most significant tension involves the Universe's expansion rate. Late-universe measurements of H 0 clash with observations of early-universe processes interpreted within ΛCDM. Late probes include distance ladder [3] and lensing time delays [4,5]. They are direct and largely independent of the cosmological model. Probes based on early-Universe processes (or early probes) rely on Planck's cosmic microwave background (CMB) plus baryon acoustic oscillations (BAO) data [1]. Early probes are indirect and rely on the predictions of the ΛCDM model. Unless the Hubble problem is due to unknown systematics, its significance demands physics beyond the simple ΛCDM model [6,7].
New-physics solutions to the Hubble problem reflect the conflict between the early and late universe. Lateuniverse solutions rely on new astrophysical effects [8] or dark energy (DE) beyond Λ. Adjusting H 0 for fixed CMB angular scale requires that the density of DE grows in time instead of remaining constant, i.e. a phantom equation of state where E, P are the energy and pressure density of DE, respectively. While disfavoured by combining BAO and type Ia supernovae (SNe) [9][10][11], other analyses favour late-time solutions to the Hubble problem [4,[12][13][14]. Galileon gravity [15] once provided a late-universe solution to the Hubble problem. Simple models with Λ = 0 were compatible with Planck and BAO [16] and unambiguously predicted a value of H 0 in agreement with distance ladder, well before the Hubble problem was troubling. Latter investigations showed that the only Galileon compatible with CMB×LSS cross-correlation modify the speed of gravitational waves (GWs) [17]. The observation of coincident gravitational and electromagnetic radiation from the neutron star merger GW170817 [18] swiftly ruled out Galileons as a solution to the Hubble problem, along with many other theories of gravity [19][20][21][22].
Early-universe solutions invoke new physics before recombination to "recalibrate" the comoving acoustic scale which depends on the ratio of the sound speed and the expansion rate up to the redsfhift of baryon drag. These solutions work because BAO measure the dimensionless quantity H 0 r s . A larger value of H 0 requires decreasing r s by increasing H. Consistency between BAO and SNe (aka inverse distance ladder) introduces a relation between H 0 and r s which is largely insensitive to lateuniverse physics [9,[23][24][25]. Combined inverse and direct distance ladder prefer a shorter acoustic scale than Planck+BAO in ΛCDM (figure 1). This hints at an earlyuniverse solution to the Hubble problem. Early-solutions rely on new sources of energy contributing to the expansion rate before recombination, cf. Eq. (2). Possible scenarios include additional radiation [14,26], neutrinos with enhanced interactions [27], variation of fundamental constants [28] and non-standard dark matter [29][30][31]. Another idea is based on early dark energy, an analog to time-dependent DE but active in the early universe. Early DE can be studied via timedependent parameterizations of the energy density in the Friedmann equation [10,32,33] (see [34,35] for earlier works) and/or the effective gravitational constants in the evolution of perturbations [36].
Dynamical early DE models introduce a canonical scalar field to solve the Hubble problem [37][38][39]. A potential V (φ) with a minimum is required to combine the phenomenology of thawing quintessence [40] and damped oscillations [41]: the scalar field is initially subdominant and frozen by Hubble friction. It thaws as the energy density of matter becomes comparable to the potential. Then it begins rolling down the potential and oscillating around the minimum, losing energy in the process until it becomes subdominant again. Data requires that the scalar starts evolving around the era of equality, setting the initial condition φ i so V (φ i ) ∼ eV 4 . The scalar's energy density is constant before equality, with V (φ i )/ρ r (z BBN ) ∝ (T eq /T BBN ) 4 ∼ 10 −24 when compared around Big-bang nucleosynthesis (BBN). Canonical scalar fields require very fine-tuned initial conditions to solve the Hubble problem.
Studies of dynamical early DE models have been restricted to canonical scalars with different potentials and simple extensions [42,43]. This covers but a narrow sliver in the space of known gravitational theories [44]. It is plausible that novel signatures and interesting features (e.g. reduced fine-tuning) can be found among extensions of early DE. The goal of this work is to explore novel solutions to the Hubble problem in viable theories beyond GR, focusing on novel early Modified Gravity mechanisms and their phenomenology.

A. Summary and guide for the busy reader
This work considers three potential solutions to the Hubble problem in scalar-tensor theories of gravity: 1. Imperfect Dark Energy at Equality (IDEE): the scalar kinetic energy dilutes faster than matter but more slowly than radiation, naturally peaking in the era of equality (III B). Minimal IDEE can not reconcile Planck+BAO and distance ladder (IV B). In IDEE-only models (dashed) the stringent constraints limit the impact on rs. Coupled EEG models (solid) relax the bounds considerably, extending the degeneracy across the BAO+SNe direction. Uncoupled/coupled LUPE models with Λ = 0 (red/orange) predict a high central value of H0 compared to the canonical Λ = 0 cases (purple, dark green), but have a worse fit and' are ruled out by other observations. LUPE models with Λ = 0 (magenta) provide an intermediate case. (Figure adapted from [25]).

Enhanced Early Gravity
Ricci coupling and can increase the expansion rate at early times (III C). Planck+BAO allow EEG to accommodate higher values of H 0 , closer to the distance ladder measurement (IV C). Local tests of gravity strongly constrain EEG (V B).
3. Late-Universe Phantom Expansion (LUPE): the scalar energy density increases with time at low redshift, w φ < −1 (III D). LUPE models with Λ = 0 are ruled out, but coupled LUPE with Λ = 0 can ease the H 0 tension (IV D).
The constraints and evolution in each scenario are summarized in figures 1 and 2, respectively. While IDEE, EEG and LUPE are general mechanisms, the results refer to a coupled cubic Galileon scalar-tensor theory of gravity (II A). In this theory IDEE relies on the initial field velocityφ i and EEG on the initial field value φ i modulating the effective Planck mass (i.e. gravitational constant) via a coupling to curvature. LUPE requires a negative sign of the quadratic kinetic term (accelerating), causing the scalar energy density to grow in time. The three mechanisms can operate together or independently. Readers interested in either of the above mechanisms are directed to visit the sections cited above, in whatever order they consider appropriate. The main findings are summarized in the conclusions VI, along with ideas for further model building and observational tests. Section II introduces the class of viable Galileon theories. Their cosmological dynamics are presented in Section III. Section IV presents the cosmological constraints (Planck, BAO, distance ladder) and Section V discusses the challenges faced by coupled models, including BBN, local gravity tests and GWs. The appendices contain additional discussions.

II. GALILEON GRAVITY AFTER GW170817
This section presents the gravity theories studied as potential solutions to the Hubble problem. I will describe Galileon gravity theories and their status, focusing on constraints from cosmology and GW observations. In section II A I will narrow down to theories compatible with the GW speed and present (exponentially) coupled cubic Galileons, the class of models used here to investigate IDEE, EEG and LUPE.
Most Galileon gravity theories are specific realizations of the Horndeski class [45][46][47], the most general action for a tensor and a scalar field, generally covariant, Lorentz invariant and leading to second order equations of motion in 4 space-time dimensions: where L m = L m [g µν , ψ M ] is the matter Lagrangian density, minimally coupled to the Jordan frame metric g µν and the gravitational interaction is given by: The four Lagrangians L i encode the dynamics the scalar field φ of the Jordan-frame metric g µν . They contain four arbitrary functions G i (φ, X) of the scalar field and its canonical kinetic term, 2X ≡ −∂ µ φ∂ µ φ. Subscripts φ, X to denote partial derivatives, e.g. G iX = ∂Gi ∂X . I will follow the conventions of the hi class code [48,49], including natural units (c = h = 1) and mostly-plus signature of the metric (−, +, +, +).
The uncoupled covariant Galileon [15,50,51] is is the most general Horndeski completion of a theory realizing the Galilean symmetry φ → φ + c + b µ x µ in flat spacetime. The theory is defined by the following Horndeski functions where M P = 1/ √ 8πG is the reduced Planck mass and M 3 = (H 2 0 M P ) 1/3 . While the linear potential ∝ c 1 φ is compatible with Galilean symmetry, it does not lead to interesting phenomenology and it is common to set it to zero (see Ref. [52] for an analysis including c 1 ).
Uncoupled covariant Galileons (9) have interesting cosmological solutions, including an unambiguous prediction of H 0 compatible with distance ladder. If the quadratic kinetic term has negative sign (c 2 < 0), the theory predicts LUPE, accelerating solutions without the need of Λ. These Λ = 0 models require a sizeable neutrino mass m ν ∼ 0.6 eV to fit CMB+BAO observations [16], a value well within the range of laboratory experiments [53] but that excluded by cosmological data if assuming ΛCDM. The high value of the neutrino mass is also necessary to solve the H 0 problem: the m ν ≈ 0 models do not only yield a worse fit, but also predict a value of H 0 above the distance ladder measurement [54].
The evolution of the metric potentials constrains the parameter space of Λ = 0 uncoupled Galileons (9). The minimal theory (c 4 , c 5 = 0) always predicts growing potentials at low redshift, instead of decaying as in ΛCDM. The growth of the potentials leads to an anti-correlation between CMB temperature and the low redshift galaxies (CMB×LSS) via the Integrated Sachs-Wolfe (ISW) effect, in stark disagreement with current measurements [55]. General Galileons (c 4 , c 5 = 0) can accommodate decaying potentials in some regions of the parameter space [17].
Galileon theories beyond-Horndeski [59][60][61] can be made compatible with the GW speed but are ruled out by other GW observations [62,63] and cosmology [64]. Beyond-Horndeski theories in the Gleyzes-Langlois-Piazza-Vernizzi (GLPV) class can be constructed in which c g = 1 on any space-times [19,20,65]. GLPV Galileons with c g = 1 have identical cosmological expansion than their Horndeski analogs, potentially providing a late-time solution to the Hubble problem compatible with the GW speed. However, GLPV theories predict a very rapid decay of GWs into scalar field excitations [62,63], and the deviations from Horndeski need to be very suppressed for any GW signal to be detected. The remaining beyond-Horndeski term compatible with GW speed and decay does not have the Galileon form [59,Eq. 40] GW speed and decay bounds allow Horndeski theories with general G 2 , G 3 , but restrict L 4 , L 5 to G 4 (φ), G 5 = 0. Galileon theories are equipped with the Vainshtein screening mechanism [66] suppressing small-scale deviations from GR, including effects in the emission of GWs [67][68][69] (although see Ref. [70] for a possible counterexample). Still, these theories are still subject to GW constraints (instabilities induced by GWs [71] and standard sirens) which I will discuss in section V C.
A. Coupled Cubic Galileon I will explore the coupled Cubic Galileon [72][73][74][75][76][77][78][79], a variant of the Galileons described above, restricted to be compatible with GW observations but extended through a non-minimal coupling between the scalar field and and the Ricci scalar. The coupling is introduced via a φdependence of G 4 , the coefficient of the Ricci scalar in the Horndeski action or equivalently, (4)(5)(6). The main effect of the coupling is to modify the strength of gravity, now depending on the value of the field. 1 Coupled Cubic Galileon theories admit a binary classification into Canonical or Accelerating models. An arbitrary redefinition of the scalar field by a constant factor in the Lagrangian (10) fixes one of the coefficients without any loss of generality [54]. For a real-valued α, the above transformation always preserves the sign of c 2 , the quadratic kinetic term. Canonical and accelerating models models correspond to a positive and negative sign of c 2 , respectively. The differences between both will be explored in sections III D and IV D. The literature often refers to self-accelerating models, in which the universe's acceleration is supported to the conformal coupling. This is defined as the acceleration condition being satisfied in the Jordan frame (used here), but not in the Einstein Frame (see footnote 1). I will not consider self-acceleration further.
For further simplicity, I will consider only an exponential form of the coupling 1 The coupled theory (10) is a minimal extension of the uncoupled Galileons (9). One may also introduce the coupling directly into the matter action Lm(gµν , ψ M ) → Lm(C(φ)gµν , ψ M ) via the the so-called Einstein frame metricC(φ)gµν . This theory is not equivalent to Eq. (10) as G 2 is corrected by some terms depending onC(φ), see [80] for general expressions and Ref. [78] for an explicit example.
The exponential form is particularly simple to study.
All couplings with C φ = 0 break the shift-symmetry φ → φ + C, but the exponential coupling introduces only a constant term in the scalar field equation. Thus, there is no dependence on φ in the scalar field equations. The dependence on the scalar field value is thus limited to the gravitational sector, as φ modulates the strength of gravity. Compared to other choices of the coupling function, the exponential form leads to convenient simplifications in the analysis of the cosmology described in the next section.

III. COSMOLOGICAL DYNAMICS
In this section I will discuss the cosmological dynamics of coupled cubic Galileons (10), specializing to the exponential form of the coupling (12). Section III A introduces the dynamical equations and important concepts related to the theory. The following subsections detail how solutions of the coupled cubic Galileon lead to IDEE (III B), EEG (III C) and LUPE (III D). The early time dynamics are further discussed in Appendix C.

A. General Considerations
Let us start by presenting the general equations for the background metric and scalar field for coupled cubic Galileons. I will then review the classification of Galileons into canonical and accelerating and some properties of the exponential coupling.

Equations & Definitions
The expansion history is governed by the modified Friedmann equation where ρ m is the total matter density in CLASS units [Mpc −2 ] [81]. The effective Planck mass modulates the strength of gravity on the cosmological background and the kinetic energy densitŷ represents the remaining contributions of the scalar field to the expansion rate. Note that all the terms inÊ are proportional toφ, while M 2 * depends only on φ. The Galileon energy fraction today is then where the kinetic contribution (15) readŝ and the dimensionless field velocity [16] ξ ≡φ provides a convenient variable. The scalar field equation can be written in a current conservation formJ Here the shift-charge density (or shift-charge) is the time-component of a Noether current J µ associated to shift symmetry φ → φ + C. The kinetic term i.e. the coefficient ofφ in Eq. (19) determines the stability of the theory. It needs to be positive for the stability of both the background and linear perturbations. Finally, the source term is given by and is proportional both to the coupling strength β = C ,φ /C and the Ricci scalar evaluated on the cosmological background.

Canonical vs Accelerating Galileon
Let us now examine the kinetic structure and solutions of cubic Galileons, starting with the uncoupled case. Solutions to Eq. (19) with P φ = 0 correspond to the shift-charge density diluting with the Universe's volume. The scalar field is thus drawn towards J ∝ c 2 ξ − 6c 3 ξ 2 → 0, corresponding to two solutions The sign of the kinetic term c 2 determines which solution is stable via the no-ghost condition D > 0, Eq. (21). The above solutions reveal a binary classification of Galileons  (20) and the field derivative Eq. (18) is shown for canonical (c2 > 0, thick) and accelerating (c2 < 0, thin) models. Absence of ghosts requires a positive slope for the curve (21), with the minimum of J corresponding to the transition to instability. Stable accelerating/canonical models tend to ξ = 0, ξ = 0 respectively (24) A positive coupling strength β > 0 sources J , delaying the approach to the asymptotic solution. Negative coupling strength β < 0 drive the field towards the ghost region.
• Canonical Galileons c 2 > 0 are driven towards the trivial solution.
As Ω φ,0 → 0, a cosmological constant is necessary for these models to be viable.
• Accelerating Galileons c 2 < 0 are driven towards the tracker solution.
As Ω φ,0 > 0, accelerating models produce LUPE and can accelerate the Universe without a cosmological constant.
The approach to these solutions is determined by the relationship between the shift charge and the scalar velocity, described in figure 3. The above classification is robust against rescalings of the scalar field, which preserve the sign of c 2 . In contrast, the sign of either β, c 3 can be fixed by a field rescaling that preserves the form of the action, see Eq. (11).
Canonical and accelerating Galileons are indistinguishable at early times, either because the cubic Galileon term dominates the dynamics (|φ| |c 2 /6c 3 |) (e.g. IDEE) or because the kinetic energy is negligible (Ω φ ∼ 0). The differences occur at late times, leading to the different values ofΩ φ,0 in the asymptotic solutions (24) and will be discussed in section III D. A non-zero coupling sources the shift-charge density, driving the solution away from J → 0, Eq. (24). This is shown in figure  3 and discussed below.

Coupling & Vainshtein Mechanism
The coupling to curvature introduces a source to the shift-charge density (22) where the above expression uses the Friedman (13) and acceleration equation, andP is the scalar field pressure removing the effect of the strength of gravity (analog tô E, cf. Ref. [82,Eq. (3.5)]). The contribution of radiation and ultra-relativistic matter to the coupling is negligible since ρ rad = 3p rad . This follows from the coupling involving the Ricci scalar, which is sourced by the trace of the energy momentum tensor T ∝ 1 − 3w m . Sources to the coupling in the matter era will be discussed in section III C. Early-universe processes in the radiation era are presented in appendix C.
Analytic expressions exist for exponential coupling when the Galileon kinetic energy is negligible. For an exponential coupling β ≡ C ,φ /C is constant and the source term (25) is independent of the field value. Then the field equation (19) can be integrated directly where the kick function reads This solution accounts for the effects of M 2 * = 1 on the expansion (13) but neglectsÊ,P ∼ 0, a very good approximation at early times. Note that the kick function also affects the integrand via It is possible to decompose the solution for the shiftcharge 26 as where J 0 describes a general initial condition, J M is the contribution from the fraction of non-relativistic matter and J Σ represents the contribution from deviations from radiation domination in the early universe. The contribution from non-relativistic matter Σ M ≈ ρ mat /ρ rad = a/a eq leads to a shift-charge Appendix C describes additional sources J Σ in the early Universe. No realistic source is able to contribute significantly to the scalar field initial conditions due to the non-linear derivative interactions.
The non-canonical nature of the cubic Galileon leads to the cosmological Vainshtein screening [72,78], an efficient suppression of the coupling at early times. If the cubic Galileon term dominates, the scalar energy fraction is related to the shift-charge aŝ where Eqs. (15,20) have been used. In contrast, if the quadratic term dominates, the equivalent expression readsΩ The ratio between the energy scales associated to the cubic Galileon and canonical kinetic term iŝ so the cubic term dominates for large dimensionless field velocities. The cosmological Vainshtein screening stems from the H 0 /( |c 3 |H) factor in Eq. (31), suppressing the effects of the coupling on the shift-charge J at early times. While the cosmological Vainshtein screening may be circumvented by reducing the value of c 3 , such a reduction will incur in constraints from local gravity tests in the late universe, unless the coupling is reduced accordingly (see section V B). The effects of the screening will be shown explicitly in section III C and appendix C.

B. Imperfect Dark Energy at Equality
Imperfect dark energy at equality (IDEE) is a distinct form of early dark energy beyond GR characterized by a contribution to the expansion history that peaks around matter-radiation equality. IDEE is sourced by the cubic Galileon term, which effectively modifies gravity and changes the evolution of the perturbations (e.g. CMB). In order to affect the acoustic scale IDEE requires a significant kinetic energy of the Galileon: an initial scalar field kinetic energyΩ φ,i ∼ 10 −4 around the nucleosynthesis era evolves into a ∼ % level contribution at equality, sufficiently to reconcile early and late measurements of H 0 .
To understand the dynamics of IDEE I assume that the contribution from L 3 → c 3φ 3 H dominates the energy budget (15)Ê  FIG. 4: Imperfect Dark Energy at Equality (IDEE) in canonical uncoupled models. The initial energy density of the scalar field dilutes faster than radiation but more slowly than matter (left panel). By virtue of this scaling, the relative scalar field abundance peaks around the era of matter-radiation equality (middle panel), lowering rs and increasing H0 for fixed θ . Energy contributions of additional relativistic particles and an early quintessence model [38] are shown for comparison. The equation of state of the scalar remains in the range w φ ∈ (0, 1/3) until the kination phase at low z (right panel).
where the final expression uses the dimensionless field velocity (18). These equations can be used to set the initial condition for the field derivativeφ. Note that the energy density scaling of IDEE relies only on the domination of the cubic term. It is otherwise independent of the Galileon energy scale |c 3 |, provided that the initial field velocity is sufficiently high, as prescribed by Eq. (34). The characteristic scaling of IDEÊ . (36) follows from substituting the solution for ξ from the off-tracker evolution (19) and neglecting the coupling P φ ∼ 0, n ∝ a −3 determines the scaling of the energy density. This particular evolution, diluting faster than matter but more slowly than radiation, allows IDEE to emerge around matter-radiation equality. Figure 4 shows the scaling of IDEE for different initial conditions, along with its effects on the acoustic scale (2). Values of the initial field derivative such thatΩ φ,i ∼ 10 −4 at z = 10 10 (around the BBN epoch) grow into sizeable early dark energy contributions ∼ 5% at the epoch of equality, sufficient to lower the acoustic scale at the level needed to reconcile CMB+BAO and distance ladder inferences of the Hubble parameter. IDEE models also induce deviations from general relativity. These are best parameterized by the dimensionless braiding function [82] where the second equality applies to the limit in which the cubic Galileon dominates the energy density, Eq. (34). α B describes the kinetic mixing between the scalar field perturbations and the gravitational potentials on the cosmological background (see [83] for a covariant description). The function α B also parameterizes the deviation from the uncoupled cubic Galileon from behaving as a perfect fluid [84,85]. The last equality shows that this deviation from GR is as important as the contribution to the expansion history. The deviations from GR induced by IDEE turn out to be very restrictive for IDEE models when compared with Planck data, as I will show in section IV B.
Non-cubic covariant Galileon theories (9) dilute more slowly with the expansion, restricting their early-universe dynamics. If the quartic Galileon G 4 ∝ X 2 term dominates, its energy density scales as E 4 ∝ a wm−3 , diluting faster than matter in the radiation era and tracking the matter density afterwards. The quintic Galileon G 5 ∝ X 2 always dilutes more slowly than matter, as E 4 ∝ a 3 8 (wm−7) , corresponding to w φ = − 1 4 , − 1 8 in the radiation & matter eras respectively. Note that c 4 , c 5 = 0 may contribute to the early-universe dynamics, provided that their effect on the speed of GWs is suppressed at late times. This could happen if the field velocity kinates away, Eq. (49), soon after matter-radiation equality. While possible, this type of early modified gravity requires much more fine-tuning than the cubic Galileon implementation of IDEE.
The properties of IDEE in the simple cubic Galileon model were first discussed in Ref. [86], where it was also pointed out that the initial kinetic energy of the field would grow until the epoch of equality and could lower the acoustic scale. Previous works analyzing Galileons with general initial conditions focused on the general model [54,75,[87][88][89][90], in which the cubic and quintic terms scale faster than matter, leading to tight constraints onΩ φ,i . A more recent analysis considered the cubic Galileon separately [91], but used the same priors as in previous models and did not explicitly discuss the relevant region in which early dark energy modifies the acoustic scale.

C. Enhanced Early Gravity
Enhanced Early Gravity (EEG) consists of a time modulation of the effective Planck mass due to the scalar field dynamics and its coupling to curvature. At early times the strength of gravity is enhanced by a constnat factor, as the cosmological Vainshtein mechanism prevents any significant evolution of φ. At late times, the scalar's time variation weakens the strength of gravity, with potentially detectable signatures in local gravity and the large-scale structure of the Universe.
In EEG models, the initial effective Planck mass affects the expansion rate at early times. This effect changes the expansion rate by where M * = C(φ) and other contributions to the Galileon energy density, including IDEE, have been neglected Ω φ ∼ 0. At a fixed matter content, reducing the Planck mass M 2 * < 1 increases the expansion rate, in turn reducing the acoustic scale r s , Eq. (2). A successful EEG model requires the strength of gravity to decrease between the early and the late universe. The effective Planck mass affects all scales in the homogeneous universe, including cosmological distances, e.g. the comoving angular diameter distance Thus, if M 2 * were constant throughout, the angular diameter distance would be modified by the same multiplicative factor as the acoustic scale. This constant factor would cancel on the angular scale leaving the value of H 0 obtained from the CMB unchanged. Decreasing r s (z * ) relative to D M (z * ) requires a positive coupling β > 0. Ultimately, EEG works because the same sign of the coupling strength β required to increase H 0 drives the field away from the ghost region, cf. figure 3. The cosmological Vainshtein mechanism prevents M 2 * from evolving at early times. Assuming matter domination the shift-charge density solution (26) is where I have neglected any initial shift-charge (or equiva-lentlyΩ φ ∼ 0). Cosmological Vainshtein screening occurs when the cubic Galileon term dominates, in which case the above shift-charge density translates tȯ The scalar evolution is very suppressed compared to the characteristic evolution scale of other species, set by H H 0 . For this reason the coupling is extremely ineffective in giving the scalar field an initial velocity, as discussed in appendix C. In contrast, the unscreened regime for canonical kinetic term corresponds tȯ In that case the scalar evolves at a rate ∝ H set by cosmic expansion. Note that the above expression applies to canonical models: in accelerating models the derivative of the field is set by the non-trivial tracker solution (24).
The evolution of the scalar field leads to a running of the effective Planck mass where the second equality corresponds to the exponential coupling. Matter-domination solution in the screened regime (42) leads to a negligible running at early times, as expected.
The unscreened regime (42) for canonical models leads to a constant running of M 2 * in the matter era. α M is a standard parameterization of the impact of deviations from GR on cosmic structure formation. Just as a constant M 2 * has no effect on background observables (cf. Eq. 40), a constant M 2 * can be compensated by rescaling the abundances of all matter species so that Ω i /M 2 * is constant, leading to no net effect on the perturbations [82]. A running of the Planck mass produces deviations from GR in structure formation, potentially observable on the LSS of matter and the CMB.
Unscreened evolution (46) is expected at intermediate and low redshifts, leading to effects in LSS and secondary CMB anisotropies. Allowing α M to affect early evolution and primary CMB requires very low values of c 3 10 −9 for |c 2 | = 1 (cf. figure 14). Since this work is focused mainly on the CMB, I will set c 3 = −1 in the canonical models with Λ = 0 (in accelerating models it is set by Ω φ,0 ).
The value of the scalar field is also related to the strength of gravity measured on small scales, including the Solar System. The potential to test EEG using precision tests of GR as well as the difficulties in modeling the connection between cosmological and small scales will be discussed in Section V B. While a full investigation of these issues is beyond the scope of this work, I remind the reader that all expressions in this section refer to the cosmological evolution of the effective Planck mass.

D. Late-Universe Dynamics
The late-time dynamics of Galileons are determined mostly by the sign of the quadratic kinetic term c 2 . In accelerating models c 2 < 0 the stable solution (24) corresponds to a growingΩ φ and leads to Late-Universe Phantom Expansion (LUPE). In canonical models c 2 > 0 the stable solution (24) corresponds to a trivial configuration Ω φ → 0.
Accelerating models are very efficient at producing dark energy. The non-trivial tracker solution withΩ φ,0 > where the coupling has been neglected (β ∼ 0). With this solution the scalar kinetic energŷ rapidly dominates the energy budget (see figure 6).
No cosmological constant is needed in accelerating models. Instead, the dark energy fraction today can be obtained by choosing the ratio of c 2 , c 3 corresponding to the tracker solution in Eq. (24), corrected by contributions due to the coupling, cf. Eq. (16). Because the dark energy density grows (instead of being constant) w φ < −1, a larger value of H 0 can be obtained for fixed distance to the last-scattering surface. This is the reason why Galileon models with Λ = 0 predict a Hubble constant well above typical ΛCDM values, requiring sizeable neutrino masses m ν ∼ 0.6eV to both give a good fit and avoid a too-high value of H 0 (see appendix B). Because of their interest as DE models, the late-time dynamics of accelerating Galileons have been studied extensively in previous works, e.g. Refs. [16,51,87,91].
In canonical models the energy density of the scalar field decreases very fast once the quadratic term dominates the dynamics. This dynamical regime, known as kination, is characterized by a rapid loss of kinetic energy of the fieldφ where the coupling has been neglected (β ∼ 0). This loss of energy will continue until the coupling term becomes dominant. In uncoupled models it will evolve towards the trivial vacuumφ = 0, Ω φ,0 = 0, as anticipated in the solution Eqs. (24). Uncoupled canonical models can thus provide only negligible amounts of dark energy in the late universe, requiring an additional cosmological constant to produce acceleration. If the coupling is nonzero, the field will stabilize at a non-zero value of the shift-charge as described in section III C. Canonical models with Λ = 0 retain the freedom to set c 2 /c 3 even after using up the scalar field rescaling (11). This ratio determines onset of the kination phase, which begins when Lowering c 3 allows for a conformal coupling (β = 0) to play a role at earlier times, by weakening the cosmological Vainshtein screening (cf. section III C). Values c 3 10 −9 ensure kination occurs after recombination, and thus that the primary CMB is only affected by IDEE and EEG, as described in the above sections. For these reasons, I will set c 3 = −1 in this analysis. Some of the consquences of varying c 3 are shown in figure 14, but a more detailed study of the role of c 3 is left for future work.

IV. COSMOLOGICAL CONSTRAINTS
This section presents tests of different solutions to the H 0 problem, as implemented in coupled cubic Galileon theories. Section IV A contains an overview of the models, data and methods used. Section IV B presents the limits on IDEE and uncoupled models. Section IV C) discusses EEG in canonical coupled models. Section IV D addresses the status LUPE in accelerating models and the role of the coupling at late-times. Appendix B discusses uncoupled LUPE models.

A. Overview of Models, Datasets and Analysis
The models under study can be classified along two separate properties: • By coupling, into uncoupled β = 0 and coupled β = 0. Uncoupled models can impact r s only via IDEE, coupled models produce EEG (cf. sections III B, III C).
I will consider several combinations of models and datasets, as shown in Table I. Uncoupled, LUPE, Λ = 0 models are discussed in appendix B. IDEE is produced by the initial field velocityφ i . This is specified via a flat prior on the initial dark energy abundance log 10 (Ω φ,i ) ∈ [−8, 0], cf. Eq. (35) evaluated at z i = 10 10 (around the BBN era). The lower limit in the logarithmic prior of Ω φ,i is indistinguishable from ΛCDM, while the upper limit corresponds to the scalar field dominating the energy budget in the radiation era. The initial field velocity will be varied freely for all models presented below.
EEG relies on the initial value of the scalar φ i , which is approximately constant at early times, cf. section C. I will set φ i through a flat prior on the initial Planck , a prior on M 2 * ,i is equivalent to a prior on the initial condition for small deviations in the strength of gravity. In uncoupled models I set φ i = 0, as the initial value is irrelevant due to shift symmetry. M 2 * ,i and the coupling strength β will be varied freely for all coupled models.
The coupling strength is varied in the range β ∈ [−0.5, ∞). 2 Ghost instabilities can occur for negative coupling β < 0 (figure 3), the prior allows the data to explore that region as well. Note that β could be set instead by fixing the final effective Planck mass M 2 * ,0 . In this analysis I will not be concerned about M 2 * ,0 , deferring the issue to the discussion of local gravity tests and GW-induced instabilities in sections V B, V C.
Galileon coefficients govern the low redshift Galileon dynamics, including LUPE. In canonical models the scalar field is normalized to c 2 = +1 and the cubic coupling is fixed to c 3 = −1 to simplify the prerecombination dynamics, cf. sections III C and III D. In accelerating models the values of the Galileon coefficients c 2 , c 3 are fully fixed by normalization of the field and fix- 2 The Planck+BAO analysis of canonical models included an upper limit β ∈ [−0.5, 0.5]. This was removed in the Planck+BAO+H 0 analysis for which β ∈ [−0.5, ∞). Both analyses yield very similar bounds on β, suggesting that the more restrictive prior was broad enough, cf. section IV D.  In EEG models the coupling increases the uncertainty on H0 by a factor ∼ 3 in coupled models with Λ = 0 (canonical and accelerating) slightly increasing the central value as well. Λ = 0 accelerating models (coupled and uncoupled) predict a high value of H0, but have a bad fit and are ruled out by other observations [17]. Note that the central value of the initial effective Planck mass is M 2 * ,i ∼ 1: CMB+BAO data has no preference in the absence of late-universe information. The last line shows the best fit log-likelihood for the reference ΛCDM and differences for each Galileon: all models with Λ = 0 have a slightly better fit, while accelerating models with Λ = 0 are disfavoured.  ing the scalar field abundance today Ω φ,0 . Accelerating Λ = 0 LUPE models require sizable neutrino masses [16]. In those cases I will vary m ν ∈ (0, ∞) assuming a degenerate hierarchy. Neglecting the neutrino mass splittings has negligible differences in cosmological predictions, note that the total mass required in LUPE Λ = 0 models, m ν ≈ 0.6eV [17] is significantly larger than both the minimal mass and the mass allowed assuming ΛCDM [92,93] (see Ref. [94] for analysis of uncoupled Galileons using different hierarchies). All other cases will assume a single massive neutrino with minimal mass m ν = 0.06eV.
Other cosmological parameters were chosen following the Planck analyses [1]. I will assume the universe to have zero spatial curvature, with the fraction of scalar field energy density Ω φ given by the closure relation i Ω i = 1. The standard cosmological parameters 100θ * (or H 0 ), ω cdm , ω b , ln(10 10 A s ), n s and τ reio ∈ [0.04, ∞) are varied with flat priors unless explicitly stated. By default I will consider the Helium fraction Y He to be set by BBN given ω b and the expansion rate at early times. I will discuss constraints from light element abundances in section V A.
To test solutions to the Hubble problem, I will consider CMB data from Planck (P), distances from Baryon Acoustic Oscillations (B) and a prior on the H 0 from the SH0ES collaboration (S) in the following combinations: • Planck + BAO (PB), as the default combination. This determines the model-dependent earlyuniverse inference of H 0 and the room to accommodate late-universe measurements. PB results are summarized in table II.
• Planck + BAO + SH0ES (PBS), including a distance ladder prior on H 0 . This analysis will serve to find the global best fit. I will consider this combination in few selected cases. PBS results are summarized in table III.
The CMB data choice follows the Planck 2018 baseline analyses [1]. It includes high-temperature TT, Emode polarization EE, their cross correlation (TE) as well as low-TT and EE spectra [95]. I will not consider the Planck lensing likelihood to focus on testing primary anisotropy effects, as much as possible. Omitting CMB lensing will not significantly impact uncoupled IDEE results, which are strongly constrained by temperature and polarization alone. Coupled models can be further constrained by CMB lensing, as late-time dynamics of the scalar field will modify the lensing potential via non-zero α M , Eq. (44). This analysis will be left for future work.
BAO data is necessary for a precise inference of H 0 , anchoring r s as determined by the CMB to the late-universe expansion. For BAO data I will use the measurements from galaxy samples from the Baryon Oscillation Spectroscopic Survey (BOSS) data release 12 [96] and the low-z sample combining the 6dF survey [97] and the main galaxy sample from SDSS data release 7 [98]. I will use the galaxy BAO data as given, including density field reconstruction. This methodology is conservative for canonical uncoupled models where the late-time dynamics is indistinguishable from Λ+GR. However, the use of reconstructed data for coupled or accelerating models assumes the validity of reconstruction. This has been tested in simple extensions of ΛCDM which assume GR [99,100]. However, modified gravity can enhance nonlinear effects, including the shift of the BAO scale [101].
As direct H 0 measurement I will use the SH0ES project 2019 measurement H 0 = 74.03 ± 1.42 [km s −1 Mpc −1 ] [3], resulting in a 4.4σ tension with Planck+BAO and ΛCDM. This value relies on a distance-ladder measurement of the expansion rate with improved Cepheid variable star measurements from the Large Magellanic Cloud. The methodology has been shown to be robust by other groups [102][103][104][105]. Other late-universe measurements of the Hubble parameter exist tend to produce larger values of H 0 than Planck+BAO within ΛCDM that are either in tension (lensing time delays [4]) or compatible (standard sirens [106], tip of the red-giant branch [107]), see Ref. [7] and [1, section 5.4] for recent overviews. Adding a prior on H 0 serves to find the best-case scenario and its goodness-of-fit in light of all available (although possibly discrepant) datasets.
Type Ia SNe data will not be included in this analysis, but left for future work. In coupled Galileons the interpretation of SNe data requires modeling the variable strength of gravity on small scales and its effect on the intrinsic SNe luminosity, as discussed in section V B 2. Note that a time-variation of SNe luminosity invalidates the inverse standard ladder method (BAO+SNe) of inferring the acoustic scale (figure 1). SNe modelling issues are absent in uncoupled Galileons, but SNe will not qualitatively change the conclusions of this analysis. In canonical uncoupled models (Λ = 0) the expansion history is indistinguishable from ΛCDM at low redshift. Accelerating uncoupled models (LUPE) are disfavoured by SNe, but the tension can be read directly by comparing the contours with the inverse-distance ladder (BAO+SNe) in figure 1.
To obtain the theoretical predictions I used the hi class code 3 [48,49,81], where the exponentially coupled cubic Galileon model (section II A) was implemented using the covariant Lagrangian approach developed in version 2.0 (see Ref. [49] for details). The parameter space of each model and dataset combination was sampled using a Markov chain Monte Carlo (MCMC) analysis with a Metropolis-Hastings proposal distribution. The sampling relied MontePython (version 3) [108,109], modified to record errors whenever model predictions can not be computed, such as unstable regions of the parameter space. To ensure convergence the MCMC runs until the variance across chains over in-chain variance (Gellman-Rubin convergence ratio) is smaller than 0.05. The resulting chains were analyzed with MontePython, Getdist [110] and CosmoSlik [111].  have a period of fast growth that affects the CMB spectrum in a characteristic scale-dependent manner.
The initial fraction of dark energy Ω φ,i is constrained by Planck+BAO to the point where its effect the acoustic scale and H 0 is negligible. Figure 7 shows the posteriors marginalized over the initial energy density and the Hub-ble parameter for uncoupled models. The relationship Ω φ,i − H 0 in the absence of constraints is shown for fixed θ * and other cosmological parameters. Using IDEE to solve the Hubble problem would require Ω φ,i 10 −4.2 in the canonical model (cf. section III B), while CMB+BAO bounds are at the level of Ω φ,i 10 −5.2 at 95% c.l. (see table I).
The bounds on IDEE make canonical uncoupled models indistinguishable from ΛCDM. The impact of IDEE in the late universe is bound to be smaller than on the primary CMB, due to the IDEE scaling in the matter era and the late kination phase (cf. figure 6). As a consequence, the parameter bounds are almost identical to those of the reference ΛCDM analysis. The role of IDEE is also negligible on accelerating uncoupled models. Those cases are indistinguishable from settingφ i to the tracker value, Eq. (24). The fit best fit likelihood to Planck+BAO is significantly worse than the canonical case. Moreover, other analyses rule out the model including CMB×LSS cross-correlations [17,112] and a combination of late-universe datasets [94].
Including a distance ladder prior on H 0 in does not alter these conclusions significantly. The bound on IDEE becomes slightly higher log 10 (Ω φ,i ) < −4.94 95% c.l. for the canonical model. Trying to fit datasets in tension leads to larger shifts on the remaining cosmological models, with a change ∆H 0 /σ H0 = 1.51 driven mainly by ∆n s /σ ns = 1.05, ∆ω b /σ ω b = 0.80, ∆ω cdm /σ ω cdm = −0.82, as can be seen comparing tables II and III.
An autopsy of IDEE shows that the strong limits on Ω φ,i originate from the growth of the scalar field perturbations around horizon crossing. This can be understood by examining the mass-squared for the field fluctuations where c s is the scalar sound speed and the time and scale-dependence of the mass is shown in the left panel of figure 8. A consequence of cubic Galileon domination is that the scale-independent contribution is negative m 2 φ < 0, a feature known as tachyon instability. Tachyons are associated with growing scalar field perturbation V X ≡ δφ φ ∼ e ±µt (i.e. imaginary frequency) on scales larger than the scalar field sound horizon k < c s /(aH). For perturbations at a scale k, the growth begins around horizon crossing. The rate of tachyonic growth is modulated by |m 2 φ |. This is proportional to Ω φ and thus enhances the scalar field perturbations significantly before recombination in models able to affect r s .
The growth of scalar perturbations is tamed by the scalar-field pressure, i.e. the scale-dependent term in the effective mass (51). Stability on small scales requires c 2 s > 0, or equivalently, that small-scale perturbations in the field have oscillatory solutions. This oscillatory regime begins once the sound speed dominates the effective mass, with scalar field perturbations decaying by virtue of the Hubble friction. This transition corresponds to the red/blue border in the left panel of figure 8. Different physical scales undergo growth and oscillations at different times, leading to different amplitudes at recombination ( figure 8, right panel). The scale dependence of the field perturbations is transferred to the gravitational potentials and other species via the modified Einstein's equations.
The interplay between growth and oscillations leave a characteristic imprint on the CMB spectra. The tachyonic growth is largest for modes that enter the horizon soon before recombination and correspond to the first peaks and troughs of the CMB spectra (figure 8). These modes have no time for the oscillatory phase to stabilize their growth, leading to a larger impact on relatively low multipoles in the CMB, as shown in figure 9. The differences are strongest for the TT spectrum on the larger angular scales ( 1000), particularly on the first peak and through. Overall there is an enhancement of the odd peaks (1, 3 and 5) and a suppression of the even peaks (2 and 4). However, this effect can not be compensated adjusting the value ω b . This odd/even pattern is overlaid with an overall suppression of intermediate angular scales (1000 2000), and an enhancement of small angular scales 2000. The EE polarization spectrum (right panel in figure 9) shows a similar trend, with deviations becoming smaller in higher multipoles.

C. Coupled Models: Viability of EEG
Now I will discuss EEG models, focusing on the initial effective Planck mass M 2 * ,i , its effect before recombination and its impact on the Hubble rate and other cosmological parameters. I will also describe the main features of EEG and the differences to IDEE and other early dark energy scenarios. Due to the Vainshtein mechanism both canonical and accelerating models have the same early-time behavior. For this reason, I will focus on canonical models and leave the discussion of both accelerating models and constraints on the coupling β for section IV D below.
A non-zero coupling introduces a significant degeneracy between the initial effective Planck mass and the Hubble parameter with the potential to accommodate high values compatible with late-universe constraints. Figure 10 shows the marginalized posteriors on the M 2 * ,i -H 0 plane, exhibiting the anti-correlation between both quantities, as anticipated in figure 1 and described in section III C. This relation can be understood as follows: a weakening of gravity at early times M 2 * ,i < 1 increases the early expansion rate before recombination and reduces the acoustic scale. Then the same projected CMB scales correspond to a larger value of H 0 , as long as the late time effective Planck mass is larger than the prerecombination value (e.g. today M 2 * ,0 < M 2 * ,i ), which requires a positive coupling constant β > 0. While the effect of the initial effective Planck mass is the basis of EEG other parameter degeneracies also play a role in relieving the Hubble tension.
EEG introduces new degeneracies with cosmological parameters, contributing to accommodate larger values of H 0 in EEG with respect to ΛCDM and IDEE. These degeneracies and the resulting enlarged posteriors are apparent from a triangle plot, figure 11. For the purpose of the Hubble problem, the most important is the anticorrelation between M 2 * ,i and both the baryon density ω b and spectral index n s . As both ω b , n s are themselves correlated with H 0 , decreasing M 2 * ,i leads to a higher Hubble rate by virtue of these degeneracies. These effects are on top of the direct reduction in the acoustic scale caused by EEG. There is an additional, mild anti-correlation between M 2 * ,i and the amplitude of perturbations σ 8 (or equivalent A s ). Interestingly, the dark matter abundance ω cdm has no apparent correlation with M 2 * ,i , although it correlates weakly with the coupling β. The introduction of a coupling increases the limits on IDEE Ω φ,i only slightly and does not allow it to play any role on the Hubble tension (cf. table II).
The anti-correlation between ω b and M 2 * ,i is driven by the BBN relation between the baryon and Helium abundances assumed in the analysis (cf. figure 11, upper right panel). This relation limits the damping scale by linking the helium fraction to ω b and M 2 * ,i . The degeneracy is analog to bounds on additional relativistic species (cf. [1, figure 39]), as both increasing N eff or decreasing M 2 * ,i lead to a faster expansion rate in the BBN era, although the CMB is independently sensitive to relativistic species via perturbations (see section V A). Lifting the standard BBN assumption will thus increase the range of allowed values for M 2 * ,i and H 0 , but considering the measured primordial helium abundances will limit this range, cf. [1, figure 41].
Note that for a constant M 2 * the equations for the cosmological expansion and linear perturbations depend only on the ratioω i ≡ ω i /M 2 * of the different matter components [82]. Therefore, the predictions remain unchanged if all physical densities ω i are rescaled, leaving ω i invariant. This might suggest a correlation ω b ∝ M 2 * , very different from the anti-correlation observed in the data, roughly ω b ∝ M −1/2 * ,i . This apparent contradiction can be explained by noting that 1) M 2 * is not constant, with M * ,0 > M * ,i due to the coupling and 2) unlike ω b , ω cdm , the radiation density is fixed by the CMB temperature today and can not be rescaled. Moreover, as discussed above the baryon fraction degeneracy is set by standard BBN and the effect of helium on the CMB damping.
While Planck+BAO constraints on EEG models allow a ∼ 3× wider range of values for H 0 , they show no significant preference towards either lower or higher values relative to ΛCDM. Including a prior from distance ladder shifts the posteriors towards high values of H 0 . This analysis results in a ∼ 2σ preference towards stronger gravity at early times M 2 relative to the Planck+BAO case. The shift of other cosmological parameters when including the H 0 prior follows the same trends as for uncoupled IDEE models, with the dominant shifts being ∆ω b /σ ω b = 0.93, ∆ω cdm /σ ω cdm = −0.43, ∆n s /σ ns = 0.97, ∆σ 8 /σ σ8 = 0.50 and ∆β/σ β = 0.51, leading to change in the Hubble rate by ∆H 0 /σ H0 = 1.37. While the relative shifts caused by the H 0 prior on cosmological parameters are similar as in uncoupled IDEE models, the larger uncertainties of coupled models lead to a stronger net shift.
The main effect of enhanced early gravity is to lower the amplitude of the CMB spectra when M 2 * ,i < 1. Figure 12 shows the impact of M 2 * ,i on the temperature and polarization power for fixed cosmological parameters and choosing the coupling β so the effective Planck mass today is M * ,0 = 1 (this is not assumed in the MCMC analysis). The lower temperature power is caused by an enhanced Sachs-Wolfe effect: stronger gravity deepens the gravitational potentials, increasing the redshift of photons emitted from overdense regions. The dependence of this effect on the angular scale is mild, allowing small shifts in cosmological parameters to partially compensate for the differences. These degeneracies are not accounted for in the solid lines of figure 12, leading to a seemingly worse fit than if other cosmological parameters had been varied.
EEG models present important differences relative to IDEE and canonical models for early dark energy. In EEG models the scalar field modulates the strength of gravity, but because of the Vainshtein mechanism the value of the field remains approximately constant in the early universe, cf. section III C. This is equivalent to a constant energy density contribution, which does not affect the ratio of energy densities of all matter and radi- ation species before recombination. This is in sharp contrast with both IDEE and quintessence models of early dark energy in which scalar energy density evolves before recombination (cf. figure 4). The constancy of the scalar field before recombination also prevents deviations from GR to affect the perturbations (i.e. α M ∼ α B ∼ 0), whose dynamics is the same as in standard cosmology but with abundances rescaled by M 2 * [82]. In contrast, IDEE models induce deviations from GR proportional to Ω φ , including the tachyonic growth described in section IV B.

D. Late-time dynamics of coupled models
Let us now examine the late-time dynamics of coupled cubic Galileons. I will discuss the constraints on the coupling and the status of canonical (EEG) and accelerating (LUPE+EEG) models. Uncoupled, LUPE-only models are discussed in appendix B, including the role of Λ and m ν . The constraints on the coupling strength β and the initial effective Planck mass M 2 * ,i are shown in figure  13.
The effect of M 2 * ,i , described in the previous subsection, is very similar across all coupled models. The preferred values of M 2 * ,i depend only mildly on the model, although including a distance-ladder prior on H 0 shows a preference for M 2 * ,i < 1 corresponding to EEG cf. figure  10. There is a significant widening of the H 0 posteriors due to EEG and the parameter degeneracies already discussed in the previous subsection. The main difference is the central value of the Hubble parameter, which is sensitive to the late-universe expansion and differs in accelerating models via LUPE. In the case of Λ = 0 accelerating Galileons, that central value is much closer to the distance ladder measurement than in canonical models with Λ. The accelerating model with Λ = 0 is an intermediate case between the two.
The coupling is constrained by stability criteria and the late-time evolution of the model. Negative values are mostly excluded as they drive the field evolution towards a ghost instability (cf. section III A). Very small negative values may be supported by the initial field velocity, but this is related to the initial energy density Ω φ,i and very limited by the analysis of IDEE models (section IV B). The Vainshtein mechanism prevents β from playing any role before recombination (cf. section III C). Therefore, the coupling is constrained by late-universe physics, including low redshift expansion history and secondary CMB anisotropies (ISW effect, CMB lensing). As late-time dynamics depends greatly on the presence of Λ and the accelerating or canonical nature of the model, each sub-class has different limits on β, as evident from figure 13.
The strongest constraints on β occur in accelerating models with Λ = 0, where the field time derivative is largest. The absence of a cosmological constant requires a largeφ at late times to support Ω φ,0 ≈ 0.7, Eq. (24). This variation translates on a sizeable running of the effective Planck mass α M ∝ βφ, Eq. (44), which is severely constrained by the ISW effect's impact on the CMB's large angular scales. Thus, a coupling is very constrained and tends to exacerbate the problems of accelerating Galileons. Coupled models fare no better than the uncoupled version. They are disfavoured by Planck+BAO (the best-fit likelihood is even worse for the coupled model, despite being an extension cf.  the uncoupled version they are also strongly ruled out by other observations, such as LSS×CMB cross-correlations [17]. The exponential coupling can thus not save accelerating Λ = 0 Galileons, offering no solution to the H 0 problem. Accelerating models with a cosmological constant have milder bounds on β. Allowing for Λ = 0 eliminates the burden of cosmic acceleration from the Galileon field, which becomes a sub-dominant contribution to the energy density. The velocity of the field is not tied anymore to the requirement of cosmic acceleration and can be lowered significantly. Note that Table II shows Ω φ,0 ≈ 0.11, but this includes the contribution from the effective Planck mass today M 2 * ,0 , which does not contribute to eitherφ 0 or α M . Coupled accelerating models with nonzero Λ give a reasonably good fit to CMB+BAO, while increasing the allowed value of H 0 both due to EEG and sub-dominant LUPE contribution. Other cosmological data (such as Ia SNe and LSS×CMB) may place further limits on this scenario.
Canonical models have the loosest constraints on β. In this case, cosmic acceleration is entirely supported by the cosmological constant, and thus the contribution from the scalar field to the energy density can be arbitrarily small. Planck+BAO prefer a very subdominant contribution Ω φ,0 ∼ 0.02, which is further split into the strength of gravity ∝ M 2 * ,0 − 1 and the kinetic termŝ E ∝φ, cf. Eq. (15). The kinetic termsÊ are typically small, as the field derivative is sourced by the coupling ∝ βρ mat ∝ a −3 , which reduces the value of α M cf. Eq. (46). In contrast, accelerating models are driven by the non-trivial solution, Eq. (24), associated with larger field derivatives.
While this analysis has focused on Planck, BAO and SH0ES, all coupled models can be further constrained by additional cosmological probes and tests of gravity. In the next section I will outline some remaining challenges for coupled models, including big-bang nucleosynthesis, precision tests of gravity and gravitational waves.

V. CHALLENGES FOR COUPLED GALILEONS
In this section I describe further observational constraints that may challenge coupled models implementing Enhanced Early Gravity and/or Late Universe Phantom Expansion. I will first discuss the effect of the effective Planck mass on primordial nucleosynthesis (section V A). Then I will address the issue of local tests of gravity, including scalar fifth forces, the value of the Planck mass and its time variation (section V B). Finally, I will discuss how GWs may induce instabilities in the scalar perturbations, pushing the theory beyond its regime of validity (section V C).

A. Primordial Element Abundances
The primordial abundance of light elements is sensitive to the expansion rate in the era of Big-bang nucleosynthesis (BBN). It can be used to place constraints on the initial effective Planck mass M 2 * ,i independent of the CMB. I will explain how to translate known BBN limits on the expansion history to EEG models and discuss their implications.
BBN limits are often quoted in terms of additional light particles, such as the number of neutrino species N ν . By comparing the Hubble law in coupled models (13) to the effects of additional radiation in the expansion history, one can derive a relation between the initial effective Planck mass and extra radiation 1 where the second equality uses ρ = π 2 30 2 + 7 2 + 7 4 N ν accounting for photons, electrons and neutrinos active in the BBN era [113]. Note that ∆N ν is defined relative to a fiducial value N ν = 3, neglecting the small correction from the energy injected by positron annihilation that lead to the difference between N ν and the more widely used N eff [114]. Note that the equivalence between N ν and M 2 * ,i can be applied to nucleosynthesis constraints because BBN is sensitive to additional components only through the expansion history. In contrast, CMB anisotropies are sensitive to perturbations in the additional species, including a phase shift due to the super-sonic propagation of neutrinos [115][116][117].
Limits on the initial effective Planck mass from primordial abundances of deuterium and helium can be translated using Eq. (52) using no CMB data. Because the CMB responds differently to relativistic particles and modified gravity, I will only quote values not involving any input from Planck (see Ref. [1, section 7.6] for the Planck implications on BBN and N eff ). From more to less conservative, several 95% c.l. limits with no CMB information are • M 2 * ,i > 0.860 (∆N ν < 1) for helium only [113,118], • M 2 * ,i > 0.911 (∆N ν < 0.6) including the degeneracy with ω b , [113, Fig. 10] and • M 2 * ,i > 0.939 (∆N ν < 0.4) marginalized over ω b , which is the value quoted in the most recent review of particle properties [119].
Besides avoiding CMB data, the above bounds assume that the only effect on the expansion history is from a constant effective Planck mass. This is an excellent assumption in IDEE models, for which Ω φ (z BBN ) 10 −5 (Planck+BAO limits) and even Ω φ (z BBN ) ∼ 10 −4 is required to reconcile H 0 values. The assumption remains valid in the presence of a coupling thanks to the cosmological Vainshtein screening, which prevents the scalar field to vary significantly at early times (section C). BBN constraints can be more stringent for coupled theories without cosmological screening, as the effective Planck mass can vary during the BBN era [122][123][124][125].
It is worth emphasizing that the BBN predictions have been included in the CMB+BAO constraints on EEG , and play an important role by relating the helium fraction, ω b and M 2 * ,i . As discussed in section IV C and emphasized in figure 11 (top right), lifting the assumption of standard BBN will lead to looser constraints. In that case, including bounds on helium and deuterium abundances will be particularly important to supplement Planck data (see Ref. [1, figure 41] for the case of additional relativistic particles). While current CMB+BAO places more stringent limits than BBN on M 2 * ,i and EEG, any improvement on the measured primordial abundances can be used as a further test.

B. Gravity on Small Scales
Deviations from Einstein's general relativity are very well constrained by local gravity tests. There are at least three effects that may be used to constrain coupled Galileons and limit EEG and LUPE solutions to the Hubble problem: • Scalar forces • Local strength of gravity • Time-variation of the local strength of gravity Reliable constraints based on these effects require solutions of the coupled Galileon theory connecting the cosmological solution to very small-scales. In this section I will discuss the challenges to connect cosmological and local dynamics of coupled Galileons, limits on the above effects from lunar laser ranging (LLR) and other precision gravity tests and the challenges in interpreting type Ia supernovae (SNe) observations in coupled models.

Scalar Forces
The Galileon scalar field mediates an additional, attractive interaction. While the scalar force is very suppressed within the Vainshtein radius [66] it leads to a small deviation form the 1/r 2 dependence of the gravitational force that causes a small shift in the orbital phase of bound objects and can be probed by sensitive enough measurements (e.g. [126]). Lunar laser ranging (LLR) measurements set the following bounds for the phase shift of the Moon to be within the observational limits [127]. Note that screening implies that the coefficient of the cubic kinetic term c 2 does not enter this bound. The upper bound (54) leads to β 2.15 for the canonical EEG models studied here (fixed c 2 = 1, c 3 = −1), a value well well above the cosmological bounds, cf. figure  13. The above bounds implicitly assume i.e. the Moon's orbit is confined within the Vainshtein radius of the Earth (53) The lower limit on the coupling (55) is only relevant only if c 3 ∼ 0. Even in the lack of screening |c 3 |/c 2 2 → 0, a change in the gravitational strength is proportional to and thus negligible if β is sufficiently small. An estimate for the bounds in the uncoupled regime (56) can be obtained by comparing the coupled free theory to the Brans-Dicke Lagrangian [128, section 3.1], where one can identify β ∼ 1, c 2 ∼ ω BD 4 · 10 4 and the lower limit is required for compliance with Solar System tests. Comparing the theory-dependent coefficient of Eq. (56) with the Brans-Dicke case suggests that β 2 · 10 −2 √ c 2 is in agreement with observations even in the lack of screening.
The above limits on the scalar force follow from an expansion around the Minkowski solution. However, the time evolution of the scalar field modifies the Galileon terms. These corrections have been computed only on de-Sitter backgrounds in which the field evolution is stationaryφ = constant, cf Ref. [15, appendix B]. For cubic Galileon only the quadratic term is affected. This can be seen from expanding the action for the total field locally φ loc = φ + ϕ as where φ = φ(t) is the cosmological solution and ϕ = ϕ(t, x) is the local correction. Since cosmological evolution is slow compared to dynamical time-scales of the Solar System, the term in brackets can be taken as constant and the difference between the connections δΓ ∼ H can be neglected. Then, cosmological evolution amounts to a redefinition of c 2 , which does not affect the constraints on the scalar force in the screened regime (54). The local time-variation will be most decisive on the coupling function, as it affects the strength of gravity measured in the Solar System.

Local Strength of Gravity & Supernovae
The scalar field coupling modulates the local value of the Newton's constant, which depends on the local value where φ = φ(t) is the cosmological solution and ϕ = ϕ(t, x) is a local correction. For exponential couplings the measured value is recovered if βφ loc ≈ 0 in the Solar System today. Because of the space-time dependence of the field, this condition does not necessarily reduce to fixing M 2 * ,0 = 1 (or φ(t 0 ) = 0) on the cosmological solution. A detailed calculation of φ loc needs to account for vastly different scales, including how the cosmological solutions adapt to the local dark matter halo, how that solution adapts to the galaxy, and so forth, all the way to the Solar System. In addition, it is necessary to model the evolution of the scalar field over the timescales in which those structures form. While such an analysis is well beyond the scope of this work, I will discuss possible outcomes for the local solution. 4 Shift-symmetry φ loc → φ loc + C guarantees the existence of solutions where the field evolves at the cosmological rate around a matter source. For a spherically symmetric configuration such a stationary solution takes the form [133] φ loc,s =φ(t 0 ) · t + ϕ s (r) , where the only differences with Eq. (58) is that the field velocityφ is constant and the local correction is staticφ s = 0. Introducing the above ansatz in the dynamical equations explicitly neglects field accelerations φ loc ,φ, losing any information about how that solution is reached. Stationary solutions (59) are a likely endpoint for the dynamical evolution near a matter source, once the local value values of the scalar field reaches equilibrium with the cosmological evolution. The main question is whether this occurs before the present time or in the cosmological future. Two scenarios are possible depending on the relation between the local and cosmological evolution: 1. Homogeneous evolution: ifφ loc ≈φ in the Solar System, then M 2 * ,0 ≈ 1 is a necessary condition. In this case α M is constrained directly by the variation of Newton's constant (section V B 3).
2. Inhomogeneous evolution: ifφ loc φ then φ loc (t 0 , 0) ∼ 0 requires M 2 * ,0 > 1. The latter condition is compatible with EEG, but requires a sizeable value of the coupling β, which enhances the scalar force (54). In this caseφ loc (t 0 , x 0 ) could be small enough to satisfy bounds on the time variation of Newton's constant (section V B 3), but cosmological effects will be larger.
The caseφ loc φ is both nonviable and unlikely. Inviable because it would yield a large time variation of Newton's constant. And unlikely because the Vainshtein mechanism slows down the field evolution in screened regions.
The analysis of the time-dependent Galileon equation in a screened region suggest that inhomogeneous evolution could happen in small scales. The starting point is a spherically symmetric field configuration in Minkowski space, where the field evolution is governed by [131] 4 See Ref. [129] for a detailed analysis of this issue in non-local gravity theories without the Vainshtein mechanism and Ref. [130] for a study of the interplay between local solutions and cosmological time dependence in cubic Horndeski theories. The time-evolution and stability of spherically symmetric systems approaching the Vainshtein screened solution was studied in Ref. [131]. Ref. [132] finds a suppression of the local field velocity in Chameleon models.
where primes denote radial derivatives, t, r are in units of M = (H 2 0 M p ) 1/3 andφ(r, t) ∝ φ loc , ρ have been made dimensionless. The kinetic coefficient is given by and screened region is characterized by Z tt c 2 . While one does not expect the field to evolve in a strictly static space-time, this simple configuration was used to model the approach to the static, Vainshtein-screened solution. A similar analysis might shed light on the interplay between the small-scale and the cosmological solution, which would appear as a boundary condition in this situation.
The evolution timescale in a screened region can be estimated evaluating Z tt on the static screened solution. Forφ Thus the characteristic timescale for the field evolution is rescaled by a factor ∼ √ Z tt ∝ (r/r V ) 3/4 , slowing down evolution in screened regions. This is analog to the cosmological screening mechanism discussed in section III C.
SNe observations need to be reinterpreted in coupled Galileons to reflect the variable strength of gravity. The intrinsic luminosity of SNe is expected to depend on the Chandrasekhar mass, the threshold for a white dwarf to be supported by electron degeneracy pressure. Its dependence on the strength of gravity implies that the intrinsic luminosity of a SNe will vary with φ loc . Early works on the subject argued that stronger gravity (lower M ch ) leads to dimmer SNe (lower ejected mass) [134][135][136]. However, a more recent study based on a semi-analytical model for SNe light curves concludes that the opposite is true, with stronger gravity producing brighter SNe after standardization [137]. 5 If the scalar field evolves homogeneously, the variation of M 2 * produces a redshift-dependent correction to the luminosity distance observed by SNe (this can be tested even independently of the specific model [140]). If the evolution is inhomogeneous, the SNe luminosity will also depend on the properties of the host galaxy/halo, leading to an additional scatter in the Hubble diagram. If the scatter is significant, it can be probed by methods used to search for lensing signatures of compact dark matter [141].
Any other observation that may rely on the strength of gravity needs to be reinterpreted in coupled models. One example is the measurement of the Hubble rate inferred from lensing time delays, for which the Hubble parameter inferred from a lens at redshift z L scales as H 0,z L = H 0,true /M 2 * (z L ) [142]. GW standard sirens observations need to be reinterpreted along similar lines once they become available over a larger redshift range. Unlike for SNe, in both gravitational lensing and GWs the relationship between the luminosity and the strength of gravity is well understood, and the only challenge is modelling the connection between the cosmological evolution and the relevant scales.

Time Variation of Newton's Constant
The variation of the local scalar field value is equivalent to a time-varying Newton's constant, an effect that can be constrained via precision tests of gravity. The most precise current bound is based on LLR are [143]Ġ N G N = (7.1 ± 7.6) · 10 −14 yr −1 , or equivalently [129] C Cφ where it has been assumed thatG N = 0 and the result is quoted in terms of the scalar field variation using G N ∝ C(φ loc ) −1 . The time variation of Newton's constant is strongly correlated with the vector describing the rotation of the Moon's core [143, section 4.1]. Since the core rotation is poorly constrained independently, the above limits assume the core rotation vector obtained from the standard withĠ N = 0. More conservative assumptions about the Moon's inner structure might lead to weaker constraints.
Note that the central value ofĠ N /G N corresponds to a growing strength of gravity. In contrast, coupled cubic Galileons predict a decreaseĠ N < 0. This is a theoretical an observational requirement. Theoretically, it follows from the need fo positive coupling constant β > 0, necessary required both to prevent ghosts (cf. III A). Observationally it is required for solving the Hubble problem via EEG and the need to increase the strength of gravity.
The impact on the model parameters and the H 0 tension requires understanding the connection between the global and local dynamics of the scalar. If the scalar field evolves homogeneously, the bound on the variation of the gravitational constant (64) can be translated directly into a stringent limit on the effective Planck mass running today [79] βφ Comparison with the approximate expressions in section III C indicate that only very small values of the coupling would be allowed, ruling out EEG.
Reducing the cubic coupling c 3 reduces α M for fixed value of today's effective Planck mass M 2 * ,0 = 1. This slowdown works by weakening the cosmological Vainshtein mechanism, allowing the field to start evolving earlier, as shown in figure 14 for models where β is adjusted so M 2 * ,0 = 1. Values c 3 10 −9 correspond to evolution before recombination, potentially impact primary CMB anisotropies. The slowdown achievable is not enough to prevent LLR constraints on Newton's constant variation, at least assuming a sizeable EEG (M 2 * ,i ∼ 0.95, assuming homogeneous evolution (65) and assuming the latest result with standard value of the Moon's core rotation (64). In addition, reducing the cubic coupling makes Galileons vulnerable to If the local evolution is inhomogeneous the limit needs to be satisfied for the value of the scalar field in the Solar system. A tentative order of magnitude estimate of the Vainshtein suppression suggest applying the constraint (64) to the Planck mass dressed by the kinetic term as discussed in section V B 2 If the above scaling holds, the dependence of the kinetic term with the radius in a screened region (62) indicates that the time variation of the Newton's constant could be very suppressed locally, allowing EEG to remain compatible with time-variation of Newton's constant.
The back-of-the envelope slowdown in screened regions (66) is likely an overestimation. While better modeling is needed, the true solution is likely to lie within the two limits, Eqs. (65,66). The homogeneous case is well beyond the limit (64) for the EEG , models in which the effective Planck mass evolves significantly M 2 * ,i ∼ 0.95 → M 2 * ,0 ≈ 1, since that evolution occurs mostly at low redsfhit. A very efficient suppression oḟ φ loc implies that a large coupling β is required to connect EEG at early times to the correct local value φ loc → 1 today. Large β would be problematic for both cosmology ( figure 13) and scalar force constraints Eq. (54).

C. Gravitational Waves
Coupled cubic Galileon Gravity avoids constraints from the GW speed and decay by construction. In this section I will discuss other GW tests of coupled cubic Galileons, focusing on the scalar instabilities induced by passing GW.
Cubic Galileon interactions may induce instabilities in scalar sector: a background GW can flip the sign of the kinetic term for scalar-field perturbations triggering a ghost or gradient instability [71]. This requires a GW with sufficient amplitude propagating on a non-screened region, which is estimated to have happened in a significant fraction of the universe unless where the cubic Galileon term contribution to the braiding is given by Eq. (37). The relevant quantity above is the the contribution of the cubic term to the braiding α B (the coupling also contributes to α B , but not to the instability). Whether a given model triggers the instability depends on the time variation of the field at low redshift. Models in which the field evolves rapidly are most susceptible to the instability. Figure 15 shows examples selected from the best-fit models in the cosmological analysis (cf. section IV). The instability is triggered in all accelerating models unless a cosmological constant is allowed and the contribution of the scalar field to the energy density is very subdominant. In canonical Λ = 0 models the field variation at late times is driven by the coupling. The instability is triggered only for values of β that are sizeable, yet allowed by cosmology (cf. figure 13). In canonical uncoupled models the field's kinetic energy dilutes very fast at late times and remains well below the unstable region. Note that any feasible amount of IDEE can not trigger the instability. This is due both to the stringent bounds from CMB and the fact that GWs sources with enough amplitude exist only at relatively low redshift, after IDEE peaks. While including limits from the instability improves over the CMB+BAO constraints, these improvements are milder than suggested by studies based on parameterizations of the alpha-functions [144].
While taking the GW-induced instability as a hard constraint is complementary to the cosmological analysis, it is important to remember that they are conservative  (37) for some best-fit models resulting from the analysis of section IV). Accelerating (magenta) or canonical models with sizeable couplings (dark cyan dashed) are able to trigger the instability at late times (shaded region), Eq. 67, where it has been assumed that no GW sources exist for z > 30. The role of the coupling β is shown for a canonical model (dark cyan) with the best fit-value (solid) and a value close to the excluded region (dashed). Accelerating Λ = 0 models produce even larger values αB,3(t0) ∼ 1 (not shown).
from a theoretical point of view. The fate of the theory after the instability is reached is uncertain. Specifically, it is not clear whether the instability is associated with any prediction which violates current experimental bounds, and simple models exist in which a similar instability is associated with no pathological behaviour [71]. All that can be said for sure is that a high-energy completion of the theory is needed to address the consequences of entering the unstable region.
Coupled Galileons also predict a mismatch between distances measured from GWs (standard sirens) and electromagnetic or geometric observations (e.g. SNe, BAO). This difference is produced by the effect of the conformal coupling G 4,φ = 0 on the GW propagation. Current bounds are very weak | C Cφ H0 | O(10) [145], well below the level of other probes discussed here. Upcoming GW observation campaigns and new detectors will improve these limits considerably [44,146,147]. However, it has been argued that the interpretation of standard sirens needs to be reconsidered in theories with screening mechanisms [148][149][150] (see also Ref. [151]). Because standard sirens are not yet a competitive test, I will not discuss them further.

VI. CONCLUSIONS
Discrepancies in the Hubble constant inferred by different methods could be an indication of physics beyond the simple ΛCDM model and its underlying as-sumptions. Here I have examined three different mechanisms by which gravity theories beyond Einstein's GR may alleviate the discrepancy on the H 0 values inferred via BAO+CMB and distance ladder observations. Imperfect Dark Energy at Equality (IDEE) and Enhanced Early Gravity (EEG) modify the pre-recombination expansion history to reduce the acoustic scale r s for fixed angular projection θ * . Late-Universe Phantom Acceleration (LUPE) is based on the dark energy density growing at low redshift w φ < −1. Each mechanism can operate individually or in combination with the others.
The three mechanisms exist in the coupled cubic Galileon, a simple scalar-tensor theory compatible with the speed of GWs, lack of GW decay and equipped with the Vainshtein screening mechanism. This investigation focused on an exponential form of the coupling G 4 = C(φ) ∝ e βφ , and considers the two possible signs of the quadratic kinetic term ∝ c 2 (∂φ) 2 , dubbed canonical (c 2 > 0) and accelerating (c 2 < 0). Different combinations of model properties (coupled/uncoupled × canonical/accelerating) were tested against Planck+BAO, including in some cases the SH0ES distance ladder measurement of H 0 to address the tension in extended models.
The main findings regarding the cosmology of these models can be summarized as follow: 1. IDEE relies on the scaling of the scalar field energy density, with dilutes faster than matter but more slowly than radiation ( figure 4). It requires a large initial velocity for the fieldφ i . Values corresponding to Ω φ (z BBN ) ∼ 10 −4 would lower r s enough to reconcile CMB+BAO with SH0ES for fixed θ * .
2. Planck+BAO constrain IDEE to Ω φ (z BBN ) 10 −5 , below the level necessary to solve the H 0 problem ( figure 7). The strong bounds on IDEE stem from modified gravity and expansion. A tachyon instability enhances the growth of scalarfield perturbations after Hubble crossing, impacting mainly the first CMB peaks (figures 8 and 9).
3. EEG relies on the coupling C(φ)R, allowing the scalar field to modulate the strength of gravity via the effective Planck mass M 2 * = G/G eff = C(φ) = exp(βφ). EEG requires for the field to roll in the late-universe to reduce r s at fixed θ * (figure 5). Initial conditons corresponding to M 2 * ,i ∼ 0.95 evolving to M 2 * ,0 = 1 could solve the H 0 problem.
The coupling strength is severely restricted in the Λ = 0 case to |β| < 0.05 95% c.l., preventing a coupling from improving the fit for the accelerating cubic Galileon ( figure 13). Λ = 0 models reduce H 0 tension to 2.5σ via a combination of LUPE and EEG.
It is remarkable that modified gravity solutions to the Hubble problem require far less fine-tuned initial conditions than other early dark energy models. IDEE stems from the initial field velocity and scales only mildly with the dominant matter component. EEG stems from the initial field value and its contribution remains constant at early times thanks to the cosmological Vainshtein mechanism. Generic initial conditions of the field produce some amount of IDEE and EEG in coupled cubic Galileons.
To reconcile H 0 , EEG requires φ i /M p ∼ −0.05/β, a sub-Planckian value of the initial scalar value for typical values of β. Early modified gravity solutions require Ω φ,i ∼ 0.05 (EEG), Ω φ,i ∼ 10 −4 (IDEE) around BBN to solve the Hubble problem. These are relatively small differences compared with those needed for canonical oscillating fields Ω quint,i ∼ (T eq /T BBN ) 4 ∼ 10 −24 . While the fine-tuning of initial conditions fares better than in other scenarios, the issues associated to Λ and cosmic acceleration remain.
Among the three mechanisms, EEG is the best candidate to reconcile CMB+BAO and distance ladder, although a combination of EEG and LUPE remains promising in light of those cosmological datasets. These scenarios tension between Planck+BAO and SH0ES to the ≈ 2.5σ level, comparable to other late DE solutions [152]. Analyses involving further observables will further constraint these mechanisms. Particularly, lateuniverse cosmological measurements (e.g. redshift space distortions, weak lensing from galaxy shear and CMB, LSS×CMB cross-correlations or type Ia SNe) will improve the bounds on the coupling β for EEG and Ω φ for LUPE, respectively. A variety of tests can be used to probe these mechanisms further, with precision tests of gravity posing the most outstanding challenge for coupled models.
EEG models need to match the observed strength of gravity measured in the solar system, which is given by the local field value φ loc (t 0 , x 0 ) and its derivative (section V B). The non-linear nature of the problem (Vainshtein screening) and the hierarchy of scales involved (cosmological background to Solar System) require further modeling to reliably address this issue. Two scenarios are plausible: 1) if the local field velocity is comparable to the cosmological value, then EEG is severely limited by the time variation of Newton's constant and the stringent bounds from LLR ( figure 14). 2) If the local field velocity is significantly slower than the cosmological one, a large coupling value is required to recover the correct local strength of gravity today, entering into conflict with cosmology and constraints on scalar forces. A related issue is the interpretation of SNe and other observations in models in which the strength of gravity depends on redshift and host properties.
Big-bang nucleosynthesis is sensitive to the early expansion history, allowing the observed abundance of light elements to place bounds on EEG (section V A). These bounds are by themselves weaker than the Planck+BAO, but when combined might improve limits on the initial effective Planck mass. Note also that the standard BBN relation between the baryon and helium fraction was assumed and played an important role in constraining EEG via the damping tail (figure 11, to right). Varying the helium fraction freely will likely weaken the Planck+BAO limits on EEG.
GW-induced instabilities are sensitive to the lateuniverse evolution (section V C). Avoiding the instability limits the value of the coupling beyond CMB+BAO bounds for EEG and severely limits LUPE, even for Λ = 0. While these limits are enticing, it is important to remember that instabilities signal a breakdown of the theoretical description, rather than a prediction contradicting known data. A UV-complete theory is needed to establish whether EEG & LUPE models can be ruled out by GW-induced instabilities.
EEG, LUPE and IDEE are general mechanisms that can be explored in theories beyond the simple exponentially coupled cubic Galileon. While none of the models studied here is likely to pass all tests, it is plausible that further model building may overcome these difficulties. The notion of IDEE can be generalized beyond the cubic Galileon. 6 In models with a canonical kinetic term one can advance the onset of the kination phase by increasing the hierarchy c 2 /|c 3 |, perhaps even in the prerecombination era. Early modified gravity is also com-6 A straightforward IDEE generalization, known as the nKGB model, is specified by the following Horndeski functions (n = 1 corresponds to the case studied here). A calculation analogous to the one outlined in section (III B) shows the following dependence on the equation of state for the scalar field The conditions for the scalar energy density to dilutes slower than radiation but faster than matter is simply n > 1/2, approaching the matter scaling in the limit n → ∞. A different value of n may improve the behavior of cosmological perturbations relative to the n = 1 case studied here. A generalization of this model has been studied in Ref. [153] as LUPE solution.
patible with quartic or quintic Horndeski terms, as long as a kination phase ensures that the speed of GWs is within acceptable bounds at low redshift. Other simple variations include modifying the coupling function beyond the simple exponential form or adding a potential term. These modifications may help lock up the local value of the scalar field in dense environments, in a manner analogous to the Symmetron model [154]. Needless to say, the properties leading to IDEE, EEG and LUPE (and perhaps completely different solutions to the Hubble problem) are likely to exist in extensions of GR other than Horndeski Gravity. In this sense, this work is only a first systematic exploration of the possibilities of theories beyond Einstein's GR to address the Hubble problem.
The mechanisms described here are extremely predictive. They can be tested using a wide range of observations across vastly different scales and epochs, from precision gravity tests in the laboratory and the Solar System or GW astronomy, all the way to the large-scale structure of the universe, abundance of primordial elements, primary and secondary CMB effects and the cosmic expansion. Future data on these fronts will be able to determine whether the Hubble problem and other cosmological tensions are due to new physics beyond the ΛCDM model. If cosmological tensions endure upcoming scrutiny, combining theoretical and observational insights will be key to illuminate the necessary amendments to the standard model and their fundamental implications for our understanding of nature.
The values of Ω Λ,0 and Ω φ,0 refer to the initial proposal for the sampling distribution, not to hard priors on the parameters. This distinction was necessary because both regions of the parameter space and had to be explored separately. The main results are shown in figure  16, where I have also included the coupled accelerating Λ = 0 model (with Planck 18) to compare the effects of EEG+LUPE, cf. section IV D.
The uncoupled Galileon is all or nothing. Comparison between different initial sampling distributions (a vs b) shows that only one form of energy density dominates. The secondary component is limited to Ω φ,0 < 0.022 in both Λ-dominated scenarios (b) and Ω Λ,0 < 0.019 in φdominated, free m ν (1a) and a more stringent limit Ω Λ,0 < 0.003 for φ-dominated, fixed m ν  [17] for details). Accelerating uncoupled Λ = 0 Galileons are either dominated by the cosmological constant or the LUPE energy density, but the coupled (EEG) model allows a more flexible combination of both dark energy components.
model allows a much wider mixture between the two dark energy components because EEG lifts the restrictions on the acoustic scale. These analyses confirm the role of the neutrino mass in φ-dominated scenarios. Sizeable m ν is necessary both to obtain a better fit and to reconcile the SH0ES value of H 0 , but this only happens in Λ ∼ 0 models. The Λ = 0 (0) and φ-dominated, Λ ∼ 0 (2a) scenarios show only minor differences, consistent with the limits on Λ discussed above. The model with fixed m ν (2a) predicts a Hubble parameter above the SH0ES central value, but is excluded by the analysis with variable neutrino mas (1a), of which (2a) is a particular case. In the Λ-dominated cases (b) the neutrino mass is constrained at a similar level as in standard ΛCDM, with negligible differences between variable (1b) and fixed (2a) m ν .

Appendix C: Coupling and early field dynamics
Let us examine the effects of a non-zero coupling on the initial conditions of the field. Using the general equations presented in section III A I will discuss the Vainshtein screening mechanism and its effect on three possible sources of initial IDEE: pressure-less matter, particles becoming non-relativistic and a hypothetical kination phase. The high efficiency of the Vainshtein mechanism makes those sources completely negligible for all practical purposes. The same suppression of early dynamics makes EEG a much more robust and simple mechanism to lower the acoustic scale. The Galileon is sourced by the trace of the matter and reduced dark energy momentum tensor (see end of section III A) Several early universe phenomena contribute to the source term-sigma Σ, cf. Eq. (29), and may affect the initial kinetic energy of the Galileon. An example is whenever the temperature in the early universe drops below the mass of a particle: for some time that particle remains important in the energy budget, while becoming partially non-relativistic and thus contributing to Σ. This phenomenon is known as "kicks" in the context of chameleon theories [124,155] (see Ref. for a study of theories with non-canonical kinetic terms of the Dirac-Born-Infeld type [156]). Phase transitions contribute similarly to the shift-charge density. At the end of this section I will also examine the effects of a hypothetical kination phase, the most favorable situation to overcome the cosmological Vainshtein screening.
We will express the contribution of a kick to the integral in Eq (26) as J ∝ daH(a)a 2 Σ(T (a)) ∼ H 0 Ω R a eΣ , assuming radiation domination and neglecting the effect of Σ on the expansion (28). The above approximating is equivalent to treating the kick as as a step function with a start a e , which is adequate to give an idea of the order of magnitude and time dependence. Typical contributions for massive standard model particles are Σ ∼ 0.05 − 0.1 (See Ref. [124] a detailed computation).
It is possible to express the shift-charge as an energy density fraction for the Galileon using Eqs. (19,15). If the cubic term dominates then where the last equality uses the simplified kick expression (C1). While this contribution dilutes slower than radiation (Ω φ,3 ∝ √ a), the initial kick is suppressed by a 2 e 1. This dependence implies that kicks at an earlier epoch are less important, making it very hard to invoke early universe physics (e.g. heavy BSM particles with m > m τ ).
The scaling of the cubic Galileon reflects the cosmological Vainshtein screening. This is very different in the case of a quadratic kinetic term, for whicĥ For a canonical kinetic term, a kick contributes a sizeable amount of kinetic energy in the fieldΩ φ,2 ∼Σβ 2 , which nonetheless kinates away rapidly asΩ φ,2 ∼ a −2 .
In contrast, the cubic Galileon is very hard to excite, but any energy injected into the field is persistent, withΩ φ,3 growing in the radiation era as characteristic of IDEE models, cf section III B. Non-luminal Galileons scale more favorably with cosmic expansion, but are equally hard to excite due to the Vainshtein mechanism. If the quartic or quintic term were to dominate the evolution (both the shift-charge and the energy density), the contribution of a kick readŝ Ω φ,4 ∼ so early kicks are suppressed by powers of the initial scale factor. Note that whileΩ φ,4−5 grows faster than in the cubic case, this dependence does not compensate for the Vainshtein mechanism, seen here as positive powers of a e , which make kicks at very early times negligible. Note also that the enhancement produced by the small coefficients c 4 , c 5 (less screening) will not lead to a large kick, but rather to the cubic or canonical term becoming the relevant one. Just for fun, let us now examine best-case scenario to generate a large IDEE fraction through a coupling. The best case to generate a large shift-charge density would be a kination phase (e.g. driven by the inflaton) with w m ≈ −1, Σ ≈ −2, H ∝ a −3 . Note that negative Σ requires a negative coupling β < 0 to produce a positive shift-charge. Then the integral in Eq. (26) reads daa 2 H(a)Σ = −2H 0 Ω r a e ln (a e /a i ) , where I assumed that kination dominates from a i to a e and the universe becomes radiation dominated at a e (hence relating H e = H 0 √ Ω r a 2 e ). To evaluate the impact of a kination phase on the Galileon density fraction one can substituteΣ → −2 ln (a e /a i ) in the expressions in the previous section. The logarithmic factor gives a mild dependence on the duration of the kination phase, which can be made arbitrarily large in the limit a i → 0, if the kination phase last long enough.
While possible, imparting a substantial initial energy to the Galileon using a kination phase is extremely un-realistic. The problem is the very rapid scaling of the energy density during a kination phase, with ρ i /ρ e = (a e /a i ) 6 . The most favorable scenario to affect the acoustic scale via IDEE requires kination to end right before nucleosynthesis, a e ∼ 10 10 , while at the same time producingΩ φ (z BBN ) ∼ 10 −4 . This would require the kick to be as large asΣβ = −2β log(a e /a i ) ∼ 2 · 10 9 (Ω φ,i /10 −4 ) 2/3 (10 −10 /a e ), corresponding to an initial energy density at the beginning of kination given by ρ i /ρ e = e |Σ/β| = e |β| −1 2·10 9 ∼ 6 · 10 868588963 , where the last value assumes β ∼ 1. 7 Needless to say, this energy scale is deeply trans-planckian, well beyond the range of validity of the theory as well as the range of validity of classical gravity.
It is clear from the above discussion that the cosmological Vainshtein screening precludes any early universe process to produce a sizeable contribution to IDEE. Inflation would dilute the initial energy density of the scalar field very efficiently, requiring a mechanism to produce a sizeable amount of IDEE at reheating or later. This necessarily involves physics beyond the classical coupled Galileon theory, possibly through an ultraviolet completion. This may happen in scenarios of Galilean genesis [157], a variant of the coupled cubic Galileon in which the scalar field is responsible for setting the initial conditions in the early universe. In this scenario, reheating is conjectured to occur when the field configuration exits the effective field theory regime of validity. While a high-energy completion of the theory is necessary for a first principle calculation, it is plausible that the Galileon field producing IDEE might be generated with a sizeable kinetic energy (note that this scalar field might be different from the one causing Galilean genesis).
The Vainshtein mechanism ensures that the initial effective Planck mass M 2 * (φ) is robust against physical processes in the early universe. The smallness of the relative variation of the fieldφ/(Hφ) guarantees that M 2 * (φ) will remain approximately constant until the Hubble rate decreases to a value H ∼ H 0 / |c 3 |. Thus whatever the initial condition φ i set in the early universe, its effect on the strength of gravity is robust by virtue of the same physics that prevent the generation of IDEEΩ φ . It is interesting that, already at the theoretical level, Enhanced Early Gravity is much more robust.