Variational Closure for Causal Bulk Viscosity in Teleparallel f(T,TG) Gravity with Unified Inflation and Late-Time Acceleration

1. Introduction

Bulk-viscous cosmology and modified gravity offer two complementary routes to explain accelerated phases of the Universe without invoking additional particulate sectors  [1,2,3,4]. In the viscous approach, an effective single fluid acquires a negative bulk pressure that can mimic dark-energy-like behavior, while causal frameworks such as the truncated Israel–Stewart theory avoid the non-hyperbolic pathologies of Eckart-type laws  [5,6,7,8]. In teleparallel gravity, torsion-based extensions of General Relativity—and, in particular, the teleparallel Gauss–Bonnet class $f(T,T_G)$—modify the background expansion and the relation between metric potentials and matter inhomogeneities [5,6,7,8,50]. A unified scenario where early and late acceleration emerge within the same framework requires both a consistent gravitational sector and a thermodynamically sound, causal description of dissipation.

A persistent gap in this program is the absence of a closure that links causal bulk-viscous relaxation to teleparallel geometric invariants at the level of first principles [1,3,10,13,14]. Most existing studies either adopt phenomenological ansätze for the viscous sector [2,6,8] or treat geometry and dissipation as independent ingredients. This work addresses that gap.

Motivation and Related Work

Relativistic bulk viscosity evolved from first-order (Eckart) closures—non-hyperbolic and acausal—to causal Israel–Stewart (IS) and extended irreversible thermodynamics (EIT). In cosmology, many studies adopt $\Pi =-3\zeta H$ with ad hoc $\zeta \left(H\right)$, detached from a variational principle. In parallel, teleparallel gravity—especially $f\left(T\right)$ and $f(T,T_G)$—modifies background expansion and linear couplings via torsion invariants. However, these two lines have largely progressed in parallel: teleparallel sectors are often treated without an explicit causal viscous closure, and viscous cosmologies rarely couple dissipation to teleparallel invariants in a controlled way. This motivates a closure that is variationally founded, hyperbolic, and geometry-aware [51].

Positioning Relative to Prior Approaches

Versus Eckart/first order: our closure is hyperbolic and causal by construction, avoiding non-hyperbolic pathologies. Versus standard IS/EIT: we retain IS relaxation but add a geometry-informed drive $\alpha \zeta K_G(T,T_G)$ from a local variational functional; entropy production and characteristic speeds satisfy simple sufficient bounds enforced as hard priors. Versus teleparallel reconstructions without viscosity: the closure enables a unified inflation–acceleration background and yields testable linear-level signatures (growth damping and $E_G$ shifts) beyond distances.

Contributions. We formulate a variationally closed causal bulk-viscous law for a single cosmological fluid evolving in $f(T,T_G)$. The closure, denoted YKGC (Yıldız-Kaykı-Güdekli-Chattopadhyay), is derived from a local Rayleigh–Onsager dissipation functional that couples the viscous pressure to a dimensionless teleparallel kernel ${\mathcal{K}}_G(T,T_G)$. The resulting relaxation equation is hyperbolic and thermodynamically consistent. We establish two sufficient conditions: (i) non-negative entropy production under explicit bounds on the kernel coupling, and (ii) subluminal characteristics for finite relaxation time, which together ensure causal, well-posed evolution. On the background, we employ a two-plateau $H\left(N\right)$ profile to reconstruct $f(T,T_G)$ models in which the same effective fluid accounts for quasi-de Sitter inflation and late-time acceleration  [15,16,17], while smoothly recovering the GR limit at early times. We analyze fixed points and delineate conditions that avoid finite-time singularities  [18,19]. At the linear level we derive $G_{\mathrm{eff}}(k,a)$ and the metric slip, finding a controlled, scale-dependent damping of structure growth that can lower $f{\sigma }_8\left(z\right)$ within thermodynamically allowed priors  [18,19]. These elements identify concrete observational targets in SN+BAO+$H\left(z\right)$ and RSD datasets  [23,24,25,26].

Scope and assumptions. We work in spatially flat FRW and adopt the diagonal tetrad ${e^A}_{\mu }=\mathrm{diag}(1,a,a,a)$ within the covariant teleparallel formulation; for FRW the vanishing spin connection (Weitzenböck gauge) is admissible and free of frame–dependence issues [27,28]. The matter sector is modeled as a single barotropic fluid with effective pressure $p_{\mathrm{eff}}=p+\Pi$, where $\Pi$ is determined by the YKGC closure in Section 3. On FRW the teleparallel invariants reduce to $T=-6H^2$ and $T_G=24H^2(\dot{H}+H^2)$  [13]. Limits to the truncated Israel–Stewart theory (vanishing kernel coupling) [7] and to GR (linear f with early-time constraints) are [9,10] stated explicitly in Secs. Section 3 and Section 2. Throughout we impose $\zeta \left(H\right)>0$, $\tau \left(H\right)>0$, and a finite local temperature $T_{\mathrm{loc}}>0$, and we enforce consistency with early-time bounds (BBN/CMB) [29] and with gravitational-wave speed constraints [30] via the priors specified in Section 4.

Organization.Section 2 summarizes the teleparallel $f(T,T_G)$ background equations and conventions. Section 3 introduces the YKGC variational closure, presents the dissipation functional, and proves the two sufficient conditions. Section 4 details the $H\left(N\right)$ reconstruction and parameter priors consistent with thermodynamics and early-time bounds. Section 5 analyzes fixed points and finite-time singularities. Section 6 derives $G_{\mathrm{eff}}(k,a)$ and the slip, and outlines observational targets. We conclude in Section 9.

Notation and Conventions

We use signature $(-,+,+,+)$, set $c=1$, and define ${\kappa }^2=8\pi G$. Overdots denote derivatives with respect to cosmic time t, and $N=lna$ is the number of e-folds. The expansion scalar is $\Theta \equiv {\nabla }_{\mu }{u}^{\mu }=3H$. Teleparallel derivatives of f are denoted $f_T\equiv \partial f/\partial T$ and $f_{T_G}\equiv \partial f/\partial T_G$. Entropy production density is $T{\nabla }_{\mu }{s}^{\mu }$. When needed, we write the effective viscous sound speed as $c_{\mathrm{visc}}^2\equiv \zeta /\left[\tau (\rho +p)\right]$.

2. Framework: Teleparallel $f(T,T_G)$ with a Single Viscous Fluid

2.1. Geometry, Tetrad and Background

We work in spatially flat FRW with line element $d{s}^2=-d{t}^2+a^2\left(t\right)d{\overrightarrow{x}}^2$ and adopt the diagonal tetrad ${e^A}_{\mu }=\mathrm{diag}(1,a,a,a)$, which is compatible with homogeneity and isotropy. We employ the covariant teleparallel formulation; for FRW, the diagonal tetrad with vanishing spin connection (Weitzenböck gauge) is admissible and avoids frame–dependence issues [27,28]. We set signature $(-,+,+,+)$ and $c=1$, and define ${\kappa }^2\equiv 8\pi G$. The Hubble rate is $H\equiv \dot{a}/a$, with overdots denoting derivatives with respect to cosmic time t.

2.2. Teleparallel Invariants

In teleparallel gravity the Weitzenböck connection has vanishing curvature and nonzero torsion [9]. On the FRW background the torsion scalar T and the teleparallel Gauss–Bonnet invariant $T_G$ reduce to

$$T=-6H^2,T_G=24H^2(\dot{H}+H^2).$$

(1)

We consider the action

$$S=\int d^{4}xe\left[{\displaystyle \frac{1}{2{\kappa }^2}}f(T,T_G)+{\mathcal{L}}_{\mathrm{m}}\right],$$

(2)

with $e\equiv det\left({e^A}_{\mu }\right)=a^3$ on FRW and ${\mathcal{L}}_{\mathrm{m}}$ the matter Lagrangian [9,10].

2.3. Effective Fluid and Conservation

The matter sector is modeled as a single barotropic fluid corrected by bulk viscosity:

$$p_{\mathrm{eff}}\equiv p+\Pi ,\Pi \mathrm{is}\mathrm{dynamically}\mathrm{determined}\mathrm{and}\mathrm{can}\mathrm{be}\mathrm{negative},$$

(3)

where p is the equilibrium pressure and $\Pi$ the bulk viscous pressure. The total energy density $\rho$ obeys

$$\dot{\rho }+3H(\rho +p+\Pi )=0.$$

(4)

Causality and thermodynamic consistency for $\Pi$ will be enforced in Section 3 via the YKGC variational closure.

2.4. Background Field Equations (Compact Form)

Variation of (2) with respect to the tetrad yields modified Friedmann equations that can be written as

$$\begin{array}{cc}\hfill 3H^2& ={\kappa }^2\left(\rho +{\rho }_f\right),\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill -2\dot{H}& ={\kappa }^2\left(\rho +p+\Pi +{\rho }_f+p_f\right),\hfill \end{array}$$

(6)

where $({\rho }_f,p_f)$ are the effective geometric contributions generated by $f(T,T_G)$ and its derivatives $(f_T,f_{T_G},{\dot{f}}_T,{\dot{f}}_{T_G},\dots )$ evaluated on (1)[10,13]. Their explicit expressions are lengthy and will be listed in Appendix A. Equations (5)–(6) together with (4) fully determine the homogeneous dynamics once the viscous closure for $\Pi$ is specified.

2.5. Limits and Consistency Conditions

We will explicitly recover: (i) the GR limit for suitable choices of f (linear in T with a vanishing $T_G$-sector) [9,10], and (ii) the standard truncated Israel–Stewart theory for vanishing kernel coupling in the viscous closure (Section 3) [7]. Throughout we assume a positive local temperature $T_{\mathrm{loc}}>0$, $\zeta \left(H\right)>0$ and $\tau \left(H\right)>0$ (defined in Section 3). We further impose consistency with early-time bounds (BBN/CMB) [29] and with observational constraints on the gravitational-wave propagation speed [30], which translate into parameter restrictions specified below.

3. YKGC Variational Closure for Causal Bulk Viscosity

3.1. Rayleigh–Onsager Dissipation and Geometric Kernel

We promote the truncated Israel–Stewart bulk-viscous law to a variationally closed form by specifying a local Rayleigh–Onsager dissipation functional [31,32,33,34,35] that couples the viscous pressure $\Pi$ to a dimensionless teleparallel kernel built from $f(T,T_G)$:

$$\mathcal{R}[\Pi ;T,T_G]={\displaystyle \frac{1}{2}}{\displaystyle \frac{{\Pi }^2}{\zeta \left(H\right)}}-\alpha \Pi {\mathcal{K}}_G(T,T_G),$$

(7)

where $\zeta \left(H\right)>0$ is the bulk-viscosity coefficient and $\alpha$ is a dimensionless coupling. To ensure differentiability and regular behavior near $H\to 0$ (e.g., in bounce-like regimes), we introduce a smooth dimensionless regulator via

$${\widehat{f}}_{T_G}\equiv H_{★}^2{f}_{T_G},Q\equiv {\displaystyle \frac{T_G}{T^2+\epsilon H_{★}^4}},$$

and define

$$S_G\equiv 1+{\lambda }_1\sqrt{f_T^2+{\epsilon }^2}+{\lambda }_2\sqrt{Q^2{\widehat{f}}_{T_G}^2+{\epsilon }^2}.$$

