Cubic Galileon Theory: Dynamics & Applications

• Cubic Galileon theory is a scalar–tensor model featuring a unique cubic derivative self-coupling that maintains second-order field equations and avoids Ostrogradski ghosts.
• The theory alters cosmological dynamics by modifying Friedmann equations, enabling late-time acceleration, early universe modifications, and effective Vainshtein screening for structure formation.
• Its rich phenomenology spans linear perturbations, black hole solutions with scalar hair, and quantum corrections, highlighting its role in both astrophysical and modified-gravity scenarios.

Cubic Galileon theory denotes the class of scalar or scalar–tensor models whose defining interaction is the lowest nontrivial Galileon derivative self-coupling, conventionally written as $(\partial\phi)^2\Box\phi$ or $X\Box\phi$, added to Einstein gravity and often supplemented by a potential, matter couplings, or further modified-gravity structure. In flat space the cubic term is Galileon-invariant up to a total derivative and still yields second-order equations of motion, thereby avoiding Ostrogradski ghosts at the classical level; in curved space it appears naturally as a Horndeski sector and in a range of softly broken variants used for late-time acceleration, early modified gravity, Vainshtein screening, and compact-object solutions (Saltas et al., 2016, Babichev et al., 2016, Ye et al., 2024).

1. Covariant formulation and defining structure

The minimal flat-space cubic Galileon contains the quadratic and cubic operators

$${\cal L}_2=-\alpha_2\,\partial_\mu\phi\,\partial^\mu\phi,\qquad {\cal L}_3=-\alpha_3\,\partial_\mu\phi\,\partial^\mu\phi\,\Box\phi,$$

with boundary terms required for a well-posed variational principle under Dirichlet conditions (Sivanesan, 2011). In the covariant language used in Horndeski theory, a representative shift-symmetric action is

$$S=\int d^Dx\,\sqrt{-g}\,\Bigl[\zeta(R-2\Lambda)-\eta\,g^{\mu\nu}\partial_\mu\phi\partial_\nu\phi+\gamma\,(\Box\phi)\,g^{\mu\nu}\partial_\mu\phi\partial_\nu\phi\Bigr],$$

which corresponds to

$$G_2(X)=2\eta X-2\zeta\Lambda,\qquad G_3(X)=-2\gamma X,\qquad G_4=\zeta,\qquad G_5=0$$

with $$X=-\tfrac12 g^{\mu\nu}\partial_\mu\phi\partial_\nu\phi$$ (Babichev et al., 2016).

Cosmological realizations often adopt the Einstein–Hilbert term plus a canonical kinetic sector, a potential, and the cubic operator. Ye and Silvestri study

$$S=\int d^4x\,\sqrt{-g}\,\Bigl[\frac{M_p^2}{2}R+X-V(\phi)-\xi\,X\Box\phi\Bigr],$$

with quartic potential $V(\phi)=V_0\phi^4+V_\Lambda$; in their construction the model is a Horndeski/Galileon theory with luminal tensor speed and a single cubic-Galileon operator (Ye et al., 2024). Related late-time dark-energy models retain the same derivative structure but use a linear potential $V(\phi)=\mu^3\phi$ or $V(\phi)=V_0\phi$, so that for $X\Box\phi$0 they reduce to quintessence (Dinda, 2018, Zhang et al., 2020).

The cubic operator also appears in nonminimally coupled or matter-coupled variants. In the Einstein-frame model used to constrain the maximum turn-around radius, matter couples through $X\Box\phi$1, equivalently through the physical metric $X\Box\phi$2 to first order in $X\Box\phi$3 (Bhattacharya et al., 2015). In a disformal realization, Standard Model matter couples to

$X\Box\phi$4

while the Galileon sector contains both the quadratic and cubic operators (Lawrence et al., 2020). A further extension embeds a cubic Galileon scalar into the dRGT massive-gravity potential, where the quasi-dilaton field $X\Box\phi$5 carries a cubic kinetic coupling $X\Box\phi$6 (Aslmarand et al., 2021).

