Dynamical System and Statefinder Analysis of Cosmological Models in f(T, B) Gravity

A major puzzle in modern physics is the universe's late-time acceleration [1−3] and the observed mass discrepancy in large-scale structures [4], whose dynamics cannot be explained by visible matter [5]. These phenomena are often attributed to dark energy and dark matter – hypothetical components that influence cosmic evolution despite having no direct electromagnetic signatures.

Dark energy presents a particular enigma due to its repulsive gravitational effect, which counteracts conventional gravitational attraction. Two major theoretical approaches have been developed to address these cosmological puzzles. The first retains the geometric framework of General Relativity (GR) while introducing new matter components, such as scalar fields with negative pressure. The second approach fundamentally modifies gravitational theory itself, either by extending the geometry underlying Einstein's field equations or by reinterpreting how matter influences spacetime curvature. Such modified gravity theories generally generalize the Einstein-Hilbert action via various geometric extensions, giving rise to several well-established frameworks. These include curvature-based modifications like $f(R)$ [6−8] and $f(G)$ gravity [9, 10], higher-order polynomial extensions such as cubic gravity [11−13], topological invariants in Lovelock gravity [14], and scalar-tensor couplings exemplified by Horndeski's theory [15] and its Galileon generalizations [16, 17]. This dichotomy underscores the central debate in modern cosmology: whether dark energy originates from new material constituents or from an extension of gravitational theory itself.

An alternative formulation of gravity, dynamically equivalent to GR at the level of field equations, employs torsion rather than curvature as the fundamental geometric descriptor. This approach, known as the Teleparallel Equivalent of General Relativity (TEGR) [18, 19], describes gravity through a torsion-based geometry where gravitational interactions result from the parallel transport of tetrad fields. In TEGR, the gravitational field is characterized by the torsion tensor with dynamics governed by the torsion scalar T, a quadratic contraction of the torsion tensor. A further equivalent representation of Einstein’s theory can be formulated in a flat, torsion-free geometry where gravitation is fully encoded in the non-metricity tensor $Q_{\alpha \mu \nu}=\nabla_\alpha g_{\mu\nu}$. This framework, known as symmetric teleparallel gravity, leads to $f(Q)$ gravity [20−22]; further details on this theory and its extensions can be found in [23−27]. Although TEGR reproduces the predictions of GR, it offers a distinct geometric interpretation [28−30], casting gravity as a manifestation of spacetime torsion rather than curvature. Nevertheless, like GR, TEGR alone does not resolve large-scale cosmological issues such as dark energy or inflation. To address these limitations, modified teleparallel theories, collectively referred to as $f(T)$ gravity, have been developed [31−33], generalizing the Lagrangian to arbitrary functions of T. As an extension of TEGR, $f(T)$ gravity opens new pathways for explaining cosmic acceleration and large-scale structure formation [34−37]. Still, whether $f(T)$ theories can outperform GR on both theoretical and observational grounds remains an open question, warranting further detailed investigation.

To construct a complete teleparallel analog of $f(R)$ gravity, the $f(T, B)$ extension plays an essential role [19, 38, 39]. In this framework, the torsion scalar T and the boundary term B respectively capture the second- and fourth-order derivative contributions present in $f(R)$ gravity. Notably, $f(T, B)$ naturally incorporates $f(R)$ gravity as the specific case $f(-T+B)$, while allowing for a wider range of gravitational Lagrangians.

The $f(T, B)$ framework has been widely explored across diverse phenomenological contexts. Gravitational wave studies in this theory indicate luminal propagation speeds and the presence of polarization modes beyond the standard transverse-traceless ones of GR [40, 41]. Solar-system tests further confirm the viability of many $f(T, B)$ models, showing agreement with high-precision astronomical measurements [42]. In cosmology, such models offer promising mechanisms to alleviate the Hubble tension [43−45]. Recent theoretical developments include: (i) a rigorous establishment of the correspondence between $f(T, B)$ and $f(R)$ gravity [19, 38], along with thermodynamic and cosmological reconstruction studies [55]; (ii) exact and perturbed black hole solutions [47], extending the theory's applicability to compact objects; and (iii) comprehensive analyses of background expansion and linear perturbation growth [48−52].