The geometric kernel is then

$$\begin{array}{|c|}\hline {\mathcal{K}}_G\equiv u^{\mu }{\nabla }_{\mu }ln{S}_G\\ \hline\end{array}$$

(8)

with $f_T\equiv \partial f/\partial T$, $f_{T_G}\equiv \partial f/\partial T_G$, small $\epsilon >0$ (we set ${10}^{-12}\le \epsilon \le {10}^{-8}$ in numerics), and a characteristic Hubble scale $H_{★}$ (e.g., $H_0$ or an inflationary plateau). On a spatially flat FRW background, $T=-6H^2$ and $T_G=24H^2(\dot{H}+H^2)$. With the above definitions the logarithm argument $S_G$ is manifestly dimensionless and strictly positive for $\epsilon >0$. In the GR limit $f(T,T_G)=T$ one has $f_T=1$ and $f_{T_G}=0$, hence $S_G\to \mathrm{const}.$ and therefore ${\mathcal{K}}_G=u^{\mu }{\nabla }_{\mu }ln{S}_G=0$. Here $H_{★}$ is a fixed reference Hubble scale introduced solely for nondimensionalization.

Figure 1. Regulated geometric kernel and regulator components for the non-trivial benchmark ($\gamma =0.2$, $p=1.2$, $q=1$). Left axis (linear): ${\lambda }_1\sqrt{f_T^2+{\epsilon }^2}$ (red) and ${\lambda }_2\sqrt{Q^2{\widehat{f}}_{T_G}^2+{\epsilon }^2}$ (green), which sum to $S_G-1$. Right axis (signed logarithm): ${\mathcal{K}}_G\equiv d(ln{S}_G)/dt=Hd(ln{S}_G)/dN$ (blue). Definitions used in the panel: $T=-6H^2$, $T_G=24H^2(\dot{H}+H^2)$, and $Q=T_G/(T^2+\epsilon H_{★}^4)$ entering $f(T,T_G)=T+\mu {(-T)}^p+\nu {(-T)}^q\sqrt{1+\gamma Q}$. The inset box lists the exact parameter values (two-plateau background $H\left(N\right)$ and regulator constants) used to generate this figure.

Figure 1. Regulated geometric kernel and regulator components for the non-trivial benchmark ($\gamma =0.2$, $p=1.2$, $q=1$). Left axis (linear): ${\lambda }_1\sqrt{f_T^2+{\epsilon }^2}$ (red) and ${\lambda }_2\sqrt{Q^2{\widehat{f}}_{T_G}^2+{\epsilon }^2}$ (green), which sum to $S_G-1$. Right axis (signed logarithm): ${\mathcal{K}}_G\equiv d(ln{S}_G)/dt=Hd(ln{S}_G)/dN$ (blue). Definitions used in the panel: $T=-6H^2$, $T_G=24H^2(\dot{H}+H^2)$, and $Q=T_G/(T^2+\epsilon H_{★}^4)$ entering $f(T,T_G)=T+\mu {(-T)}^p+\nu {(-T)}^q\sqrt{1+\gamma Q}$. The inset box lists the exact parameter values (two-plateau background $H\left(N\right)$ and regulator constants) used to generate this figure.

3.2. Causal Relaxation Equation (YKGC Law)

Stationarity of (7) in the GENERIC/Onsager sense yields the constitutive evolution of the bulk pressure:

$$\begin{array}{|c|}\hline \tau \left(H\right)u^{\mu }{\nabla }_{\mu }\Pi +\Pi =-\zeta \left(H\right)\Theta +\alpha \zeta \left(H\right){\mathcal{K}}_G(T,T_G)\\ \hline\end{array}$$

(9)

where $\Theta \equiv {\nabla }_{\mu }{u}^{\mu }=3H$ on FRW and $\tau \left(H\right)>0$ is the relaxation time. Equation (9) reduces to the standard truncated Israel–Stewart law in the decoupling limit $\alpha \to 0$ or when the kernel vanishes[7]. Together with Eqs. (5)–(6) and (4), this closes the homogeneous dynamics.

Relation to Classical Closures (Eckart, IS/EIT).

For reference, the three constitutive options most commonly used in cosmology read:

$$\begin{array}{cccccc}\hfill & \mathrm{Eckart}(\mathrm{first}\mathrm{order})\hfill & \hfill & :\hfill & \hfill & \Pi =-\zeta \Theta \left(\mathrm{non}-\mathrm{hyperbolic},\mathrm{acausal}\right),\hfill \\ \hfill & \mathrm{truncated}\mathrm{Israel}--\mathrm{Stewart}\left(\mathrm{IS}\right)\hfill & \hfill & :\hfill & \hfill & \tau u^{\mu }{\nabla }_{\mu }\Pi +\Pi =-\zeta \Theta \left(\mathrm{hyperbolic}\right),\hfill \\ \hfill & \mathbf{YKGC}\left(\mathbf{this}\mathbf{work}\right)\hfill & \hfill & :\hfill & \hfill & \tau u^{\mu }{\nabla }_{\mu }\Pi +\Pi =-\zeta \Theta +\alpha \zeta K_G(T,T_G),\hfill \end{array}$$

where $\Theta \equiv {\nabla }_{\mu }{u}^{\mu }$ and $K_G$ is the teleparallel drive from Sec. 3.1. The full Müller–Israel–Stewart / EIT frameworks admit additional second-order terms (e.g. $\mathcal{O}\left(\Pi \Theta \right)$, gradients of $T_{\mathrm{loc}}$); our implementation follows the truncated IS sector but augments it by the geometry-informed drive $\alpha \zeta K_G$ obtained variationally. In the decoupling limit $\alpha \to 0$ (or $K_G\to 0$) we recover the truncated IS law, whereas $f(T,T_G)\to T$ implies $K_G\to 0$ and the GR background.

3.3. Thermodynamic Positivity (Sufficient Condition)

In natural units ($c=\hslash =k_B=1$), the local entropy production density reads

$$T_{\mathrm{loc}}{\nabla }_{\mu }{s}^{\mu }={\displaystyle \frac{{\Pi }^2}{\zeta \left(H\right)}}-\alpha \Pi {\mathcal{K}}_G(T,T_G).$$

(10)

A sufficient condition for non-negative entropy production is

$$\begin{array}{|c|}\hline \left|\alpha {\mathcal{K}}_G(T,T_G)\right|\le {\displaystyle \frac{\left|\Pi \right|}{\zeta \left(H\right)}}\\ \hline\end{array}.$$

(11)

Justification. From (10), $T_{\mathrm{loc}}{\nabla }_{\mu }{s}^{\mu }\ge {\displaystyle \frac{{\Pi }^2}{\zeta }}-\left|\alpha \right|\left|\Pi \right||{\mathcal{K}}_G|$, which is $\ge 0$ if (11) holds. In practice we enforce (11) as a prior on $(\alpha ,{\lambda }_{1,2})$ during reconstruction and parameter inference.

Sharpness and Necessity

Condition (11) is sharp in the following pointwise sense. Writing $T_{\mathrm{loc}}{\nabla }_{\mu }{s}^{\mu }=({\Pi }^2/\zeta )-\alpha \Pi {\mathcal{K}}_G$, non-negativity at a given state $(\Pi ,{\mathcal{K}}_G)$ with no sign assumptions is equivalent to $|\alpha {\mathcal{K}}_G|\le |\Pi |/\zeta$. Indeed, the quadratic form $(1/\zeta ){\Pi }^2-\alpha \Pi {\mathcal{K}}_G$ attains its minimum at $\Pi =\left(\alpha \zeta {\mathcal{K}}_G\right)/2$, yielding a non-positive minimum unless (11) holds; saturation of (11) yields equality of the production. Thus (11) is the minimal uniform bound that guarantees $T_{\mathrm{loc}}{\nabla }_{\mu }{s}^{\mu }\ge 0$ along the realized trajectory. Tighter choices $|\alpha {\mathcal{K}}_G|\le \chi |\Pi |/\zeta$ with $0<\chi <1$ provide stricter dissipation priors.

3.4. Causality and Hyperbolicity (Sufficient Condition)

Linearizing (9) around a homogeneous background yields a first-order hyperbolic equation for $\delta \Pi$ with characteristic speed

$$c_{\mathrm{visc}}^2\equiv {\displaystyle \frac{\zeta \left(H\right)}{\tau \left(H\right)(\rho +p)}}.$$

(12)

Lemma 2 (Subluminal characteristics). If $\tau \left(H\right)>0$, $\zeta \left(H\right)>0$ and

$$\begin{array}{|c|}\hline c_{\mathrm{visc}}^2\le 1\\ \hline\end{array},$$

(13)

then the viscous sector has subluminal characteristics and the coupled system is (locally) well-posed. Sketch. The principal symbol has real eigenvalues $\pm c_{\mathrm{visc}}$; the bound (13) confines the characteristic cone within the light cone.

On Necessity vs Sufficiency (Causality)

The positivity $\tau \left(H\right)>0$, $\zeta \left(H\right)>0$ is sufficient for hyperbolicity of $\Pi$ in isolation. The subluminality prior $c_{\mathrm{visc}}^2\le 1$ in (13) is not strictly necessary for hyperbolicity, but it is the natural necessary and sufficient pointwise bound for keeping the viscous characteristic cone inside the light cone when $c_{\mathrm{visc}}^2=\zeta \left(H\right)/\left[\tau \left(H\right)(\rho +p)\right]$ parameterizes the principal symbol of the $\Pi$-sector. Couplings to the metric can only lower the effective propagation speed in our setup; hence adopting (13) as a hard prior is conservative.

3.5. Limits and Consistency Checks

IS Limit

$\alpha \to 0$ (or ${\lambda }_{1,2}\to 0$) reduces (9) to the truncated Israel–Stewart relaxation $\tau u^{\mu }{\nabla }_{\mu }\Pi +\Pi =-\zeta \Theta$.

GR Limit

For $f(T,T_G)=T$ (and vanishing $T_G$ sector) one has $f_T=1$, $f_{T_G}=0$, hence ${\mathcal{K}}_G=0$ and the background reduces to GR.

Early-Time Bounds and GW Speed

Viable parameter choices are further restricted by BBN/CMB early-time limits and by consistency with gravitational-wave speed constraints; these conditions enter as priors in Section 4.

3.6. Parametrization and Priors

To minimize arbitrariness while retaining sufficient flexibility, we adopt

$$\begin{array}{cc}\hfill \zeta \left(H\right)& ={\zeta }_0+{\zeta }_{1}H,{\zeta }_0>0,{\zeta }_1\ge 0,\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill \tau \left(H\right)& ={\tau }_0+{\tau }_1/H,{\tau }_0>0,{\tau }_1\ge 0,\hfill \end{array}$$

(15)

and constrain $({\zeta }_0,{\zeta }_1,{\tau }_0,{\tau }_1)$ jointly with $(\alpha ,{\lambda }_{1,2})$ by the thermodynamic bound (11), the causal bound (13), and early-time priors. These choices keep the total number of free parameters compact and support both quasi-de Sitter inflation and a late-time de Sitter attractor.

3.7. Practical Remarks