A recurrent structural point is that the cubic derivative interaction modifies dynamics without introducing higher-order field equations. This property underlies its use in both cosmology and strong-field gravity, but the concrete phenomenology depends strongly on the chosen potential, matter coupling, and asymptotic branch. This suggests that “the cubic Galileon” is less a single model than a restricted operator content realized in several inequivalent theories.

2. Cosmological background dynamics

In a spatially flat FRW metric,

$X\Box\phi$7

the cubic operator alters both the Friedmann equations and the scalar equation. For the $X\Box\phi$8EDE model one obtains

$X\Box\phi$9

$${\cal L}_2=-\alpha_2\,\partial_\mu\phi\,\partial^\mu\phi,\qquad {\cal L}_3=-\alpha_3\,\partial_\mu\phi\,\partial^\mu\phi\,\Box\phi,$$0

while the scalar equation on FRW becomes

$${\cal L}_2=-\alpha_2\,\partial_\mu\phi\,\partial^\mu\phi,\qquad {\cal L}_3=-\alpha_3\,\partial_\mu\phi\,\partial^\mu\phi\,\Box\phi,$$1

In the slow-roll or frozen tracking regime,

$${\cal L}_2=-\alpha_2\,\partial_\mu\phi\,\partial^\mu\phi,\qquad {\cal L}_3=-\alpha_3\,\partial_\mu\phi\,\partial^\mu\phi\,\Box\phi,$$2

(Ye et al., 2024).

Late-time dark-energy constructions based on a linear potential produce modified background densities and pressures of the form

$${\cal L}_2=-\alpha_2\,\partial_\mu\phi\,\partial^\mu\phi,\qquad {\cal L}_3=-\alpha_3\,\partial_\mu\phi\,\partial^\mu\phi\,\Box\phi,$$3

and the scalar equation

$${\cal L}_2=-\alpha_2\,\partial_\mu\phi\,\partial^\mu\phi,\qquad {\cal L}_3=-\alpha_3\,\partial_\mu\phi\,\partial^\mu\phi\,\Box\phi,$$4

(Dinda, 2018). In that class, the $${\cal L}_2=-\alpha_2\,\partial_\mu\phi\,\partial^\mu\phi,\qquad {\cal L}_3=-\alpha_3\,\partial_\mu\phi\,\partial^\mu\phi\,\Box\phi,$$5 term alone cannot drive stable acceleration, hence the inclusion of $${\cal L}_2=-\alpha_2\,\partial_\mu\phi\,\partial^\mu\phi,\qquad {\cal L}_3=-\alpha_3\,\partial_\mu\phi\,\partial^\mu\phi\,\Box\phi,$$6 (Dinda, 2018).

The background dynamics are branch-dependent. In the vacuum cubic-derivative system analyzed by Leon and collaborators, a local late-time attractor representing phantom behavior appears and is associated with a Big Rip singularity, whereas self-accelerating solutions do not arise in that model (Arcia et al., 2015). By contrast, in the linearly and quadratically coupled cubic Galileon models with linear potentials, late-time acceleration occurs without a transition from accelerating phase to phantom phase in the future; because of the matter coupling, one distinguishes a native equation of state from an effective equation of state (Kim et al., 2013). In the integrable model with

$${\cal L}_2=-\alpha_2\,\partial_\mu\phi\,\partial^\mu\phi,\qquad {\cal L}_3=-\alpha_3\,\partial_\mu\phi\,\partial^\mu\phi\,\Box\phi,$$7

the minisuperspace system admits a second conservation law and therefore an integrable dynamical system with additional exact critical points (Giacomini et al., 2017).

A separate cosmological branch is the tracker cubic Galileon used in Planck-era analyses, defined by

