The $\mathrm{\Lambda }$-cold-dark-matter ($\mathrm{\Lambda }\mathrm{CDM}$) model together with the flat Friedmann-Lemaître-Robertson-Walker (FLRW) background.

For clarity, the references cited under each keyword are limited to those addressing a specific type of solution to the Hubble tension. Our aim here is not to provide a complete literature review on the origins of the controversy, but rather to offer illustrative examples. Such a comprehensive review is urgently needed for the community, but appears to be a formidable challenge to accomplish.

While ordinary matter can exert an indirect influence on scalar field dynamics through Hubble friction, this does not provide a mechanism to control the emergence time of EDE.

Zumalacárregui [285] constructed a trigger-free solution, in which the key ingredient is an energy density that dilutes faster than matter but slower than radiation. This naturally forces the peak of its relative energy density to coincide with the epoch of matter-radiation equality. However, the slower dilution rate compared to radiation (after the peak) appears disfavored by cosmological constraints. Furthermore, this scenario requires a fine-tuning of the initial energy density of the EDE component in the earlier universe to achieve ${\mathrm{\Omega }}_{\mathrm{EDE}}\approx 10\%$ near equality.

Carrillo González et al. [286] dispute these arguments. However, they do not provide a concrete and viable parameter setting (in their Sec. 2 C) to directly address the issue, particularly the $\mathrm{\Delta }$ criterion proposed in Ref. [68]. Separately, Kamionkowski & Mathur [287] proposed a new neutrino-triggered EDE model, in which the dependence on initial conditions leaves the coincidence problem fundamentally open.

As stated by Lin et al. [63], the model proposed by Karwal et al. [62] does not correctly resolve the EDE coincidence problem. Here we refer to the model proposed by Lin et al. [63].

See discussions above Eq. (8) in Ref. [66]. While introducing a cosmological constant can provide a temporary fix, it makes the model theoretically unappealing.

Cosmic evolution spans a wide range of energy scales. A theory that solves every observational anomaly by inserting new physics at a specific energy scale becomes increasingly cumbersome. Furthermore, a collection of phenomenological patches may conceal a more fundamental, unified principle. One should not introduce too many anthropic energy scales into a gravitational system. It could be argued that this statement is challenged by the case of dark energy. Nevertheless, despite its limited acceptance, there are dark energy models that avoid any ad hoc energy scale [288, 289].

According to the KiDS-Legacy result [290], the $S_8$ tension disappears. The DES Y3 Blue Shear analysis also suggests that source-galaxy selection and intrinsic alignments may affect the inferred $S_8$ tension [291]. The existence of this tension and the possible systematic errors in the observations remain an ongoing debate [28].

A realistic theory of gravity should be compatible with observed black-hole phenomenology, as supported by gravitational-wave detections [292] and Event Horizon Telescope observations [293]. Nevertheless, the absence of identified black-hole solutions beyond Schwarzschild and the scarcity of asymptotically flat branches in Sec. 5 should not simply be taken as negative results. The key question is not only which static exact solutions exist, but which local gravitational configurations are dynamically realized in the real Universe. In particular, the final exterior state may still be the Schwarzschild branch with $S_{\mu }=0$. Gravitational collapse and stability analyses are needed to answer this question.

Unfortunately, both Starobinsky [6] and Sato [96] failed to recognize that the exponential expansion could resolve the horizon and flatness problems.

The saddle-type indicates a critical point that attracts trajectories in one direction and repels in the other. The attractor in this context specifically refers to the stable manifold—a special trajectory along the repulsive direction that attracts nearby orbits. Figure 1(b) in Ref. [114] provides an illustration.

https://www.xact.es.

Conformal matter coupling can preserve the usual matter conservation law in the corresponding physical metric [294]. However, it operates through a metric rescaling rather than direct couplings between a new field and fluid thermodynamic variables, which limits its applicability to our EDE purpose. A recent Herglotz-type extension of $f(R,\text{Matter})$ gravity may offer another route by introducing an additional closed one-form whose contribution cancels the non-conservation current [295]. This is distinct from our direct-construction approach.