(i) The regularized kernel (8) is well defined for all finite H and admits a finite $H\to 0$ limit via the regulator; equivalently, one may take the limit using $d/dN$ with $N\equiv lna$. (ii) Because ${\mathcal{K}}_G$ enters multiplied by $\alpha \zeta \left(H\right)$ in (9), the entropy bound can be imposed on the composite parameter $\alpha \zeta$ if desired. (iii) The background and linear-perturbation equations with (9) remain first order in time for $\Pi$, preserving hyperbolicity under (13).

Physical Picture

The YKGC law augments the truncated Israel–Stewart relaxation by a geometry-informed drive $\alpha \zeta {\mathcal{K}}_G(T,T_G)$ built from teleparallel invariants. Intuitively, ${\mathcal{K}}_G$ measures how rapidly a geometry-sourced scalar changes along the flow; when ${\mathcal{K}}_G>0$ during accelerated phases, it drives a more negative bulk pressure $\Pi <0$, acting as a controlled thermodynamic brake that can suppress structure growth while keeping the background close to $\Lambda$CDM under early-time and GW-speed consistency. The couplings $({\lambda }_1,{\lambda }_2)$ set the relative weight of $f_T$ and $f_{T_G}$ channels inside the kernel, while $\alpha$ controls how strongly geometry feeds into dissipation. In the viscous sector, $({\zeta }_0,{\zeta }_1)$ set the bulk-damping scale and its H-dependence, and $({\tau }_0,{\tau }_1)$ set relaxation and hence $c_{\mathrm{visc}}^2=\zeta /\left[\tau (\rho +p)\right]$.

4. Background Reconstruction and Parameter Priors

4.1. Two-Plateau $H\left(N\right)$ Profile

We reconstruct $f(T,T_G)$ from a prescribed Hubble history with two plateaus,

$$H\left(N\right)={\displaystyle \frac{H_{\inf }}{1+e^{\alpha (N-N_{★})}}}+H_{\Lambda }\left[1-{\displaystyle \frac{1}{1+e^{\beta (N-N_0)}}}\right],$$

(16)

with $\alpha ,\beta >0$, $N_{★}$ near inflation exit and $N_0$ near the matter–acceleration transition. On FRW,

$$T\left(N\right)=-6H^2\left(N\right),T_G\left(N\right)=24H^2\left(N\right)\left(H^2\left(N\right)+H\left(N\right)H^{\prime }\left(N\right)\right).$$

(17)

Figure 2. Background scalars with exact parameter values. Left axis: $H/H_0$ (blue). Right axis: signed logarithms $sgn·{log}_{10}\left(1+|-T/H_0^2|\right)$ (red) and $sgn·{log}_{10}\left(1+|T_G/H_0^4|\right)$ (green). The dashed line marks the transition center $N_{★}$. Definitions: $T=-6H^2$, $T_G=24H^2(\dot{H}+H^2)$, and $N\equiv lna$. All quantities are nondimensionalized by $H_0$. The inset box lists the exact parameter values used in this figure.

Figure 2. Background scalars with exact parameter values. Left axis: $H/H_0$ (blue). Right axis: signed logarithms $sgn·{log}_{10}\left(1+|-T/H_0^2|\right)$ (red) and $sgn·{log}_{10}\left(1+|T_G/H_0^4|\right)$ (green). The dashed line marks the transition center $N_{★}$. Definitions: $T=-6H^2$, $T_G=24H^2(\dot{H}+H^2)$, and $N\equiv lna$. All quantities are nondimensionalized by $H_0$. The inset box lists the exact parameter values used in this figure.

4.2. Reconstruction Ansatz for $f(T,T_G)$

To capture the background with minimal parameters we take

$$f(T,T_G)=T+\mu {(-T)}^p+\nu {(-T)}^q\sqrt{1+\gamma {\displaystyle \frac{T_G}{T^2+\epsilon H_{★}^4}}},$$

(18)

with $(\mu ,\nu ,p,q,\gamma )$ and the same regulator as in Section 3. The GR limit is $(\mu ,\nu )=0$. The geometric sector $({\rho }_f,p_f)$ evaluated on (17) enters Eqs. (5)–(6).

Figure 3. PRD ansatz for $f(T,T_G)$ and its derivatives. Left axis: $f(T,T_G)$ (blue). Right axis (signed logarithms): $f_T$ (red solid), $f_{TT}$ (red dashed), $f_{T{T}_G}$ (red dotted), $f_{T_G}$ (green solid), and $f_{T_G{T}_G}$ (green dashed). The ansatz is $f(T,T_G)=T+\mu {(-T)}^p+\nu {(-T)}^q\sqrt{1+\gamma Q}$ with $Q\equiv T_G/(T^2+\epsilon H_{★}^4)$, evaluated on $T=-6H^2$ and $T_G=24H^2(\dot{H}+H^2)$. In the approved benchmark ($\gamma =0$, $p=q=1$) the $T_G$-sector derivatives vanish as expected. The inset box lists the exact parameter values used in this figure.

Figure 3. PRD ansatz for $f(T,T_G)$ and its derivatives. Left axis: $f(T,T_G)$ (blue). Right axis (signed logarithms): $f_T$ (red solid), $f_{TT}$ (red dashed), $f_{T{T}_G}$ (red dotted), $f_{T_G}$ (green solid), and $f_{T_G{T}_G}$ (green dashed). The ansatz is $f(T,T_G)=T+\mu {(-T)}^p+\nu {(-T)}^q\sqrt{1+\gamma Q}$ with $Q\equiv T_G/(T^2+\epsilon H_{★}^4)$, evaluated on $T=-6H^2$ and $T_G=24H^2(\dot{H}+H^2)$. In the approved benchmark ($\gamma =0$, $p=q=1$) the $T_G$-sector derivatives vanish as expected. The inset box lists the exact parameter values used in this figure.

4.3. YKGC Law in N-Time and Closure Parameters

We evolve in $N=lna$:

$$\tau \left(H\right)H{\Pi }^{\prime }+\Pi =-3\zeta \left(H\right)H+\alpha \zeta \left(H\right){\mathcal{K}}_G,$$

(19)

with ${\mathcal{K}}_G$ from (8). For the viscous functions we use

$$\begin{array}{cc}\hfill \zeta \left(H\right)& ={\zeta }_0+{\zeta }_{1}H,{\zeta }_0>0,{\zeta }_1\ge 0,\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill \tau \left(H\right)& ={\tau }_0+{\tau }_1/H,{\tau }_0>0,{\tau }_1\ge 0.\hfill \end{array}$$

(21)

Figure 4. YKGC closure with an analytically evaluated kernel. Left axis (linear): $\zeta \left(H\right)={\zeta }_0+{\zeta }_{1}H$ (red) and $\tau \left(H\right)={\tau }_0+{\tau }_1/H$ (green). Right axis (signed logarithms): anisotropic stress $\Pi \left(N\right)$ (blue, solid), obtained from $\tau \left(H\right)Hd\Pi /dN+\Pi =-\alpha \zeta \left(H\right){\mathcal{K}}_G\left(N\right)$ via an integrating–factor solution with $\Pi \left(N_{min}\right)=0$, and the driving term $\alpha \zeta \left(H\right){\mathcal{K}}_G\left(N\right)$ (blue, dashed). The kernel is computed by the chain rule, ${\mathcal{K}}_G={\dot{\mathcal{S}}}_G/{\mathcal{S}}_G$, with ${\mathcal{S}}_G=1+{\lambda }_1\sqrt{f_T^2+{\epsilon }^2}+{\lambda }_2\sqrt{{\left(T_G/(T^2+\epsilon H_{★}^4)\right)}^2{f}_{T_G}^2+{\epsilon }^2}$ and the PRD ansatz $f(T,T_G)=T+\mu {(-T)}^p+\nu {(-T)}^q\sqrt{1+\gamma Q}$, $Q=T_G/(T^2+\epsilon H_{★}^4)$ (non-trivial benchmark: $\gamma =0.2$, $p=1.2$, $q=1$). Definitions: $T=-6H^2$ and $T_G=24H^2(\dot{H}+H^2)$. The inset box lists the exact parameter values used in this figure.

Figure 4. YKGC closure with an analytically evaluated kernel. Left axis (linear): $\zeta \left(H\right)={\zeta }_0+{\zeta }_{1}H$ (red) and $\tau \left(H\right)={\tau }_0+{\tau }_1/H$ (green). Right axis (signed logarithms): anisotropic stress $\Pi \left(N\right)$ (blue, solid), obtained from $\tau \left(H\right)Hd\Pi /dN+\Pi =-\alpha \zeta \left(H\right){\mathcal{K}}_G\left(N\right)$ via an integrating–factor solution with $\Pi \left(N_{min}\right)=0$, and the driving term $\alpha \zeta \left(H\right){\mathcal{K}}_G\left(N\right)$ (blue, dashed). The kernel is computed by the chain rule, ${\mathcal{K}}_G={\dot{\mathcal{S}}}_G/{\mathcal{S}}_G$, with ${\mathcal{S}}_G=1+{\lambda }_1\sqrt{f_T^2+{\epsilon }^2}+{\lambda }_2\sqrt{{\left(T_G/(T^2+\epsilon H_{★}^4)\right)}^2{f}_{T_G}^2+{\epsilon }^2}$ and the PRD ansatz $f(T,T_G)=T+\mu {(-T)}^p+\nu {(-T)}^q\sqrt{1+\gamma Q}$, $Q=T_G/(T^2+\epsilon H_{★}^4)$ (non-trivial benchmark: $\gamma =0.2$, $p=1.2$, $q=1$). Definitions: $T=-6H^2$ and $T_G=24H^2(\dot{H}+H^2)$. The inset box lists the exact parameter values used in this figure.

4.4. Matching and Hard Priors

Parameters are fixed by: (i) slow-roll exit with finite $H^{\prime }$ at $N\sim N_{★}$; (ii) radiation/matter eras with $|{\rho }_f|/3H^2\ll 1$; (iii) late-time de Sitter. We impose the hard priors  [29,30]

$$\begin{array}{cc}\hfill & \left|\alpha {\mathcal{K}}_G\right|\le {\displaystyle \frac{\left|\Pi \right|}{\zeta \left(H\right)}},\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill & c_{\mathrm{visc}}^2={\displaystyle \frac{\zeta \left(H\right)}{\tau \left(H\right)(\rho +p)}}\le 1,\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill & {\left.{\displaystyle \frac{|{\rho }_f|}{3H^2}}\right|}_{\mathrm{BBN}/\mathrm{CMB}}\ll 1,{\left.{\displaystyle \frac{\left|\Pi \right|}{\rho +p}}\right|}_{\mathrm{CMB}}\ll 1,\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill & \mathrm{GW}\mathrm{speed}\mathrm{consistency}\mathrm{encoded}\mathrm{as}\mathrm{bounds}\mathrm{on}(\mu ,\nu ,\gamma ).\hfill \end{array}$$

(25)

Figure 5. Tensor-sector viability of the $F(T,T_G)$ benchmark used in Figs. 3–4. Red (left axis): effective Planck mass $M_{*}^2\left(N\right)\equiv -F_T+4H{\dot{F}}_{T_G}$; positivity signals absence of tensor ghosts (shaded if negative). Blue (right axis): extra GW friction $\Delta \gamma /H=\left[-{\dot{F}}_T+4\dot{H}{\dot{F}}_{T_G}+4H{\ddot{F}}_{T_G}\right]/\left[H(-F_T+4H{\dot{F}}_{T_G})\right]$. Green: $c_T^2-1=0$, indicating luminal propagation of gravitational waves. All background and model functions correspond to the same non-trivial set as in Figs. 3B–4; definitions: $T=-6H^2$ and $T_G=24H^2(\dot{H}+H^2)$. The inset lists the exact parameter values used in this figure.