This study investigates the cosmological dynamics of $f(T,B)$ gravity within an isotropic, homogeneous FLRW universe using the dynamical system approach. The strong nonlinearity of the field equations renders exact solutions intractable and obstructs direct observational tests, owing to the complex coupling among terms. To address this, the dynamical system method [53−55], which has been used to investigate cosmological models in various modified gravity, including Hořava-Lifshitz gravity [56−59], $f(R)$ gravity [60, 61], $f(R, T)$ gravity [62, 63], Gauss-Bonnet gravity [64], Einstein cubic gravity [65], and other cosmological scenarios [66, 67]. This method reformulates the cosmological equations into an autonomous system, enabling the examination of critical points, phase-space trajectories, and stability, yielding global dynamical insight (including attractors and transient states) without relying on exact solutions. Such a global perspective is essential for comparing theories with observations, as it reveals the full dynamical landscape rather than individual solutions. Recently, using the dynamical system method, Kritpetch et al. [68] clarified the interaction mechanism between dark sector components in dark energy models incorporating both quintessence and phantom fields via a switching parameter. In parallel, Halder et al. [69] identified new stable accelerating scaling attractors within interacting phantom dark energy frameworks. These attractors offer a potential mechanism for alleviating the cosmic coincidence problem.

This study is structured as follows. Section 2 provides a brief review of teleparallel gravity and its extensions. In Section 3, we derive the cosmological dynamical system for the $f(T,B)$ model. Section 4 presents a dynamical analysis of the two considered models, and Section 5 examines their statefinder diagnostics. Finally, Section 6 summarizes the key findings and offers a concluding discussion.

This section provides a brief review of Teleparallel Gravity (TG) and its extension to $f(T, B)$ gravity. GR describes gravity using the Levi-Civita connection $\Gamma^a_{\mu\nu}$, which is characterized by non-zero curvature, zero torsion, and metric compatibility. In contrast, Teleparallel Gravity adopts the Weitzenböck connection $W_\mu{}^a{}_\nu$ – a curvature-free, metric-compatible connection that captures gravitational effects entirely through torsion [70]. This shift in geometric foundation has profound implications: while GR and its modifications employ the Riemann tensor to measure spacetime curvature, the Riemann tensor vanishes identically in TG due to the flatness of the Weitzenböck connection. As a result, TG necessitates a reconstruction of gravitational quantities from the ground up, offering novel theoretical possibilities while preserving dynamical equivalence with GR at the level of field equations.

In the TG framework, the fundamental dynamical variables are the tetrad fields (vierbeins) $e^a_\mu$, which form an orthonormal basis for the tangent space at each spacetime point $x^\mu$. The tetrads $e^a_\mu$ and their inverse fields $E^\mu_a$ satisfy the orthonormality conditions

where Latin indices $(m, n)$ refer to coordinates in the tangent space, and Greek indices $(\mu, \nu)$ denote spacetime coordinates. The metric tensor $g_{\mu\nu}$ is reconstructed from the tetrad fields via the relation

where $\eta_{ab}$ is the Minkowski metric on the tangent space. The Weitzenböck connection is defined as [19]

The torsion tensor $T^a_{\mu\nu}$ is given by the antisymmetric part of the Weitzenböck connection

Two key tensors in TG are the contorsion tensor $K_\mu{}^a{}_\nu$ and the superpotential $S_a{}^{\mu\nu}$. The contorsion tensor $K_\mu{}^\lambda{}_\nu$ is defined as

and plays a significant role in establishing the equivalence between TG and GR at the level of field equations. The superpotential $S_\sigma{}^{\mu\nu}$ is given by

The torsion scalar T is constructed through the complete contraction of the torsion tensor with its superpotential

This scalar serves as the Lagrangian density in the TEGR

where $\kappa^2=8\pi G$ and ${\cal{L}}_m$ is the matter Lagrangian. This quadratic combination encodes the teleparallel equivalent of the Ricci scalar R satisfying the identity

where $e=\mathrm{det}\left( e^a_\mu \right) \sqrt{-g}$, and the boundary term is defined as $B=\dfrac{2}{e}\partial_\mu (eT^\mu)=2\nabla_\mu T^\mu$.

A natural generalization of the TEGR action is obtained by promoting the torsion scalar T to an arbitrary function $f(T)$, leading to