For simplicity, in Sec. 2, we can assume that the metric is dimensionless and $[\mathrm{d}{x}^{\mu }]=\text{length}$. This assumption breaks down for the FLRW metric, which would assign a different dimension to $S_0$ in Sec. 3. Nevertheless, this inconsistency does not affect the discussion in Sec. II.

More specifically, the ${\nabla }_{\mu }{S}_{\nu }$, $c_2$, $c_4$ and $c_6$ terms in Eq. (2.1) are unable to ensure the stability of the critical point $F_1$ in Fig. 2.

Generally, a rescaled $G$ is required in modified gravity [231]. However, Sec. 4 shows that it is not needed for GCT.

It is unnecessary to introduce a free coefficient for the ${\nabla }_{\mu }{S}_{\nu }+{\nabla }_{\nu }{S}_{\mu }$ term in $S_{\mu \nu }$ due to a scaling property—if $S_{\mu \nu }=c_0({\nabla }_{\mu }{S}_{\nu }+{\nabla }_{\nu }{S}_{\mu })+\cdots$, then the transformation $\{{\tilde{S}}_{\mu }=c_0{S}_{\mu },{\tilde{c}}_1=c_1,{\tilde{c}}_2=c_2/c_0,{\tilde{c}}_3=c_3/c_0,{\tilde{c}}_4=c_4/c_0^2,{\tilde{c}}_5=c_5/c_0,{\tilde{c}}_6=c_6/c_0^2,{\tilde{c}}_7=c_7,{\tilde{c}}_8=c_8\}$ absorbs the coefficient and leaves the effective gravitational tensor of the vector field unchanged.

One may recall the linearized Einstein field equations.

Here we do not consider the case of zero eigenvalues.

This follows directly from the definition of $x_3$.

Our initial motivation for introducing the Poincaré sphere was to clarify the behavior of trajectories in the upper-right corner of Fig. 2 Left, which ultimately leads the analysis to the critical points at infinity.

Roughly speaking, $x_3=2$ is an asymptote of the red and green trajectories in Fig. 2. Extending the range of the phase portrait makes the trend clear [167].

The continuous infinity arc can serve as an attractor or a source in the phase portrait (see Example 1 on p. 20 in Ref. [119] and its Poincaré map for an illustration).

The results or plots that support this statement are not shown in this paper.

Based on the definition of $x_4$, the case where the equality holds can be understood from a limiting perspective.

As an illustration, in the $z_i$ panel of Fig. 6, the gray dashed curves represent trajectories with $z_0<0$, while the arrows still indicate the direction of increasing $\sigma$.

For the zero-eigenvalue case, we simply assume that the stability could be determined by the quadratic terms, without proof from center manifold theory.

Our analysis concerns only two-dimensional real systems, so the eigenvector associated with the zero eigenvalue can always be chosen real. For higher-dimensional systems, the term eigenvector should be replaced by eigenspace for mathematical accuracy.

Note that Fig. 2 is plotted for $w_{\mathrm{F}}=0$ (matter). Taking $w_{\mathrm{F}}=1/3$ (radiation) instead does not alter the qualitative structure of the phase portrait [167].

If $I_3$ is an unstable node, the following ESHU could still be realized in principle, but only at the expense of highly fine-tuned initial conditions. This scenario is not considered in the present work.

This is nontrivial, as the $n=3$ case will be shown below to fail to solve the horizon problem.

This is a rough description. Section 3c3 presents a model, in which $a$ can change by a large factor.

One exception is warm inflation, in which radiation is produced continuously during inflation, allowing thermal fluctuations to be important [297, 298].

A caveat should be noted. In the emergent Universe, the Einstein-static phase is imposed as a past-asymptotic initial configuration rather than being obtained as a generic dynamical attractor. It therefore represents a fine-tuned initial state [197]. In this sense, the horizon problem is not solved dynamically, but shifted to the assumption of a fine-tuned homogeneous initial state. By contrast, ESHU corresponds to a saddle-type attractor and can emerge from an ordinary expanding epoch.

The $z_1$ axis corresponds to the circular arc in Fig. 2. This expectation hinges on the eigenvectors being precisely aligned with the coordinate axes. If such an alignment is absent, it may be sufficient to require only that $I_3$ be a saddle to realize ESHU.