Figure 5. Tensor-sector viability of the $F(T,T_G)$ benchmark used in Figs. 3–4. Red (left axis): effective Planck mass $M_{*}^2\left(N\right)\equiv -F_T+4H{\dot{F}}_{T_G}$; positivity signals absence of tensor ghosts (shaded if negative). Blue (right axis): extra GW friction $\Delta \gamma /H=\left[-{\dot{F}}_T+4\dot{H}{\dot{F}}_{T_G}+4H{\ddot{F}}_{T_G}\right]/\left[H(-F_T+4H{\dot{F}}_{T_G})\right]$. Green: $c_T^2-1=0$, indicating luminal propagation of gravitational waves. All background and model functions correspond to the same non-trivial set as in Figs. 3B–4; definitions: $T=-6H^2$ and $T_G=24H^2(\dot{H}+H^2)$. The inset lists the exact parameter values used in this figure.

5. Dynamical System and Avoidance of Finite-Time Singularities

5.1. Autonomous Variables and Background Equations

We recast the homogeneous dynamics in the e-fold time $N\equiv lna$ [36]. Define the compact variables

$$x\equiv {\displaystyle \frac{H}{H_0}},y\equiv {\displaystyle \frac{\Pi }{3H^2}},{\Omega }_f\equiv {\displaystyle \frac{{\kappa }^2{\rho }_f}{3H^2}},{\Omega }_m\equiv {\displaystyle \frac{{\kappa }^2\rho }{3H^2}},$$

(26)

so that the first Friedmann equation reads

$$1={\Omega }_m+{\Omega }_f.$$

(27)

The effective equation-of-state parameter is

$$w_{\mathrm{eff}}\equiv -1-{\displaystyle \frac{2}{3}}{\displaystyle \frac{H^{\prime }}{H}}=-1-{\displaystyle \frac{2}{3}}{\displaystyle \frac{x^{\prime }}{x}}.$$

(28)

Using Eqs. (5)–(6) and the conservation law (4), the evolution for $(x,y)$ becomes

$$\begin{array}{cc}\hfill x^{\prime }& =-{\displaystyle \frac{3}{2}}\left(1+w_{\mathrm{eff}}\right)x=-{\displaystyle \frac{3}{2}}\left[{\Omega }_m(1+w)+{\Omega }_f(1+w_f)+y\right]x,\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill y^{\prime }& =-\left(3+{\displaystyle \frac{H^{\prime }}{H}}\right)y-{\displaystyle \frac{1}{H}}{\displaystyle \frac{y}{\tau \left(H\right)}}+{\displaystyle \frac{1}{H}}\left[-{\displaystyle \frac{\zeta \left(H\right)}{H}}+{\displaystyle \frac{\alpha \zeta \left(H\right)}{3H^2}}{\mathcal{K}}_G\right]3,\hfill \end{array}$$

(30)

where we used the YKGC law in N-time, cf. (19), and $w\equiv p/\rho$ is the barotropic index for the equilibrium fluid. The effective geometric pressure ratio $w_f\equiv p_f/{\rho }_f$ is calculable from the chosen $f(T,T_G)$ ansatz (Section 4). Explicit forms are delegated to Appendix A; for the dynamical analysis we only require their regularity on finite H.

For numerical stability we rewrite (30) in fully dimensionless form using (14)–(15) and $H=H_{0}x$:

$$y^{\prime }=-\left(3+{\displaystyle \frac{x^{\prime }}{x}}\right)y-{\displaystyle \frac{y}{\tilde{\tau }\left(x\right)}}-\tilde{\zeta }\left(x\right)+\tilde{\alpha }\left(x\right){\mathcal{K}}_G,$$

(31)

with

$$\tilde{\zeta }\left(x\right)\equiv {\displaystyle \frac{3\zeta \left(H\right)}{H}}=3\left({\displaystyle \frac{{\zeta }_0}{H_{0}x}}+{\zeta }_1\right),\tilde{\tau }\left(x\right)\equiv {\displaystyle \frac{\tau \left(H\right)H}{1}}={\tau }_0{H}_{0}x+{\tau }_1,\tilde{\alpha }\left(x\right)\equiv {\displaystyle \frac{\alpha \zeta \left(H\right)}{H^2}}.$$

(32)

5.2. Fixed Points and Acceleration

A fixed point $(x^{*},y^{*})$ satisfies $x^{\prime }=0=y^{\prime }$. Using (28) one has

$$x^{\prime }=0\iff w_{\mathrm{eff}}^{*}=-1.$$

(33)

Hence de Sitter fixed points obey

$$y^{*}=-{\Omega }_m^{*}(1+w)-{\Omega }_f^{*}(1+w_f^{*}),$$

(34)

and from (31) at the fixed point,

$$0=-3y^{*}-{\displaystyle \frac{y^{*}}{\tilde{\tau }\left(x^{*}\right)}}-\tilde{\zeta }\left(x^{*}\right)+\tilde{\alpha }\left(x^{*}\right){\mathcal{K}}_G^{*}.$$

(35)

Equations (27), (34), and (35) determine $(x^{*},y^{*},{\Omega }_m^{*})$ given $(\mu ,\nu ,p,q,\gamma )$ and viscous parameters. Acceleration requires $w_{\mathrm{eff}}<-1/3$, equivalently $x^{\prime }/x>-1$; near the fixed point this is automatic.

5.3. Linear Stability

Linearizing $(x,y)\to (x^{*}+\delta x,y^{*}+\delta y)$ yields

$$\left(\begin{array}{c}\delta x^{\prime }\\ \delta y^{\prime }\end{array}\right)={\mathbf{J}}_{*}\left(\begin{array}{c}\delta x\\ \delta y\end{array}\right),{\mathbf{J}}_{*}\equiv {\left(\begin{array}{cc}-{\displaystyle \frac{3}{2}}(1+w_{\mathrm{eff}})-{\displaystyle \frac{3}{2}}x{\partial }_x{w}_{\mathrm{eff}}& -{\displaystyle \frac{3}{2}}x{\partial }_y{w}_{\mathrm{eff}}\\ ★& -3-{\displaystyle \frac{x^{\prime }}{x}}-{\partial }_y\left[{\displaystyle \frac{y}{\tilde{\tau }}}\right]\end{array}\right)}_{*},$$

(36)

where “★’’ collects derivatives of (31) with respect to x (including those of $\tilde{\zeta },\tilde{\tau },\tilde{\alpha }$ and ${\mathcal{K}}_G$). At a de Sitter point $w_{\mathrm{eff}}^{*}=-1$ and $x_{*}^{\prime }=0$, which simplifies (36). Stability requires $\mathrm{Tr}{\mathbf{J}}_{*}<0$ and $det{\mathbf{J}}_{*}>0$. In practice, we evaluate ${\mathbf{J}}_{*}$ numerically; analytically, a sufficient (not necessary) condition is

$${\displaystyle \frac{1}{\tilde{\tau }\left(x^{*}\right)}}>0,3+{\displaystyle \frac{1}{\tilde{\tau }\left(x^{*}\right)}}>{\partial }_y{\left[\tilde{\alpha }\left(x\right){\mathcal{K}}_G\right]}_{*},$$

(37)

which is compatible with the YKGC causality prior $\tau >0$.

5.4. Classification of Finite-Time Singularities and Avoidance

We adopt the standard classification for finite-time singularities [37]: Type I (Big Rip): $a\to \infty$, $\left|H\right|\to \infty$, ${\rho }_{\mathrm{eff}},\left|p_{\mathrm{eff}}\right|\to \infty$; Type II (Sudden): $a\to a_s$, $\left|H\right|<\infty$, ${\rho }_{\mathrm{eff}}<\infty$, $|p_{\mathrm{eff}}|\to \infty$; Type III: $a\to a_s$, $\left|H\right|\to \infty$, ${\rho }_{\mathrm{eff}},\left|p_{\mathrm{eff}}\right|\to \infty$; Type IV: $a\to a_s$, $H,{\rho }_{\mathrm{eff}},p_{\mathrm{eff}}$ finite, but higher derivatives diverge.

Sufficient Avoidance Criteria

Let $H\left(N\right)$ be continuously differentiable and the kernel regularized as in (8). Then:

(i) No Type I/III if there exists $M>0$ such that for all N,

$$\left|\tilde{\zeta }\left(x\right)\right|+\left|\tilde{\alpha }\left(x\right){\mathcal{K}}_G\right|\le M,\tilde{\tau }\left(x\right)\ge {\tau }_{min}>0,$$

(38)

and the geometric sector satisfies ${\Omega }_f$ bounded. Under these, $\left|y\right|$ remains bounded by Grönwall-type estimates applied to (31), which in turn bounds $H^{\prime }$ via (29); hence H cannot blow up in finite N.

(ii) No Type II if $\tilde{\zeta },\tilde{\tau },\tilde{\alpha }{\mathcal{K}}_G$ are continuous and

$$\underset{N\to N_s}{lim}y\left(N\right)\mathrm{finite},\underset{N\to N_s}{lim}{\partial }_N\left[\tilde{\alpha }\left(x\right){\mathcal{K}}_G\right]\mathrm{finite},$$

(39)

so that $p_{\mathrm{eff}}=p+3H^{2}y$ remains finite.

(iii) No Type IV if $H,H^{\prime }$ and the regulated kernel are $C^1$ in a neighborhood of $N_s$. This is guaranteed by the regulator in (8) together with finite $(\zeta ,\tau )$.

Bounce Option

A non-singular bounce at $N_b$ requires $H\left(N_b\right)=0$ and $H^{\prime }\left(N_b\right)>0$. With the regulator in (8) the kernel is finite at $H=0$, and (29)–(31) remain regular provided $\tilde{\tau }\left(0\right)={\tau }_1>0$ and $\tilde{\zeta }\left(0\right)=3{\zeta }_0/\left(H_{0}x\right)\to \infty$ is avoided by choosing ${\zeta }_0=0$ or by starting the bounce analysis at finite x in t-time. In practice we use the t-time YKGC form near $H=0$ or set ${\zeta }_0=0$ to keep $\tilde{\zeta }$ bounded. Summary and pointer. In the baseline analysis we keep ${\zeta }_0>0$ and ${\tau }_0>0$. The non-singular bounce is handled separately in Appendix C, where we (i) recast the YKGC law in cosmic time t, (ii) implement the controlled limit ${\zeta }_0\to 0$ with $\zeta \left(H\right)={\zeta }_{1}H$ and $\tau \left(H\right)={\tau }_1/H$ so that $0<c_{\mathrm{visc}}^2\le 1$, (iii) verify that the regulated kernel $K_G$ of Eq. (8) remains finite across $H=0$, (iv) check hyperbolicity and non-negative entropy production $T_{\mathrm{loc}}{\nabla }_{\mu }{s}^{\mu }\ge 0$, and (v) confirm that the tensor-sector priors $M_{*}^2>0$ and $c_T^2=1$ are maintained for the reported benchmarks.

5.5. Acceleration Domain and Observationally Viable Attractor