Unlike the $f(R)$ framework, which leads to fourth-order field equations, $f(T)$ gravity retains second-order equations of motion. This distinction arises from the relaxed constraints of Lovelock's theorem in teleparallel geometry, where a torsion-based description allows for modifications that avoid ghost instabilities [29]. Nevertheless, the full structure of the theory involves two fundamental geometric scalars: the torsion scalar T and the boundary term B. Their relation to the Ricci scalar via the identity $R=-T+B$ motivates the generalization to $f(T, B)$ gravity, which not only encompasses $f(R)$ gravity as a special case but also provides a minimal extension incorporating both second- and fourth-order derivative terms.

In the present work, we consider the $f(T, B)$ action in the form [50]

Note that the TEGR action is recovered when $f(T,B)=0$. Varying this action with respect to the tetrad yields the field equations [19, 40]

where $f_T = \partial f /\partial T$, $f_B = \partial f /\partial B$, and $ \Theta^{\lambda}_\nu=e^a_\nu \Theta^\lambda_a $ denotes the standard energy-momentum tensor for matter. These field equations are derived under the assumption of a vanishing spin connection, which is a consistent choice in the context of a flat FLRW cosmology [19, 39, 40].

The choice of tetrad is

where $a(t)$ is the scale factor. This tetrad yields the flat FLRW metric

from which the torsion scalar T and the boundary term B are obtained as

where the overdot represents a derivative with respect to cosmic time t. The corresponding Ricci scalar R for this metric is thus recovered as $R=-T+B=6(\dot{H}+2H^2)$.

With the FLRW metric and the chosen tetrad, the field equations reduce to the modified Friedmann equations

where $\rho_m$ and $p_m$ respectively represent the energy density and pressure of the matter (baryons and dark matter) whose equation of state $\omega_m$ is defined as $p_m=\omega_m \rho_m$, while the geometric dark energy density $\rho_{\mathrm{gd}}$ and the corresponding pressure $p_{\mathrm{gd}}$ are given by

The conservation equation of the matter is

where we have assumed $p_m=0$, then the geometric dark energy also observes the conservation equation [39]

The equation of state for the geometric dark energy $\omega_\mathrm{gd}$ is defined as

while the total equation of state is given by

The deceleration parameter q is related to $\omega_\mathrm{tot}$ via

implying that the universe undergoes accelerated expansion when $q <0$, or equivalently, when $\omega_\mathrm{tot}<-\dfrac{1}{3}$.

In metric teleparallel gravity, with the affine connection $\Gamma^\alpha_{\mu \nu}$ vanishing in FLRW geometry, the boundary term satisfies $C = B$. Under this condition, the modified Friedmann equations of the present $f(T, B)$ model become equivalent to those of the $f(Q, C)$ theory [71, 72]. This equivalence indicates that the two models may share fundamental features in their cosmological dynamics. Moreover, the dynamical system analysis developed in the following sections for $f(T, B)$ gravity can be directly applied to the $f(Q, C)$ formulation when the affine connection is set to zero.

In this section, we perform a qualitative analysis of the cosmological dynamics in $f(T, B)$ gravity using the dynamical systems approach. By introducing suitable dimensionless variables, the modified field equations can be reformulated as an autonomous dynamical system. Such a system is characterized by two fundamental components: a state space comprising all possible configurations, and a set of differential equations governing the evolution of trajectories within this space.

The fixed points of the system, defined by the condition $\dot{X} = 0$ for a system of the form $\dot{X} = f(X)$ with $X = (x_1, x_2, \cdots, x_n)$, correspond to equilibrium solutions in the cosmological context. These critical points represent distinct cosmological epochs within the $f(T, B)$ framework. Their stability, determined via linear perturbation analysis, dictates the global evolutionary behavior of the universe: stable points act as cosmological attractors characterizing late-time asymptotic states, unstable points correspond to transient phases, and saddle points represent metastable regimes that temporarily influence the dynamics before the system evolves toward an attractor. Through this approach, key epochs in cosmic history – such as radiation domination, matter domination, and late-time acceleration – naturally arise as specific critical points in the phase space, revealing how different eras of universe evolution are embedded within the structure of $f(T, B)$ gravity.

To construct the dynamical system for the cosmological model, we introduce the following dimensionless variables:

Here, $\Omega_m$ quantifies the relative density of matter in the total effective cosmic fluid, while $\Omega_{\mathrm{gd}} = -x + 2u + \sigma - z$ represents the contribution from geometric dark energy within the $f(T,B)$ framework. The quantities x, u, σ, and z characterize different aspects of the geometric dark energy sector. Their relative dominance in the phase space can signal transitions between distinct dark energy-dominated regimes. At fixed points of the dynamical system, the values of these variables help identify the nature of the corresponding cosmological epoch, clarifying the physical role of each component in the evolution governed by the critical points.