$${\cal L}_2=-\alpha_2\,\partial_\mu\phi\,\partial^\mu\phi,\qquad {\cal L}_3=-\alpha_3\,\partial_\mu\phi\,\partial^\mu\phi\,\Box\phi,$$8

In that setting the background is fully fixed by $${\cal L}_2=-\alpha_2\,\partial_\mu\phi\,\partial^\mu\phi,\qquad {\cal L}_3=-\alpha_3\,\partial_\mu\phi\,\partial^\mu\phi\,\Box\phi,$$9 once $$S=\int d^Dx\,\sqrt{-g}\,\Bigl[\zeta(R-2\Lambda)-\eta\,g^{\mu\nu}\partial_\mu\phi\partial_\nu\phi+\gamma\,(\Box\phi)\,g^{\mu\nu}\partial_\mu\phi\partial_\nu\phi\Bigr],$$0 is chosen, yielding

$$S=\int d^Dx\,\sqrt{-g}\,\Bigl[\zeta(R-2\Lambda)-\eta\,g^{\mu\nu}\partial_\mu\phi\partial_\nu\phi+\gamma\,(\Box\phi)\,g^{\mu\nu}\partial_\mu\phi\partial_\nu\phi\Bigr],$$1

(Barreira et al., 2014). The coexistence of tracker, potential-driven, coupled, and early-modified-gravity realizations is central to the literature: the same cubic operator can support qualitatively different background histories, including quasi-$$S=\int d^Dx\,\sqrt{-g}\,\Bigl[\zeta(R-2\Lambda)-\eta\,g^{\mu\nu}\partial_\mu\phi\partial_\nu\phi+\gamma\,(\Box\phi)\,g^{\mu\nu}\partial_\mu\phi\partial_\nu\phi\Bigr],$$2CDM behavior, de Sitter-like acceleration, early dark-energy episodes near recombination, or phantom late-time vacuum asymptotics.

3. Linear perturbations, braiding, screening, and structure formation

In the Horndeski parameterization of Bellini and Sawicki, the $$S=\int d^Dx\,\sqrt{-g}\,\Bigl[\zeta(R-2\Lambda)-\eta\,g^{\mu\nu}\partial_\mu\phi\partial_\nu\phi+\gamma\,(\Box\phi)\,g^{\mu\nu}\partial_\mu\phi\partial_\nu\phi\Bigr],$$3EDE model has only two nonzero $$S=\int d^Dx\,\sqrt{-g}\,\Bigl[\zeta(R-2\Lambda)-\eta\,g^{\mu\nu}\partial_\mu\phi\partial_\nu\phi+\gamma\,(\Box\phi)\,g^{\mu\nu}\partial_\mu\phi\partial_\nu\phi\Bigr],$$4-functions,

$$S=\int d^Dx\,\sqrt{-g}\,\Bigl[\zeta(R-2\Lambda)-\eta\,g^{\mu\nu}\partial_\mu\phi\partial_\nu\phi+\gamma\,(\Box\phi)\,g^{\mu\nu}\partial_\mu\phi\partial_\nu\phi\Bigr],$$5

so the characteristic effect is kinetic braiding (Ye et al., 2024). On subhorizon scales the effective Newton constant for growth is

$$S=\int d^Dx\,\sqrt{-g}\,\Bigl[\zeta(R-2\Lambda)-\eta\,g^{\mu\nu}\partial_\mu\phi\partial_\nu\phi+\gamma\,(\Box\phi)\,g^{\mu\nu}\partial_\mu\phi\partial_\nu\phi\Bigr],$$6

while the scalar sound speed satisfies

$$S=\int d^Dx\,\sqrt{-g}\,\Bigl[\zeta(R-2\Lambda)-\eta\,g^{\mu\nu}\partial_\mu\phi\partial_\nu\phi+\gamma\,(\Box\phi)\,g^{\mu\nu}\partial_\mu\phi\partial_\nu\phi\Bigr],$$7