Combining (35) with the priors (22)–(25) yields a compact region in parameter space where the late-time de Sitter fixed point is stable and $w_{\mathrm{eff}}\to -1$. We map this domain numerically and find that typical trajectories starting from inflationary initial data flow into this viable basin. Within this region, the effective matter fraction ${\Omega }_{m0}$ and $H_0$ extracted from the background are consistent with late-time distance indicators; linear-growth viability is addressed in Section 6.

6. Linear Perturbations and Observables

6.1. Setup and Gauge

We work in Newtonian gauge,

$$d{s}^2=-(1+2\Phi )d{t}^2+a^2\left(t\right)(1-2\Psi )d{\overrightarrow{x}}^2,$$

(40)

see Refs. [38,39] for general treatments of cosmological perturbations. We consider scalar perturbations of the single effective fluid with density contrast $\delta \equiv \delta \rho /\rho$ and velocity divergence $\theta \equiv {\partial }_i{v}^i$ in Fourier space. The effective pressure is $p_{\mathrm{eff}}=p+\Pi$, and at first order

$$\delta p_{\mathrm{eff}}=c_s^2\delta \rho +\delta \Pi ,c_s^2\equiv {\left.{\displaystyle \frac{\partial p}{\partial \rho }}\right|}_s,$$

(41)

where $c_s^2$ is the adiabatic sound speed of the equilibrium sector. The expansion scalar $\Theta =3H$ perturbs as

$$\delta \Theta =-{\displaystyle \frac{k^2}{a^2}}v+3(\dot{\Psi }+H\Phi ),$$

(42)

with $\theta \equiv k^{2}v$.

6.2. Linearized Conservation and YKGC Closure

Energy–momentum conservation yields (for a barotropic background $w\equiv p/\rho$)

$$\begin{array}{cc}\hfill \dot{\delta }& =-(1+w)\left(\theta -3\dot{\Psi }\right)-3H\left(c_s^2-w\right)\delta -3H{\displaystyle \frac{\delta \Pi }{\rho }},\hfill \end{array}$$

(43)

$$\begin{array}{cc}\hfill \dot{\theta }& =-H(1-3w)\theta +{\displaystyle \frac{k^2}{a^2}}\Phi -{\displaystyle \frac{k^2}{a^2}}\left({\displaystyle \frac{c_s^2\delta \rho +\delta \Pi }{\rho +p}}\right).\hfill \end{array}$$

(44)

The YKGC constitutive law (9) linearizes to

$$\tau \dot{\delta \Pi }+\delta \Pi =-\zeta \delta \Theta -\delta \left(\tau \right)\dot{\Pi }-\delta \left(\zeta \right)\Theta +\alpha \zeta \delta {\mathcal{K}}_G+\alpha \delta \zeta {\mathcal{K}}_G,$$

(45)

where all background functions are time-dependent and evaluated on FRW. The metric and tetrad perturbations also modify the geometrical sector of $f(T,T_G)$, producing a Poisson-like equation and a slip relation of the form

$$\begin{array}{cc}\hfill -{\displaystyle \frac{k^2}{a^2}}\Psi & =4\pi G\mu (a,k)\rho \Delta ,\hfill \end{array}$$

(46)

$$\begin{array}{cc}\hfill \Phi -\Psi & =\Sigma (a,k){\displaystyle \frac{8\pi G\rho \Delta }{k^2/a^2}},\hfill \end{array}$$

(47)

where $\Delta \equiv \delta +3aH(1+w)\theta /k^2$ is the comoving density contrast. The functions $\mu (a,k)$ and $\Sigma (a,k)$ encode departures from GR due to $f(T,T_G)$ and are algebraic combinations of $(f_T,f_{TT},f_{T_G},\dots )$; their explicit expressions are lengthy and relegated to Appendix A. The slip parameter is $\eta (a,k)\equiv \Phi /\Psi -1$ and relates to $\Sigma$ and $\mu$ by $\eta =\Sigma /\mu$ in the quasistatic limit (see below).

6.3. Quasistatic Subhorizon Limit

On scales $k\gg aH$, time derivatives of the potentials are subleading, and (42) reduces to $\delta \Theta \simeq -\theta$. Retaining the leading terms in $k^2/a^2$ and using (46)–(47), the continuity and Euler equations become

$$\begin{array}{cc}\hfill \dot{\delta }& \simeq -(1+w)\theta ,\hfill \end{array}$$

(48)

$$\begin{array}{cc}\hfill \dot{\theta }& \simeq -H(1-3w)\theta +{\displaystyle \frac{k^2}{a^2}}\Psi -{\displaystyle \frac{k^2}{a^2}}\left({\displaystyle \frac{c_s^2\delta \rho +\delta \Pi }{\rho +p}}\right).\hfill \end{array}$$

(49)

The linearized YKGC law (45) simplifies to

$$\tau \dot{\delta \Pi }+\delta \Pi \simeq \zeta \theta +\alpha \zeta \delta {\mathcal{K}}_G,$$

(50)

where $\delta {\mathcal{K}}_G$ depends on perturbations of $(f_T,f_{T_G})$ and is finite by construction (Section 3).

6.4. Growth Equation and Effective Drag

For late-time matter-dominated growth we set $w\simeq 0$ and $c_s^2\simeq 0$ for the equilibrium sector. Combining (48)–(50) and eliminating $\theta$ and $\delta \Pi$ yields a second-order growth equation for the matter contrast,

$${\delta }^{\prime \prime }+\left[2+{\displaystyle \frac{H^{\prime }}{H}}+{\Gamma }_{\mathrm{visc}}(a,k)\right]{\delta }^{\prime }-{\displaystyle \frac{3}{2}}\mu (a,k){\Omega }_m\left(a\right)\delta ={\mathcal{S}}_{\mathrm{visc}}(a,k)\delta ,$$

(51)

where prime is $d/dlna$, and ${\Gamma }_{\mathrm{visc}}$ is an effective drag coming from the viscous sector,

$${\Gamma }_{\mathrm{visc}}(a,k)\equiv {\displaystyle \frac{1}{aH\tau \left(H\right)}}{\displaystyle \frac{\Pi }{\rho }}{\mathcal{F}}_1(a,k;\alpha ,\zeta ,\tau ,{\mathcal{K}}_G),{\mathcal{S}}_{\mathrm{visc}}\equiv {\displaystyle \frac{k^2}{a^2{H}^2}}{\mathcal{F}}_2(a,k;\alpha ,\zeta ,\tau ,{\mathcal{K}}_G),$$

(52)

with ${\mathcal{F}}_{1,2}$ dimensionless, order-unity functions obtained by solving (50) (explicit formulae are given in Appendix A for our $f(T,T_G)$ ansatz). The YKGC priors ensure $\tau >0$ and keep $\Pi /(\rho +p)$ bounded, so ${\Gamma }_{\mathrm{visc}}\ge 0$ acts as a friction term that tends to dampen ${\delta }^{\prime }$; ${\mathcal{S}}_{\mathrm{visc}}$ represents a small scale-dependent correction that vanishes when $\alpha \to 0$. The baseline GR growth equation can be found in [40]

Observable Dictionary

At the linear level we track three knobs: (i) a geometric modification $\mu (a,k)$ in the Poisson sector, (ii) a gravitational slip $\eta (a,k)$, and (iii) an effective viscous drag ${\Gamma }_{\mathrm{visc}}(a,k)$ from $\Pi$-dynamics. Schematically, the growth equation acquires a friction shift and an effective source,

$${\delta }^{\prime \prime }+\left(2+{\mathcal{H}}^{\prime }/\mathcal{H}+{\Gamma }_{\mathrm{visc}}\right){\delta }^{\prime }-{\displaystyle \frac{3}{2}}\mu (a,k){\Omega }_m\left(a\right)\delta =0,$$

so that $f{\sigma }_8\left(z\right)$ is damped primarily by ${\Gamma }_{\mathrm{visc}}$ while lensing and the $E_G$ statistic respond to the combination $\mu (1+\eta )$. In YKGC, $(\alpha ,{\lambda }_{1,2})$ dominantly feed $\mu$ and $\eta$ through ${\mathcal{K}}_G$, whereas $({\zeta }_i,{\tau }_i)$ dominantly control ${\Gamma }_{\mathrm{visc}}$ via $c_{\mathrm{visc}}^2=\zeta /\left[\tau (\rho +p)\right]$.

6.5. Effective Newton Constant and Slip

We define the effective Newton coupling and slip as

$$G_{\mathrm{eff}}(a,k)\equiv G\mu (a,k),\eta (a,k)\equiv {\displaystyle \frac{\Phi }{\Psi }}-1.$$

(53)

In our $f(T,T_G)$ reconstruction (18), $\mu (a,k)$ and $\eta (a,k)$ are background-determined algebraic functions that reduce to $(1,0)$ in the GR limit. The viscous sector modifies growth only via ${\Gamma }_{\mathrm{visc}}$ and ${\mathcal{S}}_{\mathrm{visc}}$, not through a direct change of the Poisson kernel, so background–perturbation degeneracy is broken.

6.6. Observables: $f{\sigma }_8\left(z\right)$, ISW and $E_G$

The linear growth rate is $f\equiv dlnD/dlna$, where $D\left(a\right)$ solves (51) with $D(a\to 0)=a$. The RSD observable is

$$f{\sigma }_8\left(z\right)=f\left(z\right)D\left(z\right){\sigma }_{8,\mathrm{ini}},$$

(54)

with ${\sigma }_{8,\mathrm{ini}}$ fixed by the chosen normalization at $z=0$. The late-time Integrated Sachs–Wolfe (ISW) signal probes $\dot{\Phi }+\dot{\Psi }$, which depends on $(\mu ,\eta )$ and on viscous drag through the time evolution of $\delta$. A scale-robust discriminator is the $E_G$ statistic,

$$E_G(k,z)\equiv {\displaystyle \frac{k^2(\Phi +\Psi )}{3H_0^2{a}^{-1}f\delta }}\simeq {\displaystyle \frac{\mu (a,k)[1+\eta (a,k\left)\right]}{2f(a,k)}},$$

(55)

in the quasistatic regime. Our framework predicts $E_G$ shifts driven by $\mu ,\eta$ and a suppressed f via ${\Gamma }_{\mathrm{visc}}$, providing a multi-probe handle beyond background distances.

The model parameters $\alpha$ and ${\lambda }_{1,2}$ determine the coupling between the teleparallel geometric kernel $K_G(T,T_G)$ and the effective fluid dynamics, which in turn directly modifies the effective Newtonian coupling $\mu (k,a)$ and the gravitational slip parameter $\eta (k,a)$. For representative values $\alpha =0.1$ and ${\lambda }_1={\lambda }_2=1$, the predicted value of the lensing–growth consistency statistic $E_G\left(z\right)$ deviates from the standard GR expectation by approximately $\mathcal{O}(1\%)$ at intermediate redshifts ($z\sim 0.5-1$) [54]. A deviation of this magnitude lies within the forecast sensitivity of upcoming weak lensing surveys such as Euclid and LSST. The bulk viscosity parameters ${\zeta }_{0,1}$ and ${\tau }_{0,1}$ determine the effective viscous damping rate ${\Gamma }_{\mathrm{visc}}$ in the linear growth equation. For illustrative values ${\zeta }_0={10}^{-3}$, ${\zeta }_1={10}^{-2}$, ${\tau }_0=0.1H_0^{-1}$, and ${\tau }_1=0.05H_0^{-1}$, the resulting suppression of the growth function $f{\sigma }_8\left(z\right)$ reaches approximately $2\%$ at $z\approx 0.8$ relative to the Planck 2018 $\Lambda$CDM baseline [1,52,53]. This damping amplitude is comparable to the precision of current redshift-space distortion measurements (e.g., BOSS, eBOSS) and provides a concrete observational signature of the viscous closure mechanism. See Table 1 for a summary of the qualitative mapping between model parameters and observational quantities.