Using the dimensionless variables defined in Eq. (21), the Friedmann equation (12) takes the form

while Eq. (13) becomes

The term $ \dot{f}_T/H $ can be expanded as

where $ f_{TT} = \partial^2 f / \partial T^2 $ and $ f_{TB} = \partial^2 f / (\partial T \partial B) $.

Furthermore, from the dynamical variables y and z, we derive the following relations:

in which $ f_{BT} = \partial^2 f / (\partial B \partial T) $, $ f_{BB} = \partial^2 f / \partial B^2 $, and it is assumed that $ f_{BB} \neq 0 $.

Finally, the dark energy equation of state $ \omega_{\mathrm{gd}} $, the total equation of state $ \omega_{\mathrm{tot}} $, and the deceleration parameter q are given by

To construct the cosmological dynamical system, we introduce the independent variable $ N \equiv \ln a $, commonly used in expanding cosmological scenarios ($ H > 0 $) but inapplicable in bouncing models where $ H = 0 $ at the bounce epoch [50]. The field equations of $ f(T, B) $ gravity can then be expressed as the following autonomous dynamical system:

where the prime symbol $ ' $ denotes differentiation with respect to $ N = \ln a $. Using the variables defined in Eq. (21) and the relations in Eqs. (22) and (23), the system can be reduced to:

Upon specifying the functional form of $ f(T, B) $, the dynamical system presented above becomes fully autonomous, in contrast to approaches that rely on the parameterization $ \lambda = \ddot{H}/H^3 $ as used in [49, 51, 52, 73]. In the subsequent sections, we focus on two specific $ f(T, B) $ models introduced in [39]: the power-law model

and the mixed power-law model

where $ c_1, c_2, c_3, \alpha, \beta $ are constant parameters of the cosmological model. Both of these representative and well-motivated prototype functional forms are chosen because they can naturally reduce to GR or to other established modified gravity theories, such as $ f(T) $ or $ f(R) $ gravity, within specific parameter limits, thereby ensuring theoretical consistency. They also serve complementary aims: the multiplicative power-law form is intended to explore novel dynamical effects that stem from a non-trivial coupling between the torsion scalar T and the boundary term B, an interaction inherently absent in additive or pure $ f(R) $ models. In contrast, the additive power-law form enables a clear separation and comparative assessment of the individual contributions of T and B to cosmic evolution. Furthermore, the power-law ansatz yields homogeneous terms in the resulting Friedmann equations, which significantly facilitates the search for exact scaling solutions and the construction of a closed autonomous dynamical system. It should also be noted that other forms of $ f(T,B) $ may be considered in future studies, such as $ f(T,B)=A_0 + A_1 T + A_2 T^2 + A_3 B + A_4 T B $ [40, 49], $ f(T,B) = \xi T+ \alpha B\ln B $ [52], as well as $ f(T,B) = B g(T) $ and $ f(T,B) = T g(B) $ [50]. By systematically analyzing the two foundational ansätze selected here, we aim to map the key dynamical features of $ f(T,B) $ cosmology and establish a benchmark for future studies involving more complex functional dependencies.

We first consider the power law model $ f(T, B)= c_1 T^\alpha B^\beta $. For Eq. (26) to be well-defined, the condition $ f_{BB} \neq 0 $ must be satisfied, which requires $ \beta \neq 0, 1 $ and $ c_1 \neq 0 $. In this case, the dynamical variables u and v can be written as

Substituting these into Eq. (26) yields

The cosmological dynamical system for the power law model then reduces to the autonomous form

The autonomous system (32) admits a unique fixed point, denoted as $ P^1_{\mathrm{gd}} $. Its coordinates and the corresponding cosmological parameters are summarized in Table 1. The eigenvalues of the Jacobian matrix evaluated at this point, denoted as $ \{e_{a1}, e_{a2}, e_{a3}\} $, are given by

where