In the early tracking phase, no-ghost and no-gradient stability require $$S=\int d^Dx\,\sqrt{-g}\,\Bigl[\zeta(R-2\Lambda)-\eta\,g^{\mu\nu}\partial_\mu\phi\partial_\nu\phi+\gamma\,(\Box\phi)\,g^{\mu\nu}\partial_\mu\phi\partial_\nu\phi\Bigr],$$8 (Ye et al., 2024). The background effect of braiding can be summarized as

$$S=\int d^Dx\,\sqrt{-g}\,\Bigl[\zeta(R-2\Lambda)-\eta\,g^{\mu\nu}\partial_\mu\phi\partial_\nu\phi+\gamma\,(\Box\phi)\,g^{\mu\nu}\partial_\mu\phi\partial_\nu\phi\Bigr],$$9

which boosts $$G_2(X)=2\eta X-2\zeta\Lambda,\qquad G_3(X)=-2\gamma X,\qquad G_4=\zeta,\qquad G_5=0$$0 when $$G_2(X)=2\eta X-2\zeta\Lambda,\qquad G_3(X)=-2\gamma X,\qquad G_4=\zeta,\qquad G_5=0$$1, while in the Poisson equation it adds an $$G_2(X)=2\eta X-2\zeta\Lambda,\qquad G_3(X)=-2\gamma X,\qquad G_4=\zeta,\qquad G_5=0$$2 attractive fifth force for subhorizon modes (Ye et al., 2024).

The cubic interaction is also the canonical source of Vainshtein screening. In the nonrelativistic, quasi-static limit, the fifth force on a test mass is

$$G_2(X)=2\eta X-2\zeta\Lambda,\qquad G_3(X)=-2\gamma X,\qquad G_4=\zeta,\qquad G_5=0$$3

and around a spherical source of mass $$G_2(X)=2\eta X-2\zeta\Lambda,\qquad G_3(X)=-2\gamma X,\qquad G_4=\zeta,\qquad G_5=0$$4 the Vainshtein radius is

$$G_2(X)=2\eta X-2\zeta\Lambda,\qquad G_3(X)=-2\gamma X,\qquad G_4=\zeta,\qquad G_5=0$$5

For $$G_2(X)=2\eta X-2\zeta\Lambda,\qquad G_3(X)=-2\gamma X,\qquad G_4=\zeta,\qquad G_5=0$$6 the extra force is strongly screened, while for $$G_2(X)=2\eta X-2\zeta\Lambda,\qquad G_3(X)=-2\gamma X,\qquad G_4=\zeta,\qquad G_5=0$$7 it approaches the unscreened ratio $$G_2(X)=2\eta X-2\zeta\Lambda,\qquad G_3(X)=-2\gamma X,\qquad G_4=\zeta,\qquad G_5=0$$8 (Chan et al., 2018). In the turn-around analysis of bound cosmic structures, the same screening scale must exceed the maximum turn-around radius in order for the solution to lie in the stable screened region (Bhattacharya et al., 2015).