Mathematically, adopting the conventions in Appendix pp1-s2, Eq. (3.29) gives $\mathrm{d}{u}_c/\mathrm{d}\tilde{\sigma }$ identically equal to zero, where $u_c\propto 3(w_{\mathrm{F}}-1)z_0+(2-6w_{\mathrm{F}})(z_1-1/{\tilde{c}}_{35})$. This is similar to the behavior of Eq. (3.25) along its center manifold, i.e., the $z_1$ axis [see discussions below Eq. (3.26)]. These two dynamical systems impose opposite stability requirements along the center manifold. One requires the center-manifold direction to be repelling, whereas the other requires it to be attracting. This distinction makes the vanishing flow fatal for Eq. (3.29), but not for Eq. (3.25).

This conclusion applies to the case of constant $w_{\mathrm{F}}$. A time-varying $w_{\mathrm{F}}$ could lead to different dynamical behavior and deserves separate study.

This value is essentially arbitrary and does not affect the subsequent calculations or the physical conclusions.

As shown in Fig. 8, the contribution from $[N_{\mathrm{S}.\mathrm{sta}},N_{\mathrm{S}.\max }]$ is smaller than that from $[N_{\mathrm{S}.\max },N_{\mathrm{S}.\mathrm{end}}]$, and thus is also negligible.

We have not yet demonstrated this conclusion rigorously in numerical calculations due to numerical instabilities. Preliminary calculations [167] suggest that, for the illustrative trajectories shown as the orange dashed curve in Fig. 2, an extremely large $x_{3,\max }$ may be obtained with $x_{3,\mathrm{ini}}\gtrsim \mathcal{O}(1)$ and $x_{1,\mathrm{ini}}\gg 1$.

In light of the mismatch between the analytic approximation and the numerical evolution in case 1b discussed later, it is necessary to point out that the $n=3$ and $n=2$ examples displayed here are chosen on the same side in the phase portrait as case 1a, corresponding to the orange dashed curve in Fig. 2. The opposite side has not been tested for these parameter points.

This normalization is arbitrary, and the present choice is made for convenience in the discussion below.

In principle, ${\beta }_{\mathrm{c}.\mathrm{m}}$. does not depend on the explicit form of $\varphi (u_c)$. To see this, start from $\mathrm{d}{u}_c/\mathrm{d}\tilde{\sigma }=g_c(u_h,u_c)$, which follows from Eq. (a5), and expand $g_c(u_h,u_c)$ near the origin. The first nonzero terms are quadratic, of the form $u_h^2$, $u_h{u}_c$, and $u_c^2$. Since the center-manifold expansion implies $u_h\propto u_c^2$ near the origin, the terms involving $u_h$ are higher order. Therefore, $\varphi (u_c)$ does not enter the leading term of $g_c(\varphi (u_c),u_c)$.

Note that the eigenvector associated with ${\tilde{\lambda }}_{I_3,+}=0$ in the system $\{\mathrm{d}{z}_i/\mathrm{d}\tilde{\sigma }\}$ is not aligned with the $z_0$ axis. The approximation $z_0\approx u_c$ reflects the relative scaling between the $u_c$ and $z_0$ coordinates, since the transformation matrix in Eq. (3.43) includes a stretching.

Strictly, Eq. (4.6) should be understood as evaluating ${\tilde{h}}_{\mu \nu }$ at $\tilde{x}$ and all untilded quantities on the right-hand side at $x$. The reason is that, in the passive view, a gauge transformation compares the same physical point described in two coordinate systems [232].

It is worth emphasizing the distinction between $S_{\mu }$ and the electromagnetic vector potential $A_{\mu }$ regarding gauge freedom. In Maxwell theory, the internal $U(1)$ gauge transformation (not a coordinate transformation) $A_{\mu }\mapsto A_{\mu }+{\partial }_{\mu }\chi$ leaves the physical field strength $F_{\mu \nu }={\partial }_{\mu }{A}_{\nu }-{\partial }_{\nu }{A}_{\mu }$ invariant [299]. By contrast, no analogous gauge freedom has been identified for $S_{\mu }$, and thus no internal gauge-fixing condition is imposed in our analysis.