Table 2 summarizes the existence and linear stability conditions for $ P^1_{\mathrm{gd}} $, along with its acceleration behavior, where the symbols $ \wedge $ and $ \vee $ denote logical “and” and “or”, respectively. For the specific parameter values $ \alpha = 3 $ and $ \beta = -4 $, the fixed point $ P^1_{\mathrm{gd}} $ is a stable node. The phase space stream plot of the model for this parameter set is shown in Figure 1, while the evolution of the corresponding cosmological parameters is displayed in Figure 2.

Figure 1.  (color online) Phase space flow of the model $ f(T, B)=c_1 T^\alpha B^\beta $ for $ (\alpha_1, \beta_1)=(3,-4) $.

Figure 2.  (color online) Evolution of cosmological parameters for $ (\alpha_1, \beta_1)=(3,-4) $ and initial conditions $ (\Omega_{m_0}, x_0, y_0)= $$ (0.3, 0.3, 2.5) $.

The fixed point $ P^1_{\mathrm{gd}} $ corresponds to a cosmological epoch dominated by geometric dark energy, characterized by a de Sitter expansion $ a \sim e^{H_0 t} $ and a total equation of state $ \omega_{\mathrm{tot}} = -1 $. At this point, the geometric dark energy components x, u, and σ collectively sustain the accelerated expansion. Owing to its stable nodal behavior for suitable parameter choices $ (\alpha, \beta) $, this fixed point provides a viable mechanism for explaining the late-time cosmic acceleration within the power-law $ f(T, B) $ framework.

We now consider the mixed power law model $ f(T, B) = c_2 T^\alpha + c_3 B^\beta $. Here, the condition $ f_{BB} \neq 0 $ is satisfied provided that $ \beta \neq 0, 1 $ and $ c_3 \neq 0 $. In this model, the variable x and the term $ \dot{B}/(6H^3) $ can be written as

The corresponding autonomous dynamical system takes the form

This system possesses two fixed points, denoted as $ P^i_{\mathrm{gd}} $ with coordinates $ (\Omega^i_m, u^i, y^i, v^i) $ for $ i = 2, 3 $.

The fixed point $ P^2_{\mathrm{gd}} $ is characterized by the coordinates

where $ v_* \in \mathbb{R} $ and $ v_* \neq 0 $. More precisely, $ P^2_{\mathrm{gd}} $ corresponds to a line of equilibrium points. The eigenvalues of the linearized system at this point are given by

with

At this fixed point, the cosmological parameters satisfy $ \Omega_m = 0 $ and $ \Omega_{\mathrm{gd}} = 1 $, indicating a universe dominated by geometric dark energy. As a de Sitter point with a scale factor evolving as $ a \sim e^{H_0 t} $, $ P^2_{\mathrm{gd}} $ provides a potential explanation for the current cosmic acceleration within the model, provided it acts as a stable attractor.

The fixed point $ P^3_{\mathrm{gd}} $ is located at

This point belongs to the equilibrium line $ P^2_{\mathrm{gd}} $, corresponding to the specific case where $ v_* = \dfrac{\beta}{3(\beta - 1)} $. The eigenvalues of the linearization at $ P^3_{\mathrm{gd}} $ are

Like $ P^2_{\mathrm{gd}} $, this point also corresponds to a phase of geometric dark energy dominance and exponential expansion of the universe. If linearly stable, it could provide a mechanism for late-time cosmic acceleration.

The coordinates and cosmological parameters of both fixed points are summarized in Table 3, while their existence, stability, and acceleration properties are listed in Table 4, where $ \tilde{v} = \beta(\alpha - 1)/(3(\beta - 1)(\alpha - \beta)) $. In the mixed power-law model, both obtained fixed points exhibit one zero eigenvalue and are therefore non-hyperbolic. The central manifold theorem was attempted to assess their stability; however, after decomposing the system into linear and nonlinear components, it was found that the nonlinear terms do not vanish in the vicinity of the equilibrium, thus precluding definitive stability conclusions through this method. Consequently, we specify the condition that the real parts of the eigenvalues $ e_{j3} $, $ e_{j4} $ (with $ j = b, c $) must be negative in Table 4, and present the stability behavior of both points via phase portrait analysis in Figure 3.

Figure 3.  (color online) Phase space flow of the mixed power-law model under different parameters: (a), (c): $ (\alpha_2, \beta_2, v_2)= $$ (5,-1000, 0.00125) $; (b): $ (\alpha_3, \beta_3, v_3)=(5, 2/3, 1) $; (d): $ (\alpha_3, \beta_3)=(5, 2/3) $.