6.7. Initial Conditions and Normalization

For CMB-safe initial conditions we impose a GR-like regime at early times: $\mu \to 1$, $\eta \to 0$, $\Pi /(\rho +p)\to 0$, consistent with the priors of Section 4. We initialize $D\left(a\right)$ deep in matter domination and match to the inflationary exit via $H\left(N\right)$ continuity. The tensor speed constraint is enforced by our parameter priors (GW-compatible sector).

6.8. Stability and Absence of Pathologies

No ghostlike or gradient instabilities arise at the scalar level provided $c_{\mathrm{visc}}^2\le 1$ (Section 3) and the algebraic kernels $(\mu ,\eta )$ remain finite for finite H; both hold under our regulated kernel and early-time priors. Tensor modes keep luminal speed in the admissible parameter region by construction.

7. YKGC for Astrophysical Systems: Standard Form and Adapters

7.1. Canonical Standard Form

We promote the closure to a portable standard:

$$\tau u^{\mu }{\nabla }_{\mu }\Pi +\Pi =-\zeta \Theta +\alpha \zeta \mathcal{K},\Theta \equiv {\nabla }_{\mu }{u}^{\mu },$$

(56)

with $\zeta >0$, $\tau >0$, and a dimensionless kernel $\mathcal{K}$. In teleparallel cosmology we set $\mathcal{K}={\mathcal{K}}_G(T,T_G)$ (Section 3). For other contexts:

Metric GR Adapter

When torsion invariants are absent, define

$${\mathcal{S}}_{\mathrm{GR}}\equiv 1+{\lambda }_R{L}^2\left|R\right|+{\lambda }_G{L}^4\left|G\right|,\mathcal{K}\equiv u^{\mu }{\nabla }_{\mu }ln{\mathcal{S}}_{\mathrm{GR}},$$

(57)

with Ricci scalar R, Gauss–Bonnet G, a characteristic scale L, and dimensionless ${\lambda }_{R,G}$ [41].

Newtonian/MHD Adapter

For compressible flows with velocity $v^i$, $\theta \equiv {\nabla }_i{v}^i$, shear ${\sigma }_{ij}$ and vorticity ${\omega }_{ij}$, set [42,43]

$$\begin{array}{cc}\hfill {\mathcal{S}}_{\mathrm{N}}& \equiv 1+{\lambda }_{\theta }{\displaystyle \frac{\left|\theta \right|}{{\theta }_{★}}}+{\lambda }_{\sigma }{\displaystyle \frac{\sqrt{2{\sigma }_{ij}{\sigma }^{ij}}}{{\sigma }_{★}}}+{\lambda }_{\omega }{\displaystyle \frac{\sqrt{2{\omega }_{ij}{\omega }^{ij}}}{{\omega }_{★}}},\hfill \end{array}$$

(58)

$$\begin{array}{cc}\hfill \mathcal{K}& \equiv D_{t}ln{\mathcal{S}}_{\mathrm{N}},D_t\equiv {\partial }_t+v^i{\nabla }_i,\hfill \end{array}$$

(59)

so that the nonrelativistic YKGC law reads

$$\tau D_t\Pi +\Pi =-\zeta \theta +\alpha \zeta D_{t}ln{\mathcal{S}}_{\mathrm{N}}.$$

(60)

7.2. Minimal Parameterization for Adopters

We recommend

$${\Theta }_{\mathrm{YKGC}}^{\min }=\{{\zeta }_0,{\zeta }_1,{\tau }_0,{\tau }_1,\alpha ,{\lambda }_{•}\},$$

with $\zeta ={\zeta }_0+{\zeta }_1\mathcal{H}$ and $\tau ={\tau }_0+{\tau }_1/\mathcal{H}$, where $\mathcal{H}=H$ in cosmology and $\mathcal{H}=\left|\theta \right|$ in NR flows. Thermodynamic and causal priors from Section 3 apply unchanged.

8. Well-Posedness and Entropy Production: A Sufficient Theorem

Theorem 1 

(YKGC sufficiency). Consider (56) with $\zeta >0$, $\tau >0$, a dimensionless kernel $\mathcal{K}$ defined as a convective logarithmic derivative of a positive scalar (e.g. (8), (57), (59)). Assume a barotropic background with $T_{\mathrm{loc}}>0$ and smooth coefficients. If the bounds

$$\left|\alpha \mathcal{K}\right|\le {\displaystyle \frac{\left|\Pi \right|}{\zeta }},c_{\mathrm{visc}}^2\equiv {\displaystyle \frac{\zeta }{\tau (\rho +p)}}\le 1,$$

(61)

hold locally, then (i) the entropy production density satisfies $T_{\mathrm{loc}}{\nabla }_{\mu }{s}^{\mu }\ge 0$ and (ii) the Cauchy problem for $(\rho ,\Pi )$ is locally well posed with subluminal characteristics for the viscous sector  [44].

Proof 

(Sketch). Using $T_{\mathrm{loc}}{\nabla }_{\mu }{s}^{\mu }={\Pi }^2/\zeta -\alpha \Pi \mathcal{K}$ and the bound in (61), we have non-negative production. Linearizing (56) yields a first-order hyperbolic equation for $\delta \Pi$ with characteristic speed $c_{\mathrm{visc}}$; the second bound confines the cone within the light cone. Smooth coefficients ensure local existence/uniqueness.    □

Remark. Equality in (61) saturates irreversible production. Any stricter bound $\left|\alpha \mathcal{K}\right|\le \chi \left|\Pi \right|/\zeta$ with $0<\chi \le 1$ yields a Lyapunov-type decay of $\Pi$ toward the steady solution.

9. Conclusions

We have formulated a variationally closed, causal bulk–viscous law (YKGC) for a single cosmological fluid within teleparallel Gauss–Bonnet gravity. The closure follows from a local Rayleigh–Onsager functional with a dimensionless teleparallel kernel ${\mathcal{K}}_G(T,T_G)$ and yields the hyperbolic relaxation equation (9). We established two sufficient conditions—non–negative entropy production and subluminal characteristics—which we adopt as hard priors throughout the analysis. The construction reduces continuously to the truncated Israel–Stewart theory when the kernel decouples and to GR for linear f, with early–time consistency imposed.

On the background we reconstructed $f(T,T_G)$ using a two–plateau $H\left(N\right)$ profile (16), obtaining unified quasi–de Sitter inflation and late–time acceleration while recovering the GR limit at early times. The dynamical–systems analysis identifies a stable late–time de Sitter attractor and provides explicit criteria that avoid finite–time singularities under the same thermodynamic and causal bounds. The regularized kernel remains well defined around $H\to 0$, allowing nonsingular extensions when desired.

At the linear level we derived a growth equation in the quasistatic regime that separates geometric effects—encoded in an effective Newton coupling $\mu (a,k)$—from a purely viscous drag${\Gamma }_{\mathrm{visc}}$ induced by the YKGC closure. Within thermodynamically allowed priors this structure generically produces a controlled damping of $f{\sigma }_8\left(z\right)$ and correlated shifts in the $E_G$ statistic, supplying concrete observational targets beyond background distances.

The framework is compact in parameters, explicit in its limits, and consistent with early–time bounds and gravitational–wave speed constraints via priors set in Section 4. Future investigations can build on the present framework in several well-defined directions. A natural next step is to confront the model with joint SN, BAO, $H\left(z\right)$ and RSD datasets, complemented by CMB consistency checks including ISW and lensing, in order to place quantitative constraints on the viscous parameters. Extending the current linear treatment to the nonlinear regime and exploring tensor-mode signatures would help to assess the model’s phenomenology on smaller scales and its potential implications for gravitational wave propagation. Finally, the local and causal structure of the YKGC law makes it readily adaptable as a bulk-viscous module in other cosmological or astrophysical contexts, providing a clear path for future applications. These developments lie beyond the scope of the present work but offer concrete opportunities for further testing and extending the framework introduced here.

Summary

We proposed a variationally closed, causal bulk–viscous law (YKGC) for a single fluid evolving in teleparallel $f(T,T_G)$ gravity, proved two sufficient conditions that guarantee non–negative entropy production and subluminal characteristics, and showed how the same framework unifies a quasi–de Sitter inflationary phase with late–time acceleration while recovering GR at early times. On the linear side we derived $G_{\mathrm{eff}}(k,a)$ and the metric slip, finding a controlled, scale–dependent suppression of growth that leaves background distances close to $\Lambda$CDM under the priors of Sec. 4. These elements define a clean, testable target for late–time data combinations.

Testable Predictions

Within the thermodynamic and causality priors (entropy bound (11) and subluminality bound (13)), the framework predicts: (i) a late–time damping of $f{\sigma }_8\left(z\right)$ [45] governed mainly by the viscous sector $({\zeta }_i,{\tau }_i)$ via ${\Gamma }_{\mathrm{visc}}$, (ii) correlated shifts of the scale–robust statistic $E_G$ [45] driven chiefly by $(\alpha ,{\lambda }_{1,2})$ through the teleparallel kernel, and (iii) background distances that can remain $\Lambda$CDM–like when early–time and GW–speed consistency are imposed. As emphasized in the Introduction and Sec. 6, these signatures point to joint constraints from SN+BAO+$H\left(z\right)$ and RSD/lensing. We do not attempt to address the $H_0$ tension here; our baseline keeps background distances anchored while targeting the late–time growth tension.

Outlook

(1) Joint late–time analysis: fit $\{\alpha ,{\lambda }_{1,2},{\zeta }_0,{\zeta }_1,{\tau }_0,{\tau }_1\}$ under the priors of Sec. 4 to SN+BAO+$H\left(z\right)$+RSD and weak–lensing datasets, reporting posteriors for $(\mu ,\eta ,{\Gamma }_{\mathrm{visc}})$ along with $f{\sigma }_8$ and $E_G$. (2) Robustness: stress–test the priors by varying the $f(T,T_G)$ ansatz within the GR– and IS–consistent limits to quantify model dependence. (3) Cross–domain adopters: use the “canonical standard form” (Sec. 7) to port YKGC to metric GR (adapter in Eq. (57)) and to Newtonian/MHD flows (Eqs. (58)–(60)), enabling astrophysical tests with the same closure.

Appendix A.1. Building Blocks and Notation

On FRW, $T=-6H^2$ and $T_G=24H^2(\dot{H}+H^2)$ [13]. Define the teleparallel derivatives

$$f_T\equiv {\displaystyle \frac{\partial f}{\partial T}},f_{T_G}\equiv {\displaystyle \frac{\partial f}{\partial T_G}},f_{TT}\equiv {\displaystyle \frac{{\partial }^{2}f}{\partial T^2}},f_{T{T}_G}\equiv {\displaystyle \frac{{\partial }^{2}f}{\partial T\partial T_G}},f_{T_G{T}_G}\equiv {\displaystyle \frac{{\partial }^{2}f}{\partial T_G^2}}.$$