In principle, this boundary condition should be viewed as an assumption of the present Minkowski-based perturbative setup, rather than as an established default for local systems in GCT (see the discussion at the beginning of Sec. 4).

This reconstruction is not unique. Here it is used only as a diagnostic example, and our final conclusion will rely on the original second-order system rather than any particular reconstructed Lagrangian.

Note that the vacuum condition does not preclude a point mass, as in the Schwarzschild solution of general relativity.

Because the constant part of $\nu$ can be absorbed by rescaling the time coordinate in Eq. (5.1), the branch ${\nu }^{\prime }=0$ can be fixed to $\nu =0$ without loss of generality.

Solving Eq. (5.6) with Maple (Version 2022) yields this transformation (instead of an explicit solution).

Note that, for numerical integration, Eq. (5.6) is preferable to Eq. (5.7) for two reasons. It avoids the need to reconstruct $r$ and remains well behaved at points where ${\lambda }^{\prime }$ vanishes, which indeed occur in the solutions discussed below.

Section 5b2 shows that the numerical $\mathrm{\Lambda }$-like behavior does not extend smoothly to $r=0$. The expansion about $r=0$ is therefore used only as an analytic probe of the $\mathrm{\Lambda }$-like branch, not as a faithful continuation of the plotted profile to the origin.

For the regular analytic branch considered here, if odd powers are included in a more general series, the differential equation itself forces their coefficients to vanish.

If $\lambda$ is taken too large when ${\lambda }^{\prime }=0$, the approximately constant-curvature region may shrink to a very small interval, rendering the $\mathrm{\Lambda }$-like behavior barely visible. Even so, the de Sitter-related terminology is retained in the main text.

If one assumes that the asymptotic limit of $\mu$ exists and remains finite, then imposing $\mathrm{d}\mu /\mathrm{d}\lambda =0$ in Eq. (5.7) yields the candidate limits $\{0,1/2\}$ as $\lambda \to +\infty$ and $\{0,(1\pm \sqrt{2})/2\}$ as $\lambda \to -\infty$. Although this argument is not a rigorous proof, it provides a useful consistency check on the numerical behavior.

Such near-$\lambda$-axis behavior may be useful for diagnosing candidate horizons in the more general system (5.8), where $\nu$ is allowed to vary.

The boundary data are taken to be $(\lambda ,\mu ,\eta )=(2,-0.5,\pm 0.5)$ and $(2,-0.5,\pm 1.5)$ for the red trajectories, $(0,-0.2,\pm 1)$ for the green trajectories, and $(0,0.1,\pm 0.9)$ and $(2,0.5,\pm 2.2)$ for the blue trajectories.

Again, this is not a rigorous proof, since it assumes that $\mu$ and $\eta$ approach finite constants as $\lambda \to -\infty$.

The present analysis does not establish whether these are physical singularities.

The Schwarzschild metric corresponds to the first branch and is discussed in Sec. 5c2.

Equation (5.8) provides a rigorous proof of this symmetry. Note that this symmetry disappears if $\mathrm{\Lambda }\ne 0$.

This is trivial—we can set the Schwarzschild radius to be negative in the Schwarzschild metric to construct a naked singularity.

See Eq. (2.11) in Ref. [111] for an example.

A center manifold is not necessarily unique [120]. This does not affect the local stability analysis because the required series expression is unique and can be fixed order by order. Possible differences between center manifolds are flat at the critical point, e.g., of order $\exp (-1/u_c^2)$.

When transformed back to the original $x_i$ coordinates, the same center manifold may contain constant, linear, or even fractional powers. Here fractional powers can appear if the curve is parametrized by an inappropriate variable, e.g., rewriting $y\propto x^2$ as $x\propto y^{1/2}$. This is why the initial series ansatz should be imposed in the $u_i$ coordinates, rather than the $x_i$ coordinates.

One of $P_m$ and $Q_m$ may vanish if the degrees of $P$ and $Q$ differ.

This $\sigma$ is distinct from the one introduced in Eq. (a8). The same symbol is reused, with its meaning fixed by the dynamical system under consideration.