Weak-lensing and LSS signatures have been computed in several model realizations. In the linear-potential cubic Galileon studied by Dinda, the growth factor is suppressed relative to $$G_2(X)=2\eta X-2\zeta\Lambda,\qquad G_3(X)=-2\gamma X,\qquad G_4=\zeta,\qquad G_5=0$$9CDM by up to $$X=-\tfrac12 g^{\mu\nu}\partial_\mu\phi\partial_\nu\phi$$0 today, the matter power spectrum at $$X=-\tfrac12 g^{\mu\nu}\partial_\mu\phi\partial_\nu\phi$$1 is lower by $$X=-\tfrac12 g^{\mu\nu}\partial_\mu\phi\partial_\nu\phi$$2 at $$X=-\tfrac12 g^{\mu\nu}\partial_\mu\phi\partial_\nu\phi$$3 and by $$X=-\tfrac12 g^{\mu\nu}\partial_\mu\phi\partial_\nu\phi$$4 at $$X=-\tfrac12 g^{\mu\nu}\partial_\mu\phi\partial_\nu\phi$$5, and the convergence power spectrum differs from $$X=-\tfrac12 g^{\mu\nu}\partial_\mu\phi\partial_\nu\phi$$6CDM by $$X=-\tfrac12 g^{\mu\nu}\partial_\mu\phi\partial_\nu\phi$$7 at $$X=-\tfrac12 g^{\mu\nu}\partial_\mu\phi\partial_\nu\phi$$8 and $$X=-\tfrac12 g^{\mu\nu}\partial_\mu\phi\partial_\nu\phi$$9 at $$S=\int d^4x\,\sqrt{-g}\,\Bigl[\frac{M_p^2}{2}R+X-V(\phi)-\xi\,X\Box\phi\Bigr],$$0 (Dinda, 2018). In the N-body study that approximates the modified force law by injecting $$S=\int d^4x\,\sqrt{-g}\,\Bigl[\frac{M_p^2}{2}R+X-V(\phi)-\xi\,X\Box\phi\Bigr],$$1 into the tree-PM solver, the $$S=\int d^4x\,\sqrt{-g}\,\Bigl[\frac{M_p^2}{2}R+X-V(\phi)-\xi\,X\Box\phi\Bigr],$$2 matter power spectrum is suppressed by about $$S=\int d^4x\,\sqrt{-g}\,\Bigl[\frac{M_p^2}{2}R+X-V(\phi)-\xi\,X\Box\phi\Bigr],$$3 for $$S=\int d^4x\,\sqrt{-g}\,\Bigl[\frac{M_p^2}{2}R+X-V(\phi)-\xi\,X\Box\phi\Bigr],$$4, the marked density power spectrum differs by about $$S=\int d^4x\,\sqrt{-g}\,\Bigl[\frac{M_p^2}{2}R+X-V(\phi)-\xi\,X\Box\phi\Bigr],$$5 around $$S=\int d^4x\,\sqrt{-g}\,\Bigl[\frac{M_p^2}{2}R+X-V(\phi)-\xi\,X\Box\phi\Bigr],$$6, halo abundances above $$S=\int d^4x\,\sqrt{-g}\,\Bigl[\frac{M_p^2}{2}R+X-V(\phi)-\xi\,X\Box\phi\Bigr],$$7 are lower by $$S=\int d^4x\,\sqrt{-g}\,\Bigl[\frac{M_p^2}{2}R+X-V(\phi)-\xi\,X\Box\phi\Bigr],$$8, and the number of low-density cells is higher than in $$S=\int d^4x\,\sqrt{-g}\,\Bigl[\frac{M_p^2}{2}R+X-V(\phi)-\xi\,X\Box\phi\Bigr],$$9CDM (Zhang et al., 2020). The same study emphasizes that the full nonlinear Galileon field is not solved on the mesh, so a dynamically formed screened Vainshtein interior is absent (Zhang et al., 2020).

4. Compact objects and strong-gravity configurations

Shift-symmetric cubic Galileon theory admits black-hole solutions with linearly time-dependent scalar hair,

$V(\phi)=V_0\phi^4+V_\Lambda$0

which provides a controlled evasion of static no-hair theorems (Babichev et al., 2016). In the spherically symmetric sector, analytic three-dimensional BTZ-like solutions exist, and in four dimensions one finds numerical black-hole families parameterized by the scalar “velocity” $V(\phi)=V_0\phi^4+V_\Lambda$1. Their asymptotic behavior matches the same homogeneous de Sitter background for fixed branch, while near the horizon the metric functions deviate from Schwarzschild–(A)dS (Babichev et al., 2016). Babichev and collaborators identify three branches of effective cosmological constant in the de Sitter case and interpret $V(\phi)=V_0\phi^4+V_\Lambda$2 as primary hair (Babichev et al., 2016).