We use Newtonian gauge with $d{s}^2=-(1+2\Phi )d{t}^2+a^2(1-2\Psi )d{\overrightarrow{x}}^2$ and scalar velocity divergence $\theta \equiv k^{2}v/a^2$. The expansion scalar $\Theta \equiv {\nabla }_{\mu }{u}^{\mu }=3H$ perturbs as

$$\delta \Theta =-{\displaystyle \frac{k^2}{a^2}}v+3(\dot{\Psi }+H\Phi ),\delta H={\displaystyle \frac{1}{3}}\delta \Theta .$$

(A1)

Appendix A.2. Teleparallel Variations δT, δT G and δK G

Since $T=-6H^2$,

$$\delta T={\displaystyle \frac{\partial T}{\partial H}}\delta H=-12H\delta H=-4H\delta \Theta .$$

(A2)

For $T_G=24H^2(\dot{H}+H^2)\equiv 24H^{2}A$ with $A\equiv \dot{H}+H^2$,

$$\delta T_G=48HA\delta H+24H^2\delta A,\delta A=\delta \dot{H}+2H\delta H.$$

(A3)

Using (A1), $\delta \dot{H}=\frac{1}{3}\frac{d}{dt}\left(\delta \Theta \right)$ up to subleading lapse terms. Hence

$$\delta T_G=16HA\delta \Theta +8H^2{\displaystyle \frac{d}{dt}}\left(\delta \Theta \right).$$

(A4)

In the quasistatic (QS) subhorizon regime ($k\gg aH$), $\delta \Theta \simeq -\theta$, so

$$\delta T\simeq 4H\theta ,\delta T_G\simeq -16H(\dot{H}+H^2)\theta -8H^2\dot{\theta }.$$

(A5)

The regulated geometric kernel is ${\mathcal{K}}_G=u^{\mu }{\nabla }_{\mu }ln{\mathcal{S}}_G$ with

$${\mathcal{S}}_G\equiv 1+{\lambda }_1\sqrt{f_T^2+{\epsilon }^2}+{\lambda }_2\sqrt{{\left({\displaystyle \frac{T_G}{T^2+\epsilon H_{★}^4}}\right)}^2{f}_{T_G}^2+{\epsilon }^2},Q\equiv {\displaystyle \frac{T_G}{T^2+\epsilon H_{★}^4}}.$$

(A6)

Its background time derivative is

$${\mathcal{K}}_G={\displaystyle \frac{{\dot{\mathcal{S}}}_G}{{\mathcal{S}}_G}}={\displaystyle \frac{1}{{\mathcal{S}}_G}}\left({\partial }_T{\mathcal{S}}_G\dot{T}+{\partial }_{T_G}{\mathcal{S}}_G{\dot{T}}_G\right),\dot{T}=-12H\dot{H},{\dot{T}}_G=48H\dot{H}(\dot{H}+H^2)+24H^2(\ddot{H}+2H\dot{H}).$$

(A7)

Linearizing,

$$\delta {\mathcal{K}}_G=u^{\mu }{\nabla }_{\mu }\left({\displaystyle \frac{\delta {\mathcal{S}}_G}{{\mathcal{S}}_G}}\right)+\delta u^{\mu }{\nabla }_{\mu }ln{\mathcal{S}}_G,{\displaystyle \frac{\delta {\mathcal{S}}_G}{{\mathcal{S}}_G}}={\displaystyle \frac{{\partial }_T{\mathcal{S}}_G}{{\mathcal{S}}_G}}\delta T+{\displaystyle \frac{{\partial }_{T_G}{\mathcal{S}}_G}{{\mathcal{S}}_G}}\delta T_G.$$

(A8)

In QS, $\delta u^0\simeq -\Phi$ and $\delta {\mathcal{K}}_G\simeq {\displaystyle \frac{d}{dt}}(\delta {\mathcal{S}}_G/{\mathcal{S}}_G)-\Phi {\mathcal{K}}_G$, with $\delta T,\delta T_G$ from (A5).

Appendix A.3. Field Equations in Algebraic QS Form

At first order the modified teleparallel field equations yield a Poisson-like constraint and a slip relation which, in QS, can be cast as

$$\begin{array}{cc}\hfill -{\displaystyle \frac{k^2}{a^2}}\Psi & =4\pi G\mu (a,k)\rho \Delta ,\Delta \equiv \delta +{\displaystyle \frac{3aH(1+w)}{k^2}}\theta ,\hfill \end{array}$$

(A9)

$$\begin{array}{cc}\hfill \Phi -\Psi & ={\displaystyle \frac{\Sigma (a,k)}{8\pi G\rho \Delta }}{\displaystyle \frac{k^2}{a^2}}.\hfill \end{array}$$

(A10)

Solving the linearized tetrad equations in QS gives $\mu$ and $\Sigma$ as rational algebraic functions of f-derivatives and background H [47],

$$\begin{array}{|c|}\hline \mu (a,k)={\displaystyle \frac{{\mathcal{N}}_{\mu }}{\mathcal{D}}},\Sigma (a,k)={\displaystyle \frac{{\mathcal{N}}_{\Sigma }}{\mathcal{D}}},\eta (a,k)\equiv {\displaystyle \frac{\Phi }{\Psi }}-1={\displaystyle \frac{\Sigma }{\mu }}-1\\ \hline\end{array}.$$

(A11)

Here the common denominator and numerators read

$$\begin{array}{cc}\hfill \mathcal{D}& ={\mathcal{E}}_1+{\mathcal{E}}_2{\displaystyle \frac{k^2}{a^2}},\hfill \end{array}$$

(A12)

$$\begin{array}{cc}\hfill {\mathcal{N}}_{\mu }& ={\mathcal{E}}_3+{\mathcal{E}}_4{\displaystyle \frac{k^2}{a^2}},{\mathcal{N}}_{\Sigma }={\mathcal{E}}_5+{\mathcal{E}}_6{\displaystyle \frac{k^2}{a^2}},\hfill \end{array}$$

(A13)

with QS coefficients

$$\begin{array}{cc}\hfill {\mathcal{E}}_1& \equiv f_T+12H^2{f}_{TT}+24H(\dot{H}+H^2)f_{T{T}_G}+24H^2{f}_{T_G}+48H^4{f}_{T_G{T}_G},\hfill \end{array}$$

(A14)

$$\begin{array}{cc}\hfill {\mathcal{E}}_2& \equiv {\displaystyle \frac{1}{H^2}}\left[c_1{f}_{TT}+c_2{f}_{T{T}_G}+c_3{f}_{T_G{T}_G}\right],\hfill \end{array}$$

(A15)

$$\begin{array}{cc}\hfill {\mathcal{E}}_3& \equiv f_T+d_1{H}^2{f}_{TT}+d_{2}H(\dot{H}+H^2)f_{T{T}_G}+d_3{H}^2{f}_{T_G}+d_4{H}^4{f}_{T_G{T}_G},\hfill \end{array}$$

(A16)

$$\begin{array}{cc}\hfill {\mathcal{E}}_4& \equiv {\displaystyle \frac{1}{H^2}}\left[e_1{f}_{TT}+e_2{f}_{T{T}_G}+e_3{f}_{T_G{T}_G}\right],\hfill \end{array}$$

(A17)

$$\begin{array}{cc}\hfill {\mathcal{E}}_5& \equiv s_1{H}^2{f}_{TT}+s_{2}H(\dot{H}+H^2)f_{T{T}_G}+s_3{H}^2{f}_{T_G}+s_4{H}^4{f}_{T_G{T}_G},\hfill \end{array}$$

(A18)

$$\begin{array}{cc}\hfill {\mathcal{E}}_6& \equiv {\displaystyle \frac{1}{H^2}}\left[u_1{f}_{TT}+u_2{f}_{T{T}_G}+u_3{f}_{T_G{T}_G}\right].\hfill \end{array}$$

(A19)

The dimensionless numbers $\{c_i,d_i,e_i,s_i,u_i\}$ are order-unity QS coefficients fixed by the linearized tetrad equations (their explicit forms follow from your $f(T,T_G)$ ansatz and the chosen covariant teleparallel prescription; we provide them in the reference implementation and export them as a kernels.tex file to be \input’ed here). In the GR limit $f_T\to 1$ and all higher derivatives $\to 0$, so that

$$\mu \to 1,\Sigma \to 0,\eta \to 0,$$

(A20)

as required.

Appendix A.4. Growth Equation and Viscous Drag (For Completeness)

Combining the conservation equations with the linearized YKGC law $\tau \dot{\delta \Pi }+\delta \Pi \simeq \zeta \theta +\alpha \zeta \delta {\mathcal{K}}_G$ one obtains

$${\delta }^{\prime \prime }+\left[2+{\displaystyle \frac{H^{\prime }}{H}}+{\Gamma }_{\mathrm{visc}}(a,k)\right]{\delta }^{\prime }-{\displaystyle \frac{3}{2}}\mu (a,k){\Omega }_m\left(a\right)\delta =S_{\mathrm{visc}}(a,k)\delta ,$$

(A21)

with primes $d/dlna$. The effective drag and source are

$${\Gamma }_{\mathrm{visc}}(a,k)={\displaystyle \frac{1}{aH\tau \left(H\right)}}{\displaystyle \frac{\Pi }{\rho }}{F}_1(a,k;\alpha ,\zeta ,\tau ,{\mathcal{K}}_G),S_{\mathrm{visc}}(a,k)={\displaystyle \frac{k^2}{a^2{H}^2}}{F}_2(a,k;\alpha ,\zeta ,\tau ,{\mathcal{K}}_G),$$

(A22)

where $F_{1,2}$ are dimensionless $\mathcal{O}\left(1\right)$ algebraic functions determined by (A8) and the background; they vanish when $\alpha \to 0$ (decoupling).

Appendix A.5. Consistency with Thermodynamic and Causal Priors

All occurrences of ${\mathcal{K}}_G$ and $\delta {\mathcal{K}}_G$ in the YKGC law are multiplied by $\alpha \zeta \left(H\right)$. We impose, at each step,

$$\left|\alpha {\mathcal{K}}_G\right|\le {\displaystyle \frac{\left|\Pi \right|}{\zeta \left(H\right)}},c_{\mathrm{visc}}^2={\displaystyle \frac{\zeta \left(H\right)}{\tau \left(H\right)[\rho +p]}}\le 1,$$

(A23)

ensuring non–negative entropy production and subluminal characteristics.

Appendix B.1. Background Geometric Sector on FRW

We adopt the action (Sec. Section 2) with $f(T,T_G)$ given by (18). On FRW one has $T=-6H^2$ and $T_G=24H^2(\dot{H}+H^2)$, hence

$$f_T\equiv {\displaystyle \frac{\partial f}{\partial T}},f_{T_G}\equiv {\displaystyle \frac{\partial f}{\partial T_G}},{\dot{f}}_T=f_{TT}\dot{T}+f_{T{T}_G}{\dot{T}}_G,{\dot{f}}_{T_G}=f_{T{T}_G}\dot{T}+f_{T_G{T}_G}{\dot{T}}_G.$$