We perform numerical analysis for two representative parameter pairs: $ (\alpha_2, \beta_2) = (5, -1000) $ and $ (\alpha_3, \beta_3) = (5, 2/3) $, as illustrated in Figures 3 and 4. In the y-u plane, $ P^2_{\mathrm{gd}} $ is stable under both parameter choices and exhibits spiral dynamics for $ (5, -1000) $. The cosmological parameters $ \Omega_m $, $ \Omega_{\mathrm{gd}} $, $ \omega_{\mathrm{tot}} $, and $ \omega_{\mathrm{gd}} $ all asymptotically approach the same final state $ (0, 1, -1, -1) $. However, for the pair $ (5, 2/3) $, the evolution of $ \omega_{\mathrm{tot}} $ and $ \omega_{\mathrm{gd}} $ shows distinct damped oscillations before settling to the de Sitter attractor.

Figure 4.  (color online) Evolution of cosmological parameters for the mixed power-law model from initial conditions $ (\Omega_{m_0}, u_0, y_0, v_0)=(0.3, 1, 2.5, 1) $. (a): $ (\alpha_2, \beta_2)=(5,-1000) $; (b): $ (\alpha_3, \beta_3)=(5,2/3) $.

This section applies the statefinder diagnostic to the two cosmological models introduced above. The statefinder parameters r and s, first proposed in [74, 75], provide a useful tool for distinguishing between different dark energy scenarios. These dimensionless quantities are defined in terms of the scale factor and its higher-order derivatives as [74]:

The statefinder diagnostic offers a practical framework for classifying alternative dark energy models. Even when two models predict similar expansion histories, their evolutionary trajectories in the r-s plane can reveal clear differences. Characteristic values of $ \{r, s\} $ for several standard dark energy models include [74, 75]:

$ \bullet $ $ \{r = 1, s = 1\} $: Standard cold dark matter (SCDM) model

$ \bullet $ $ \{r = 1, s = 0\} $: ΛCDM model

$ \bullet $ $ \{r > 1, s < 0\} $: Chaplygin gas model

$ \bullet $ $ \{r < 1, s > 0\} $: Quintessence model

The application of statefinder diagnostics within a dynamical systems framework proves particularly effective in distinguishing between different acceleration regimes, as demonstrated in recent studies of dark energy models [76, 77]. In present work, for the power-law model $ f(T, B) = c_1 T^\alpha B^\beta $, the statefinder parameters take the form

whereas for the mixed power-law model $ f(T, B) = c_2 T^\alpha + c_3 B^\beta $, they are given by

In both cases, the deceleration parameter is expressed as $ q = 2 - y $.

The evolutionary trajectories of the statefinder parameters r and s are illustrated in the r-q and r-s planes: Figure 5 corresponds to the power-law model $ f(T, B) = c_1 T^\alpha B^\beta $, and Figure 6 to the mixed power-law model $ f(T, B) = c_2 T^\alpha + c_3 B^\beta $.

Figure 5.  (color online) Evolution of statefinder parameters in the r-q and r-s planes for the model $ f(T, B) = c_1 T^\alpha B^\beta $, with initial values $ (\Omega_{m_0}, x_0, y_0) = (0.3, 0.3, 2.5) $ and parameter choice $ (\alpha_1, \beta_1) = (3, -4) $.

Figure 6.  (color online) Evolution of statefinder parameters in the r-q and r-s planes for the model $ f(T, B) = c_2 T^\alpha + c_3 B^\beta $, with initial values $ (\Omega_{m_0}, u_0, y_0, v_0) = (0.3, 1, 2.5, 1) $. Panels (a) and (b): $ (\alpha_2, \beta_2) = (5, -1000) $; panels (c) and (d): $ (\alpha_3, \beta_3) = (5, 2/3) $.

As shown in Figure 5, the power-law model exhibits dynamical behavior analogous to that of a Chaplygin gas model before asymptotically approaching a ΛCDM-like regime. In Figure 6, although the mixed power-law model also converges to a ΛCDM-like state in the late-time limit, the evolutionary paths in the r-q and r-s planes differ significantly between the two parameter sets. For $ (\alpha, \beta) = (5, -1000) $, the trajectory remains within the region characteristic of the Chaplygin gas model before reaching the ΛCDM point $ \{r = 1, s = 0\} $. In contrast, for $ (\alpha, \beta) = (5, 2/3) $, the trajectory crosses the ΛCDM point, showing damped transitions between the Chaplygin gas and quintessence regimes. Additionally, the trajectory for $ (5, 2/3) $ displays a spiral structure in the r-q plane.