Rotating solutions are more unusual. Van Aelst and collaborators constructed asymptotically flat rotating black holes with ansatz

$V(\phi)=V_0\phi^4+V_\Lambda$3

and found that the lapse behaves as

$V(\phi)=V_0\phi^4+V_\Lambda$4

so the usual $V(\phi)=V_0\phi^4+V_\Lambda$5 mass monopole is absent and the Komar mass vanishes (Aelst et al., 2019). These rotating solutions show significant deviations from Kerr: for fixed horizon angular velocity they carry larger angular momentum, admit a “super-Kerr” regime with $V(\phi)=V_0\phi^4+V_\Lambda$6, have smaller surface gravity, and develop thinner ergoregions (Aelst et al., 2019).

A later fully consistent 3+1 spectral construction confirms the vanishing ADM mass numerically while preserving nonzero angular momentum (Grandclément, 2023). In this formulation the black holes are treated as apparent horizons in equilibrium, using maximal slicing and spatial harmonic gauge. The same computation shows that the surface gravity is not constant on the horizon, i.e. $V(\phi)=V_0\phi^4+V_\Lambda$7, and relates this to violation of the weak and dominant energy conditions by the Galileon energy density (Grandclément, 2023). A common misconception is that vanishing Komar or ADM mass would imply a trivial geometry; the numerical black-hole solutions explicitly contradict that implication, because the metric is nonflat and the scalar stores energy in a nonlinear screened configuration (Aelst et al., 2019, Grandclément, 2023).

5. Hamiltonian structure, effective metrics, and quantum corrections

At cubic order, the Hamiltonian of the flat-space single-Galileon theory contains the canonical kinetic term, the spatial gradient term, and two cubic interactions,

$V(\phi)=V_0\phi^4+V_\Lambda$8

For a static point source, the radial equation admits two branches $V(\phi)=V_0\phi^4+V_\Lambda$9 with

$V(\phi)=\mu^3\phi$0

and the corresponding energies have equal magnitude and opposite sign,

$V(\phi)=\mu^3\phi$1

The negative-energy branch is interpreted as a nonlinear manifestation of ghost instability, while the short-distance divergence of the point source is regularized by the dominant cubic term (Sivanesan, 2011).

Perturbations around nontrivial backgrounds propagate in an emergent effective metric rather than the spacetime metric. For steady-state spherical accretion onto Schwarzschild, perturbing $V(\phi)=\mu^3\phi$2 yields a massless wave equation

$V(\phi)=\mu^3\phi$3

with

$V(\phi)=\mu^3\phi$4

The sonic horizon is defined by $V(\phi)=\mu^3\phi$5, and in the Schwarzschild case regularity fixes $V(\phi)=\mu^3\phi$6 and selects $V(\phi)=\mu^3\phi$7. Using the Moncrief energy method, one finds $V(\phi)=\mu^3\phi$8 and $V(\phi)=\mu^3\phi$9, establishing linear stability of the accreting configuration (Bergliaffa et al., 2016). This effective-metric viewpoint is essential: causal propagation and stability are not determined solely by the background spacetime metric.

At the quantum level, the Vilkovisky–DeWitt covariant effective action reveals loop-generated higher-derivative operators even when the classical theory is cubic. The divergent part of the one-loop effective action includes

$V(\phi)=V_0\phi$0

with coefficients

$V(\phi)=V_0\phi$1

up to the common dimensional-regularization pole (Saltas et al., 2016). Their relevance is controlled by the ratios $V(\phi)=V_0\phi$2, $V(\phi)=V_0\phi$3, and $V(\phi)=V_0\phi$4. A plausible implication is that the classical second-order structure is an effective low-energy property rather than a radiatively closed one; near the cutoff, extra counterterms are unavoidable.

6. Observational status, constraints, and unresolved issues