The effective fluid $({\rho }_f,p_f)$ entering (5)–(6) is obtained by inserting $T\left(T_G\right)$ into the field equations and moving all non-Einstein terms to the right-hand side. In practice: (i) compute $\{f,f_T,f_{T_G},{\dot{f}}_T,{\dot{f}}_{T_G}\}$; (ii) substitute into the background equations; (iii) solve algebraically for ${\rho }_f$ in (5) and $p_f$ in (6). The resulting expressions are lengthy but uniquely fixed by (18) and (17). We will list these explicitly in the code release associated with this paper.

Appendix B.2. Regulated Geometric Kernel K G

Let ${\mathcal{S}}_G$ denote the positive regulator argument in (8):

$${\mathcal{S}}_G\equiv 1+{\lambda }_1\sqrt{f_T^2+{\epsilon }^2}+{\lambda }_2\sqrt{{\underbrace{{\left({\displaystyle \frac{T_G}{T^2+\epsilon H_{★}^4}}\right)}^2{f}_{T_G}^2+{\epsilon }^2}}_{{\displaystyle \equiv \mathcal{Z}}}},Q\equiv {\displaystyle \frac{T_G}{T^2+\epsilon H_{★}^4}}.$$

(A24)

Then ${\mathcal{K}}_G=u^{\mu }{\nabla }_{\mu }ln{\mathcal{S}}_G$ and, for any time function $X\left(t\right)$ on FRW, $u^{\mu }{\nabla }_{\mu }X=\dot{X}$. Hence

$${\mathcal{K}}_G={\displaystyle \frac{{\dot{\mathcal{S}}}_G}{{\mathcal{S}}_G}}={\displaystyle \frac{1}{{\mathcal{S}}_G}}\left({\partial }_T{\mathcal{S}}_G\dot{T}+{\partial }_{T_G}{\mathcal{S}}_G{\dot{T}}_G\right),\dot{T}=-12H\dot{H},{\dot{T}}_G=48H\dot{H}(\dot{H}+H^2)+24H^2(\ddot{H}+2H\dot{H}).$$

(A25)

Explicit Partial Derivatives

Define $A\equiv \sqrt{f_T^2+{\epsilon }^2}$ and $B\equiv \sqrt{\mathcal{Z}}$. Using

$${\partial }_{T}A={\displaystyle \frac{f_T{f}_{TT}}{A}},{\partial }_{T_G}A={\displaystyle \frac{f_T{f}_{T{T}_G}}{A}},$$

and with $\mathcal{Z}=Q^2{f}_{T_G}^2+{\epsilon }^2$ and

$${\partial }_{T}Q={\displaystyle \frac{\partial }{\partial T}}\left({\displaystyle \frac{T_G}{T^2+\epsilon H_{★}^4}}\right)=-{\displaystyle \frac{2T{T}_G}{{(T^2+\epsilon H_{★}^4)}^2}},{\partial }_{T_G}Q={\displaystyle \frac{1}{T^2+\epsilon H_{★}^4}},$$

we obtain

$$\begin{array}{cc}\hfill {\partial }_T\mathcal{Z}& =2Q\left({\partial }_{T}Q\right)f_{T_G}^2+2Q^2{f}_{T_G}{f}_{T{T}_G},\hfill \end{array}$$

(A26)

$$\begin{array}{cc}\hfill {\partial }_{T_G}\mathcal{Z}& =2Q\left({\partial }_{T_G}Q\right)f_{T_G}^2+2Q^2{f}_{T_G}{f}_{T_G{T}_G}.\hfill \end{array}$$

(A27)

Therefore

$$\begin{array}{cc}\hfill {\partial }_T{\mathcal{S}}_G& ={\lambda }_1{\displaystyle \frac{f_T{f}_{TT}}{A}}+{\lambda }_2{\displaystyle \frac{1}{2B}}{\partial }_T\mathcal{Z},\hfill \end{array}$$

(A28)

$$\begin{array}{cc}\hfill {\partial }_{T_G}{\mathcal{S}}_G& ={\lambda }_1{\displaystyle \frac{f_T{f}_{T{T}_G}}{A}}+{\lambda }_2{\displaystyle \frac{1}{2B}}{\partial }_{T_G}\mathcal{Z}.\hfill \end{array}$$

(A29)

Plugging (A28)–(A29) into (A25) yields the fully explicit background kernel ${\mathcal{K}}_G$.

Appendix B.3. Linear Perturbation of K G

To linear order,

$$\delta {\mathcal{K}}_G=u^{\mu }{\nabla }_{\mu }\left({\displaystyle \frac{\delta {\mathcal{S}}_G}{{\mathcal{S}}_G}}\right)+\delta u^{\mu }{\nabla }_{\mu }ln{\mathcal{S}}_G.$$

(A30)

In Newtonian gauge, $d{s}^2=-(1+2\Phi )d{t}^2+a^2(1-2\Psi )d{\overrightarrow{x}}^2$, keeping leading quasi-static terms ($k\gg aH$), one has $\delta u^0\simeq -\Phi$ and spatial derivatives dominate time derivatives of the potentials; then

$$\delta {\mathcal{K}}_G\simeq {\displaystyle \frac{d}{dt}}\left({\displaystyle \frac{\delta {\mathcal{S}}_G}{{\mathcal{S}}_G}}\right)-\Phi {\mathcal{K}}_G,{\displaystyle \frac{\delta {\mathcal{S}}_G}{{\mathcal{S}}_G}}={\displaystyle \frac{{\partial }_T{\mathcal{S}}_G}{{\mathcal{S}}_G}}\delta T+{\displaystyle \frac{{\partial }_{T_G}{\mathcal{S}}_G}{{\mathcal{S}}_G}}\delta T_G.$$

(A31)

QS Expressions for δT and δT G

The expansion scalar is $\Theta \equiv {\nabla }_{\mu }{u}^{\mu }=3H$. At first order

$$\delta \Theta =-{\displaystyle \frac{k^2}{a^2}}v+3(\dot{\Psi }+H\Phi ),\Rightarrow \delta H={\displaystyle \frac{1}{3}}\delta \Theta .$$

(A32)

Since $T=-6H^2$, one finds

$$\delta T={\displaystyle \frac{\partial T}{\partial H}}\delta H=-12H\delta H=-4H\delta \Theta .$$

(A33)

In the QS subhorizon limit ($k\gg aH$), $\delta \Theta \simeq -\theta$ with $\theta \equiv k^{2}v/a^2$, hence

$$\delta T\simeq 4H\theta .$$

(A34)

For $T_G=24H^2(\dot{H}+H^2)\equiv 24H^{2}A$ with $A\equiv \dot{H}+H^2$,

$$\delta T_G=48HA\delta H+24H^2\delta A,\delta A=\delta \dot{H}+2H\delta H.$$

(A35)

Using (42), $\delta \dot{H}=\frac{1}{3}\frac{d}{dt}\left(\delta \Theta \right)$ up to lapse corrections that are subleading in the QS limit. Therefore,

$$\delta T_G\simeq 16HA\delta \Theta +8H^2{\displaystyle \frac{d}{dt}}\left(\delta \Theta \right)\underset{\mathrm{QS}}{\overset{\delta \Theta \simeq -\theta }{\to }}\delta T_G\simeq -16HA\theta -8H^2\dot{\theta }.$$

(A36)

Equations (A34) and (A36), together with (A28)–(A29), give a closed algebraic expression for $\delta {\mathcal{K}}_G$ via (A31). In our numerical implementation we use (A30) directly, evaluating all background derivatives from time series and computing $(\delta T,\delta T_G)$ with (A33)–(A35).

Remark (Consistency with Priors) All occurrences of ${\mathcal{K}}_G$ and $\delta {\mathcal{K}}_G$ enter multiplied by $\alpha \zeta \left(H\right)$ in the YKGC law, and the thermodynamic and causal priors,

$$\left|\alpha {\mathcal{K}}_G\right|\le {\displaystyle \frac{\left|\Pi \right|}{\zeta \left(H\right)}},c_{\mathrm{visc}}^2={\displaystyle \frac{\zeta \left(H\right)}{\tau \left(H\right)(\rho +p)}}\le 1,$$

are imposed at each step of the reconstruction/evolution.

Appendix C.1. Entropy Positivity (Sufficient Condition)

Starting from (10) with $T_{\mathrm{loc}}>0$ and $\zeta >0$,

$$T_{\mathrm{loc}}{\nabla }_{\mu }{s}^{\mu }={\displaystyle \frac{{\Pi }^2}{\zeta }}-\alpha \Pi {\mathcal{K}}_G\ge {\displaystyle \frac{{\Pi }^2}{\zeta }}-\left|\alpha \right|\left|\Pi \right||{\mathcal{K}}_G|.$$

Thus $T_{\mathrm{loc}}{\nabla }_{\mu }{s}^{\mu }\ge 0$ is ensured by

$$|\alpha {\mathcal{K}}_G|\le {\displaystyle \frac{\left|\Pi \right|}{\zeta }},$$

which is (11). Equality saturates irreversible production. Any stricter bound $|\alpha {\mathcal{K}}_G|\le \chi |\Pi |/\zeta$ with $0<\chi \le 1$ implies a Lyapunov-type decay of $\Pi$ toward the steady solution because $\dot{\Pi }$ obeys (9) with a strictly dissipative linear term.

Appendix C.2. Causality and Hyperbolicity

Linearizing (9) about a homogeneous background and Fourier-analyzing a plane-wave perturbation along the fluid worldline yields

$$\tau {\partial }_t\delta \Pi +\delta \Pi =-\zeta \delta \Theta +\alpha \zeta \delta {\mathcal{K}}_G.$$

The principal part in $(\delta \Pi ,\delta \Theta )$ forms a first-order hyperbolic system with characteristic speed $c_{\mathrm{visc}}^2=\zeta /\left[\tau (\rho +p)\right]$. Imposing $0<c_{\mathrm{visc}}^2\le 1$ (i.e. $\tau >0$, $\zeta >0$ and (13)) confines the characteristic cone within the light cone, ensuring local well-posedness (existence/uniqueness) by standard theorems for linear symmetric-hyperbolic systems with smooth coefficients.

Appendix D.1. IMEX Update for the YKGC Law

For stiff regimes we adopt the IMEX update

$${\Pi }^{n+1}={\displaystyle \frac{{\Pi }^n+\Delta t[-{\zeta }^n{\Theta }^n+\alpha {\zeta }^n{\mathcal{K}}_G^n]}{1+\Delta t/{\tau }^{n+1}}},$$

(A37)

with ${\tau }^{n+1}$ evaluated from the updated primitives and $(\zeta ,\Theta ,{\mathcal{K}}_G)$ at level n (or $n+\frac{1}{2}$ by predictor)[48,49]. This preserves hyperbolicity and avoids overshoot when $\Delta t\gtrsim \tau$.

Appendix D.2. Benchmark Set

• NR acoustic decay: analytic attenuation rate vs. YKGC-Lite (60).
• Shock tube with bulk damping (NR/MHD): post-shock $\Pi$ and density profiles.
• FRW recovery: $H\left(N\right)$ and $\Pi \left(N\right)$ vs. quasi-steady IS limit.
• Linear growth: $f{\sigma }_8\left(z\right)$ damping and $E_G$ shifts from (51).

Appendix D.3. Reproducibility

Parameter files and scripts reproduce Figs. 1–5. The repository includes adapters (57), (59) and utilities to compute $\delta {\mathcal{K}}_G$ per (A30).