In particular, the distinct evolutionary trajectories observed in the statefinder diagnostic planes highlight the sensitivity and discriminatory power of the $ \left\{ r, q \right\} $ and $ \left\{ r, s\right\} $ diagnostics in distinguishing between multiplicative power-law and additive mixed power-law forms of $ f(T,B) $ models. The notable differences in both the paths and transitional behaviors between these two functional forms can be attributed to their distinct gravitational dynamics. The multiplicative form $ f(T,B) = c_1 T^\alpha B^\beta $ introduces a strong and inseparable coupling between torsion and the boundary term, which constrains the evolutionary path to remain within Chaplygin gas-like regimes until late times. Conversely, the additive form $ f(T,B) = c_2 T^\alpha + c_3 B^\beta $ decouples the contributions of T and B, permitting more varied interactions. This decoupling enables transitions between Chaplygin gas and quintessence behaviors, and can even lead to oscillatory or spiral approaches to the ΛCDM attractor, as evident in the r-q plane for the parameter set $ \left( \alpha, \beta \right) $ = $ \left( 5, 2/3 \right) $. Therefore, the statefinder analysis not only effectively distinguishes between the two $ f(T,B) $ ansätzes but also reveals how the structural choice (namely, whether to couple or to separate T and B) fundamentally shapes the dynamical character and potential transient phases of cosmic acceleration.

This study systematically analyzes the cosmological dynamics of two well-motivated models within the framework of $ f(T,B) $ modified gravity. Focusing on the multiplicative power-law form $ f(T, B) = c_1 T^\alpha B^\beta $ and the additive mixed power-law form $ f(T, B) = c_2 T^\alpha + c_3 B^\beta $, we examine how the coupling and decoupling of the torsion scalar T and the boundary term B shape the late-time evolution of a flat FLRW universe.

By constructing autonomous systems for both models, we identify stable de Sitter-type fixed points that act as late-time attractors, providing a purely geometric explanation for cosmic acceleration. The multiplicative power-law model shows a smooth convergence toward a ΛCDM-like state via an intermediate Chaplygin gas regime. In contrast, the additive mixed power-law model displays richer dynamical behavior, including damped oscillations and spiral trajectories in the statefinder planes, as illustrated for parameters such as $ (\alpha, \beta) = (5, 2/3) $. Moreover, statefinder diagnostics, specifically the r-s and r-q planes, effectively distinguish each model from the other and from the standard ΛCDM scenario, highlighting observationally testable features.

Methodologically, our analysis differs from several earlier dynamical studies in modified gravity, such as those in [49, 51, 52, 73], which commonly adopt the parameterization $ \lambda = \ddot{H}/H^3 $. Instead, we introduce the auxiliary variables y and v, which ensure the autonomy of the dynamical system without requiring additional phenomenological assumptions.

In comparison with $ f(T) $ and $ f(R) $ cosmologies, the present work highlights how the structural choice between multiplicative and additive coupling of T and B qualitatively influences the dynamical landscape. Unlike $ f(T) $ or $ f(R) $ models, the additive mixed power-law form of $ f(T,B) $ permits richer transitional behaviors that are less common in simpler frameworks. These findings indicate that $ f(T,B) $ gravity, especially in its additive form, provides a more flexible phenomenological framework for describing dynamical dark energy while naturally accounting for late-time cosmic acceleration.

Future research could naturally build on this foundation in several ways. Extending the dynamical analysis to explicitly include matter-dominated phases would provide a more complete description of cosmic evolution. The models should also be tested against a broader set of observational data, such as cosmic chronometers, baryon acoustic oscillations, and the growth of large-scale structure, to better constrain the parameters α, β, and $ c_i $. Furthermore, exploring more general functional forms of $ f(T,B) $ (for example, logarithmic, exponential, or piecewise-defined combinations) could help assess the robustness of the dynamical features identified in this work. Another fruitful avenue would be to examine the implications of such models for early-universe cosmology, including scenarios related to inflation and singularity avoidance.

The authors wish to express their sincere gratitude to the anonymous reviewers for their invaluable comments and constructive suggestions, which have significantly contributed to enhancing the quality and clarity of this work.