Observational conclusions depend sharply on which cubic-Galileon realization is tested. In the early-modified-gravity model $V(\phi)=V_0\phi$5EDE, Planck PR4 + BAO + Pantheon+ plus Cepheid-calibrated Pantheon+ type Ia supernovae favor a nonzero cubic contribution,

$V(\phi)=V_0\phi$6

exclude zero at about $V(\phi)=V_0\phi$7, shift

$V(\phi)=V_0\phi$8

and yield

$V(\phi)=V_0\phi$9

with Bayes factor $X\Box\phi$00 over canonical EDE despite one extra parameter (Ye et al., 2024). The same work emphasizes that the model does not resolve the $X\Box\phi$01 tension because the modified-gravity fifth force is always attractive and very short-lived (Ye et al., 2024).

The post-Planck analysis of the original tracker Cubic Galileon branch reaches a different conclusion. With Planck CMB, lensing, and BAO, the massless-neutrino version fits poorly, while the version with free neutrino mass can fit the data with best-fit $X\Box\phi$02, $X\Box\phi$03, and $X\Box\phi$04 (Barreira et al., 2014). However, the cubic branch predicts that the lensing potential deepens at late times on subhorizon scales, implying a negative ISW effect on those scales, which is at odds with the observational suggestion of a positive ISW effect (Barreira et al., 2014). This ISW issue has become one of the standard phenomenological objections to the simplest late-time self-accelerating cubic Galileon.

Direct background probes constrain another corner of parameter space. Using cosmic-chronometer $X\Box\phi$05 data with $X\Box\phi$06 and $X\Box\phi$07 priors, Bellini and Jimenez found that in the simplest cubic model with linear potential the acceleration is effectively provided by the constant part of the potential, with only about $X\Box\phi$08 of the acceleration energy density attributable to genuine Galileon dynamics; $X\Box\phi$09 stays within $X\Box\phi$10 at $X\Box\phi$11, and the Bayesian evidence gives odds $X\Box\phi$12 against the cubic Galileon relative to $X\Box\phi$13CDM (Bellini et al., 2013).

Astrophysical-scale constraints are likewise model-specific. The maximum turn-around-radius analysis gives

$X\Box\phi$14

improving earlier solar-system bounds (Bhattacharya et al., 2015). In the rotation-curve application, the Milky Way and a 28-galaxy SPARC subsample are fitted without dark matter using $X\Box\phi$15 in the approximate range $X\Box\phi$16, with $X\Box\phi$17 for the Milky Way and typical reduced $X\Box\phi$18 (Chan et al., 2018). By contrast, in the disformally coupled scenario the cubic Galileon is assumed to be subdominant precisely to evade the ISW exclusion; there the relevant bounds are collider limits

$X\Box\phi$19

from ATLAS $X\Box\phi$20 TeV data, while GW170817 does not constrain the model because $X\Box\phi$21 under the allowed parameter choices (Lawrence et al., 2020).

Cubic Galileon structure also survives in modified-gravity extensions beyond standard Horndeski dark energy. In cubic Galileon massive gravity, self-accelerating backgrounds exist, the late-time Hubble law is parameterized by

$X\Box\phi$22

and a fit to the 557-object Union2 SNIa sample gives

$X\Box\phi$23

while tensor modes obey $X\Box\phi$24 (Aslmarand et al., 2021).

The present status is therefore mixed rather than uniform. Some realizations are strongly constrained or effectively reduced to $X\Box\phi$25CDM at the background level, some remain viable only in extended parameter spaces, and some—such as the cubic-braided early-modified-gravity model—gain support by addressing the Hubble tension (Ye et al., 2024, Bellini et al., 2013, Barreira et al., 2014). Proposed next steps include extending to more general $X\Box\phi$26 and $X\Box\phi$27 couplings, imposing positivity-bound priors for UV consistency, exploring stochastic GW-background and B-mode signatures of oscillating Galileons, and releasing fully covariant Horndeski modules for Einstein–Boltzmann evolution (Ye et al., 2024).