Qualitative behaviour of higher-curvature gravity with boundary terms i.e the f(Q) gravity models by dynamical system analysis

The higher-curvature gravity with boundary terms i.e the $f(Q)$ theories, grounded on non-metricity as a fundamental geometric quantity, exhibit remarkable efficacy in portraying late-time universe phenomena. The aim is to delineate constraints on two prevalent models within this framework, namely the Log-square-root model and the Hyperbolic tangent-power model, by employing the framework of Big Bang Nucleosynthesis (BBN). The approach involves elucidating deviations induced by higher-curvature gravity with boundary terms in the freeze-out temperature ($T_{f}$) concerning its departure from the standard $\Lambda$CDM evolution. Subsequently, constraints on pertinent model parameters are established by imposing limitations on $\vert \frac{\delta T_{f}}{T_{f}}\vert$ derived from observational bounds. This investigation employs dynamical system analysis, scrutinizing both background and perturbed equations. The study systematically explores the phase space of the models, identifying equilibrium points, evaluating their stability, and comprehending the system's trajectory around each critical point. The principal findings of this analysis reveal the presence of a matter-dominated saddle point characterized by the appropriate matter perturbation growth rate. Subsequently, this phase transitions into a stable phase of a dark-energy-dominated, accelerating universe, marked by consistent matter perturbations. Overall, the study substantiates observational confrontations, affirming the potential of higher-curvature gravity with boundary terms as a promising alternative to the $\Lambda$CDM concordance model. The methodological approach underscores the significance of dynamical systems as an independent means to validate and comprehend the cosmological implications of these theories.


I. INTRODUCTION
The quest to comprehend the prevailing cosmic acceleration stands as a paramount and intricate challenge in contemporary theoretical physics, especially within the realm of cosmology.The current phase of cosmic expansion, marked by an unprecedented acceleration, presents a puzzle demanding explanation.Numerous observational avenues, starting with type Ia Supernova observations [1][2][3], have provided initial indications of this phenomenon.Subsequently, additional support has been garnered from a convergence of evidence derived from the Cosmic Microwave Background (CMB) radiation, Large Scale Structures, data from the Wilkinson Microwave Anisotropy Probe (WMAP) [4,5], and insights obtained from Baryonic Acoustic Oscillations (BAO) [6,7].Collectively, these empirical observations serve as pillars reinforcing the understanding that our universe is undergoing an accelerated expansion, a paradigm shift in our comprehension of cosmic dynamics.
There are generally two kinds to the explanations that have been put forward for the phenomenon that has been observed.One is to assume is the existence of an exotic energy with negative pressure known as dark energy.There are numerous candidates for dark energy at the moment, including the cosmological constant, quintessence, phantom, quintom, and so on.The cosmological constant model (ΛCDM), with the equation of state w Λ = −1, is the most competitive cosmic dark energy model [8][9][10][11].However, it has two significant theoretical issues, namely the coincidence problem and the cosmological constant problem [12,13].As a result, a few scalar field models are proposed, including quintessence [14] and phantom [15][16][17].It has been shown that the equation of state for models with a single scalar field cannot pass the phantom divide line (w Λ = −1).Recent observations have hinted at the plausibility of traversing the phantom divide line within cosmological models [18][19][20].This intriguing possibility has spurred the establishment of various theoretical frameworks aimed at realizing this phenomenon.Models featuring a blend of phantom and quintessence [21][22][23][24] components have emerged, seeking to elucidate the crossing of this divide.Additionally, scalar field models incorporating scalar-dependent couplings in conjunction with the kinetic term have been proposed.Furthermore, fluid models, characterized by their unique formulations, have been introduced to explore and potentially explain the crossing of this crucial cosmic threshold [25][26][27].Modifying Einstein's general relativity theory is another technique for explaining the current, accelerating cosmic expansion.These modifications often manifest through adjustments to the Einstein-Hilbert action, primarily employing altered Lagrangian densities.Beyond non-Lagrangian models such as Modified Newtonian Dynamics (MOND), a plethora of modified gravity theories currently exist.Notably, f (R), f (R, T ), f (T ), and higher-curvature gravity with boundary terms i.e the f (Q) gravity theories stand out as the most cosmically viable among these frameworks.In f (R) theory, the conventional Einstein-Hilbert action's Ricci scalar R is replaced by an arbitrary function f (R).Remarkably, the resultant effective equation of state showcases the capacity to traverse the phantom divide, transitioning from a phantom phase to a non-phantom one [28][29][30][31].Whereas, f (R, T ) theory extends the gravitational action by incorporating supplementary terms involving both the Ricci scalar (R) and the energy-momentum tensor (T ), the latter delineating the distribution of matter and energy across spacetime [32][33][34].
Alternative approaches to defining gravitational interactions encompass torsion and non-metricity, leading to the emergence of GR analogues, namely the f (T ) and f (Q) gravities in modified forms [35][36][37].These variations pivot on unconventional metric-affine connections, such as the Weitzenböck and metric-incompatible connections, divergent from the established Levi-Cevita connection within general relativity.Utilizing these connections facilitates the derivation of both f (T ) and f (Q) gravities.In teleparallelism, the Weitzenböck connection correlates with torsion possessing zero curvature, while in general relativity, the Levi-Civita connection is associated with curvature devoid of torsion [38,39].Notably, f (T ) gravity aligns as the teleparallel equivalent of the General Relativity (TEGR) adaptation, known for its comprehensibility and ease of conceptualization.Recent scrutiny has focused on the Symmetric Teleparallel Equivalent of General Relativity (STEGR) as a novel gravitational theory.STEGR delineates gravitational interactions through the construct of non-metricity (Q), characterized by both zero torsion and curvature.Noteworthy modifications to STEGR involve higher-curvature gravity with boundary terms, representing a refined iteration known as modified Symmetric Teleparallel Gravity or f (Q) gravity, serving as the focal point of exploration in this paper.
The complexities inherent in gravitational theories, stemming from non-linear elements within field equations, pose a substantial challenge in obtaining solutions-whether analytical or numerical-hindering direct comparisons with observational data.To address this, employing "Dynamical System Analysis" emerges as a valuable technique for elucidating, solving, or comprehensively examining the general behavior of these equations [40][41][42][43][44].The essence of dynamical system analysis lies in seeking numerical solutions that encapsulate the qualitative behavior of a given physical system [45][46][47].This method primarily involves identifying critical points within a set of first-order ordinary differential equations.Determining the system's behavior around these critical points involves linearizing the system post-calculation of the Jacobian matrix at each critical point and its eigenvalues, thus deriving stability conditions.This pursuit centers on investigating the stability traits surrounding a specific critical point [48][49][50].Linear stability theory, while effective for hyperbolic points (where none of the eigenvalues of the Jacobian matrix have a zero real part), encounters limitations with non-hyperbolic points.In such cases, alternative methods like Lyapunov Functions and Center Manifold Theory become necessary to explore and analyze stability properties.The present study employs "Center Manifold Theory" to investigate and comprehend stability properties surrounding non-hyperbolic critical points within the gravitational field equations, thus contributing to a deeper understanding of their dynamic behaviors.
This investigation focuses on exploring the cosmic dynamics of f (Q) gravity models by employing dynamical system analysis as a mathematical tool at background and perturbation levels, this study aims to contribute insights that can potentially support observational analyses.The structure of this work unfolds as follows: Section II offers a formal delineation of Symmetric Teleparallel gravity, elucidating the field equations of f (Q).This section provides the foundational background necessary for deriving both background and perturbed cosmological equations.Section III delves into an examination of Center Manifold Theory, elucidating its relevance within this context.Subsequently, Section IV introduces the formalism of Big Bang Nucleosynthesis (BBN) and its observational implications, essential for constraining various classes of f (Q) models.Building upon these constraints, Section V delves into a discussion regarding the application of BBN constraints specifically within the framework of f (Q) gravity.Section VI then presents an analytical exploration of the phase space for the Log-square-root model and the Hyperbolic tangent-power model.This involves identifying equilibrium points, assessing their stability, and comprehending the system's behavior within the phase space.Finally, Section VII encapsulates an extensive summary of the paper.It incorporates a comprehensive discussion concerning the results obtained, their interpretations, and outlines prospective avenues for future exploration within this domain.

II. SYMMETRIC TELEPARALLEL GRAVITY
General Relativity is founded upon Lorentzian Geometry, which is established through the selection of a link that adheres to both symmetry and metric compatibility criteria.Within this framework, the Levi-Civita connection serves as the chosen link, giving rise to non-zero curvature exclusively owing to its inherent properties, while maintaining zero values for torsion and nonmetricity [51][52][53][54][55].However, when employing geometrodynamics as the foundational mathematical theory for gravity, the utilization of novel types of connections becomes feasible.Notably, the most comprehensive connection in this context is termed metric-affine and is expressed by the formula where the Christoffel symbols of the Levi-Civita connec-tion { β µν } is defined by µν and K β µν is the deformation and the contortion defined by with the torsion tensor T β µν defined as the anti-symmetric part of the affine connection, In conclusion, the three geometric deformations { β µν }, L β µν , and K β µν provide a sort of "trinity of gravity" that includes all elements of the general connection Γ β µν [56].This effectively shows that torsion, non-metricity, or curvature can be used to express the geometry of a gravitational theory.Symmetric Teleparallel Gravity connections are those that accept only non-metricity when both curvature and torsion are zero, whereas Teleparallel connections only accept zero curvature.The non-metricity tensor is given by [51,52] introduced the concept of symmetric teleparallel gravity, also identified as f (Q) gravity.This modified gravitational framework is characterized by an action defined as follows: wherein g represents the determinant of the metric tensor g µν , while L m denotes the Lagrangian density associated with matter.The function f (Q) in this formulation stands as an arbitrary expression of Q, the non-metricity scalar that governs gravitational interactions within this context.The non-metricity scalar Q is written as [51][52][53] where the quantities Q α ≡ Q α µ µ and Qα ≡ Q µα µ represent two independent traces derived by contracting the non-metricity tensor, expressed as . This tensorial manipulation allows for the definition of the non-metricity scalar, given by Q = −Q αµν Q.The corresponding field equation is given by with the conjugate to f (Q) is defined as T µν is the matter energy-momentum tensor, whose form is and f Q = ∂f ∂Q .To find the field equations, we set the action in (3) as constant with respect to variations over the metric tensor At the background level, the analysis assumes a spatially flat Friedmann-Lemaitre-Robertson-Walker (FLRW) spacetime characterized by homogeneity and isotropy.The metric describing this spacetime takes the form: where N (t) represents the Lapse function, and t denotes cosmic time.Within Q theories, a residual time reparameterization invariance is retained, permitting the selection of a specific form for the Lapse function.In the context of the non-metricity scalar, specifically within the FLRW metric, the expression Q = 6H 2 N 2 emerges.Here, H = ȧ a signifies the Hubble function, indicative of the universe's expansion rate at time t, where '.' denotes differentiation with respect to t.The function a(t) represents the scale factor, quantifying the universe's size at time t, while x, y, z denote Cartesian coordinates.Through symmetry considerations, setting N (t) = 1 yields the expression Q = 6H 2 within the FLRW metric, aligning with the framework's inherent characteristics and simplifying the analysis within this spacetime context.Applying the FLRW metric, the corresponding field equations are [51,53,[57][58][59][60] , The quantities f Q = df dQ and f QQ = d 2 f dQ 2 represent derivatives of the function f with respect to the non-metricity scalar Q.Additionally, ρ m , ρ r , p m , and p r denote the energy densities and pressures associated with perfect fluids characterizing matter and radiation.These quantities adhere to the energy conservation equations, expressed as follows in the absence of interaction: ρ m + 3H(ρ m + p m ) = 0 (11) ρr + 3H(ρ r + p r ) = 0 (12) Here, ρ m and ρr denote the time derivatives of matter and radiation energy densities, respectively, while H represents the Hubble parameter.Furthermore, these perfect fluids adhere to a linear equation of state parameter linking their pressures to their energy densities, given by p = ρw, where w belongs to the interval [−1, 1].This equation of state parameter w characterizes the relationship between pressure and energy density for the matter and radiation components.We can re-write the equations ( 9)-( 10) as where ρ DE and p DE are the energy density and pressure of effective dark energy.From equation ( 9)- (10) and equation ( 13)-( 14), the energy density and pressure of effective dark energy sector can be defined as For a universe experiencing acceleration, the condition is necessary [51,61].Consequently, introducing energy density parameters for different sectors becomes advantageous: where the index "i" represents matter, radiation, and dark energy.To facilitate the analysis, assuming the decomposition f (Q) = Q + F (Q) and setting 8πG = 1 for simplicity, the corresponding field equations ( 9)-( 10) can be expressed as: Here, F Q = dF dQ and F QQ = d 2 F dQ 2 .These equations describe the evolution of the Hubble parameter H and the effective energy density ρ, encompassing contributions from the non-metricity scalar Q and the function F (Q), providing a comprehensive understanding of the universe's dynamics within this theoretical framework.From equation ( 18), we have Hence the Friedmann's equation( 18) can be simply written as S where and Introducing the effective total energy density ρ ef f and total energy pressure p ef f , as [55] the corresponding total equation of state w ef f can be written as In the investigation of linear perturbations, our focus lies on the matter density contrast δ = δρ ρ , wherein δ ρ denotes the perturbation in the matter energy density.The governing equation for the evolution of the matter over density in the quasi-static regime, as presented in [59], is expressed as: where the denominator on the right-hand side elucidates the emergence of an effective Newtonian constant.Notably, within the context of minute scales significantly confined within the cosmic horizon, the terms associated with temporal derivatives in the perturbation equations are neglected, thereby reducing the formulation to predominantly encompass spatial derivative terms [55,62].

III. CENTER MANIFOLD THEORY
The stability analysis of orbits within the center manifold relies on three fundamental theorems.The initial theorem establishes the existence of the center manifold, providing a foundational basis for further investigation.Subsequently, the second theorem is dedicated to addressing stability considerations within this manifold.Finally, the third theorem delineates the process of locally constructing the real center manifold and underscores its adequacy in exploring stability properties [54,57].Consider the dynamical system which can be rewritten in the form where (µ, ν) ∈ R α × R β with α and β, the dimension of E α and E β respectively , and the functions ϕ and ψ satisfies where ∇ represents the gradient operator, E α and E β denote the center and stable subspaces correspondingly.These subspaces are characterized by the eigenvectors stemming from the Jacobian matrix, where the eigenvectors associated with eigenvalues having zero real parts pertain to E α , while those with negative real parts align with E β .Furthermore, within the framework of system (25), matrices A and B exhibit square forms with dimensions α × α and β × β, respectively.Their eigenvalues correspond to zero real parts for matrix A and negative real parts for matrix B. A geometrical space is a centre manifold for (25) if it can be locally represented as for δ sufficiently small and h(µ) is (sufficiently regular) function on R β

Definition 1 (Stable Fixed Point)
A fixed point ξ 0 in the context of the dynamical system ξ ′ = F (ξ) is deemed stable if, given any arbitrary small positive value ε, it is possible to determine a corresponding positive value δ such that for any solution η(t) of the system ξ ′ = F (ξ) satisfying the condition η(t 0 ) − ξ 0 < δ, the solution η(t) exists and remains defined for all t ≥ t 0 , ensuring that η(t) − ξ 0 < ε holds for all t ≥ t 0 .
In simpler terms, a fixed point ξ 0 is considered stable when all solutions ξ(t) starting in close proximity to ξ 0 remain in the vicinity of this point throughout their evolution.

Definition 2 (Asymptotically Stable Fixed Point)
A fixed point ξ 0 within the dynamical system ξ ′ = F (ξ) attains the status of asymptotic stability when, for any given positive value ε, there exists a corresponding positive value δ such that for any solution η(t) of the system ξ ′ = F (ξ) meeting the condition η(t 0 ) − ξ 0 < δ, the limit as time tends towards infinity (lim t→∞ η(t)) equals ξ 0 .
In simpler terms, a fixed point ξ 0 is considered asymptotically stable when it possesses stability characteristics and all solutions, initiated from nearby initial conditions, ultimately converge towards this fixed point as time progresses.
In physical terms, when a system exhibits stability and demonstrates a gradual reduction of perturbations from its equilibrium state over time, it is termed "asymptotically stable."This characterization signifies that subsequent to a disturbance, the system not only reverts to its equilibrium but also undergoes a convergent process towards it.As time tends toward infinity, these deviations diminish, elucidating the system's propensity to approach and ultimately achieve its equilibrium state.

Theorem 1 (Existence)
There exists a centre manifold for equation (25).The dynamics of the system (25) restricted to the centre manifold is given by for υ ∈ R α is sufficiently small.

Theorem 2 (Stability)
Assuming the existence of a stable zero solution (whether asymptotically stable or unstable) for equation (25), it follows that the zero solution of equation ( 26) is also stable (either asymptotically stable or unstable).Additionally, if the solution (µ(t), ν(t)) satisfies equation (25) with initial conditions (µ(0), ν(0)) that are suitably small, there exists a solution υ(t) of equation ( 26) such that and as t → ∞, where δ > 0 represents a constant.

Theorem 3 (Approximation)
Let Ψ : R α → R β be a mapping with Ψ(0) = ∇Ψ(0) = 0 such that N (Ψ(µ)) = O(|µ| q ) as µ → 0 for some q > 0. Then The Big Bang Nucleosynthesis (BBN) constraints formalism is discussed in this section.The radiation period sees the realisation of BBN [63][64][65][66].In the context of Big Bang Nucleosynthesis (BBN) and standard cosmology within the framework of the Standard Model of particle physics and general relativity, an approximation is made to the first Friedmann equation as detailed in reference [68]: During the radiation-dominated epoch, crucial in BBN, the consideration of the energy density of relativistic particles becomes imperative.This density, denoted by ρ r , is expressed as: Here, T signifies the temperature, while the parameter g * , approximately valued at 10, represents the effective number of degrees of freedom accounting for relativistic particles [66,67].Hence from equation ( 27) and ( 28), we obtain where is the Planck mass.We can extract the expression between temperature and time, since the radiation conservation equation finally results in a scale factor evolution of the form a ∼ t 1 2 , namely In the context of Big Bang Nucleosynthesis (BBN), the calculation of neutron abundance involves considering the rate of conversion between protons and neutrons, as referenced in [65,66].This conversion rate λ pn (T ) is defined as a combination of various processes: The total rate, denoted as λ tot (T ), is the summation of λ np (T ) and λ pn (T ), where λ np (T ) represents the inverse of λ pn (T ).Upon straightforward calculations, an expression for λ tot (T ) is obtained, as referenced in [67,[69][70][71][72]: Here, M = m n − m p = 1.29 × 10 −3 GeV signifies the mass difference between the proton and neutron, while The assumptions underlying this calculation entail equality in temperatures among different particles (neutrinos, electrons, and photons) and the application of Boltzmann distribution, assuming sufficiently low temperatures to use this distribution instead of the Fermi-Dirac distribution.The comprehensive framework of BBN details can be found in the Appendix of [67].
The determination of the freeze-out temperature (T f ) involves a comparison between the expansion rate of the universe ( 1H ) and the total interaction rate λ tot (T ).A significant distinction arises based on the relative magnitudes: if 1  H is notably smaller than λ tot (T ), it signifies that all processes occur within a state of thermal equilibrium, as cited in [63,64].Conversely, when 1 H significantly exceeds λ tot (T ), particles lack adequate time intervals for interactions to occur.The freezeout temperature T f delineates the temperature at which particles undergo a disconnection phase, and it corresponds to the moment when H = λ tot (T ) ≈ qT 5 , with q = 4A4!≃ 9.8 × 10 −10 GeV −4 , as referenced in [67,[69][70][71][72]. Combining equations ( 29) and ( 31), the expression for T f is derived: Here, M P denotes the Planck mass, and g * represents the effective number of degrees of freedom.This determination holds significance in understanding the temperature regime where particle interactions cease, leading to the transition from equilibrium to non-equilibrium conditions in the early universe.
In the epoch of Big Bang Nucleosynthesis (BBN), Modified gravity introduces an effective dark energy density, denoted as ρ DE , which coexists alongside other components, particularly radiation density (ρ r ).Considering ρ DE to be smaller than ρ r , it is deemed a first-order deviation.Consequently, the Hubble function is modified as follows: where H GR represents the Hubble parameter in standard cosmology, and the deviation from this standard scenario is encapsulated by: This deviation, δH, induces a consequential shift in the freeze-out temperature, denoted as δT f .This alteration signifies the impact of modified gravity, particularly the introduced effective dark energy, on the critical temperature where particle interactions cease during BBN.However, the relationship H GR = λ tot ≈ qT 5 f leads to an expression for the deviation δH as . Considering ρ DE to be significantly smaller than ρ r during this period, the equation simplifies to: This theoretically derived ratio δT f T f needs comparison with an observational constraint: This constraint is established through observational assessments of the baryon mass fraction converted into 4 H e , as documented in [73][74][75][76][77][78][79].In forthcoming subsections, this formalism, specifically expressed by equation (33), is applied to impose limitations on ρ DE and consequently on specific models of modified gravity.This process aids in understanding the bounds of the effective dark energy and its implications for various scenarios within modified gravity frameworks.

V. BBN CONSTRAINTS ON THE HIGHER-CURVATURE GRAVITY WITH BOUNDARY TERMS
In this section, the models, i.e the Log-square- ) will be tested against the constraint (34).This will be realized by calculating (33) for each of the models.

Log-square-root model
In accordance with the methodology outlined in [55,80], we introduce a logarithmic-square-root model represented by f (Q): Q ) Here, n and λ > 0 are the parameters characterizing the model.Correspondingly, the associated effective energy density (designated as Equation ( 15)) is expressed as: The first Friedmann equation, evaluated at the present time, leads to a relationship determining the parameter n: Here, Ω F 0 = 1 − Ω m0 − Ω r0 , denoting the present value of a quantity, with subscripts indicating the present time.Substituting equations ( 35) and ( 36) into (33) yields: Upon incorporating equation ( 37) into (34), it's observed that within the range 0 ≤ Ω F 0 ≤ 1, the constraint (34) remains consistently satisfied.Consequently, this model aligns with the requirements of BBN and offers a potential explanation for the acceleration observed in the late-time universe.

Hyperbolic tangent-power model
Following a similar approach to [55,81], we explore the hyperbolic tangent power model represented by: Here, the parameters n and λ > 0 are introduced.The corresponding effective energy density (as denoted in Equation ( 15)) is formulated as: The expression for the parameter n is derived from the first Friedmann equation at the present epoch: Here, Ω F 0 = 1 − Ω m0 − Ω r0 denotes the present value of a quantity.Substituting equations ( 38) and ( 39) into (33) yields: Upon imposing the constraint (34) on the parameter space of the hyperbolic tangent-power model, we obtained Figure 1.Through the plot of T f against the model parameter n, it's observed that the BBN constraints (34) are satisfied within the range n ≤ 1.9.This validation indicates that this model aligns with the requirements of BBN and presents a potential explanation for the acceleration observed in the late-time universe.

VI. DYNAMICAL SYSTEM ANALYSIS
Within this section, we establish a dynamical system that encompasses both the background and perturbed equations associated with a general function F (Q).To facilitate this, we introduce specific dynamical variables, aiming to transform equations ( 18)-( 19) and ( 24) into first-order autonomous systems represented by: The variable u quantifies the expansion of disturbances in matter, while variables x and y are intricately linked to the evolutionary aspects of the universe's background.Consequently, at any given instance, the matter density contrast remains positive.Perturbations in matter are characterized as increasing when u > 0 and decreasing when u < 0. The cosmic background parameters, namely Ω m , Ω Q , and w eff , are represented by the following expressions: .
These parameters allow the cosmological equations to be reformulated into a dynamic system using variables (41) given by: where ( ′ ) denotes differentiation with respect to ln a and (.) stands for differentiation with respect to t.Additionally, Ḣ The physical system comprises a composite space involving the perturbed space denoted as P, housing the variable u, and the background phase space termed as B, encompassing variables x and y.The collective phase space of this combined system, adhering to the physical condition 0 ≤ Ω m ≤ 1, is defined as: It's noteworthy that the projection of orbits originating from the product space Ψ onto space B results in the reduction to corresponding background orbits.The identification and evaluation of the system's critical points play a pivotal role in understanding its dynamical evolution.Specifically concerning matter perturbations, the system demonstrates instability when u > 0, indicating an unbounded expansion of matter disturbances in a physical context.Conversely, a stable point with u < 0 signifies the decay of matter disturbances, portraying the system's asymptotic stability concerning perturbations.Moreover, when the system reaches a stable point with u = 0, it suggests the constancy of matter perturbations.In essence, the expansion of matter perturbations, particularly in unstable or saddle points where u > 0, doesn't imply an everlasting state.This scenario is crucial in understanding the matter epoch of the universe.Notably, these unstable or saddle points must precede a stable late-time attractor characterized by u = 0, signifying the phase of acceleration [52].
To proceed with a detailed analysis, it's imperative to define the function F and subsequently establish the term QF QQ .In doing so, we will delve into the exploration of two specific models previously discussed in section IV.These models have been recognized for their intriguing cosmic phenomenology, which will be elucidated in the subsequent subsections.
The Log-square-root model, represented by the function , yields insights into the system's behavior.Specifically, when evaluating the second derivative of F (Q) with respect to Q (QF QQ ), it follows that Considering the system described by equations ( 42)-( 44), its dynamics are delineated as follows: The state variables evolve according to the differential equations: This dynamical system gives rise to four critical points: * Critical Point (1 − y, y, 0) : The identified point on the curve denotes a solution characterized by an energy density Ω m = 0.This implies an absence of matter contribution to the overall energy density within the universe.Conversely, the energy density Ω Q = 1 signifies the complete dominance of dark energy in the total energy density of the universe.This suggests that the observed accelerated expansion of the universe is exclusively propelled by the repulsive nature of dark energy, with no substantial contribution from matter.This acceleration aligns with an effective equation of state w ef f = −1, akin to the behavior exhibited by a cosmological constant.Additionally, at the perturbation level, u = 0 signifies a constant nature of matter perturbation.The curve, being one-dimensional with one vanishing eigenvalue and other eigenvalues at −3 and −2, exhibits stability consistently across the signature of the non-vanishing eigenvalue.Moreover, the critical point (1 − y, y, 0) is nonhyperbolic due to the presence of one vanishing eigenvalue in the corresponding Jacobian matrix.Using the center manifold theory, the system of equations demonstrates asymptotic stability at this equilibrium point.In this scenario, x acts as the central variable while (y, u) represent the stable variables.The corresponding matrices A and B are characterized by A = 0 and The structure of the center manifold takes the form y = h 1 (x) and u = h 2 (x), with the approximation N comprising two components.
For zeroth approximation: ). Therefore the reduced equation gives us This analysis yields the negative linear aspect applicable for n ∈ R {0}.Consequently, the system of equations ( 46)- (48) demonstrates asymptotic stability at the equilibrium point, consistent with the central manifold theory.
In summary, the discussion characterizes a late universe predominantly influenced by dark energy, both in the background and perturbation levels.The configuration characterized by Ω m = 1 and w ef f = 0 signifies the prevalence of matter dominance on the cosmological scale.However, this description falls short in elucidating the structural evolution at the perturbation level, primarily due to the fluctuation in matter overdensity, which follows a behavior of ρ ∝ a − 3 2 when u = − 3 2 .Analysis of the Jacobian matrix reveals eigenvalues of ( . This critical point exhibits an Unstable Spiral behavior for n > 1 and an Unstable nature for n ≤ 1, shedding light on its dynamical characteristics in cosmological models.
* Critical Point (0, 0, 1): The configuration characterized by background parameters Ω m = 1, Ω Q = 0, and w ef f = 0 signifies a critical solution dominated by matter.The perturbative analysis yields u = 1, indicating a proportional variation of matter overdensity δ with the scale factor, i.e., ρ ∝ a, illustrating an increase concurrent with the universe's expansion.The Jacobian matrix associated with this critical point exhibits eigenvalues of (− Notably, the point (0, 0, 1) consistently presents as a saddle point across all values of n.Trajectories converge toward this point and subsequently deviate from it, converging toward a late-time stable point.This observation suggests that this particular critical point serves as a pivotal choice to elucidate the formation of structures during the matter-dominated epoch, effectively addressing dynamics at both the background and perturbation levels.
In this subsection we consider the hyperbolic tangentpower model As The system ( 42)-( 44) becomes In this case, the system has two curves of critical points * Critical Point (1 − y, y, 0) : The trajectory associated with this point characterizes a solution where the effective dominance of dark energy prevails, denoted by Ω Q = 1, inducing cosmic acceleration with an effective equation of state parameter w ef f = −1, akin to the behavior exhibited by a cosmological constant.This configuration suggests a critical point primarily governed by the geometric component within the framework of the f (Q) model.
Case 1: For y = 0 The curve corresponding to (1 − y, y, 0) collapses into the critical point (1, 0, 0), resulting in a one-dimensional structure characterized by one vanishing eigenvalue alongside remaining eigenvalues of −3 and −2.The stability of this point's curve can be evaluated by examining the signatures of the non-zero eigenvalues, revealing consistent stability.Moreover, the critical point (1, 0, 0) is non-hyperbolic due to the presence of one vanishing eigenvalue within its associated Jacobian matrix.Utilizing the center manifold theory approach, the system of equations is revealed to be asymptotically stable.
In this context, with x as the central variable and y and u as stable variables, the matrices involved are A = 0 and The resultant center manifold takes the form of y = h 1 (x) and u = h 2 (x), and the approximation N comprises two components: For zeroth approximation: ). Therefore the reduced equation gives us The analysis of the negative linear portion concludes that the system of equations ( 50)-( 52) exhibits asymptotic stability at the equilibrium point, aligning with the principles established in the central manifold theory.The transformation from (1 − y, y, 0) to ( 12 , 1 2 , 0) yields a one-dimensional structure characterized by one vanishing eigenvalue and remaining eigenvalues of −2 and −3( 7−4n 2 4−4n 2 ).The stability of this point's curve is consistently maintained for n ∈ R {1}, determined through the center manifold theory technique due to the nonhyperbolic nature of the critical point ( 12 , 1 2 , 0).In this scenario, considering y as the central variable and x and u as stable variables, the associated matri- . The resulting center manifold takes the form of x = h 1 (y) and u = h 2 (y), and the approximation N comprises two components.
For zeroth approximation: ). Therefore the reduced equation gives us The analysis of the negative linear portion concludes that the system of equations ( 50)-( 52) exhibits asymptotic stability at the equilibrium point, aligning with the principles established in the central manifold theory.Case 3: For y = 1 The transition from (1−y, y, 0) to (0, 1, 0) results in a onedimensional structure characterized by a vanishing eigenvalue and remaining eigenvalues of −2 and −3( 5−4n 2−4n ).The stability of this curve is consistently maintained for n ∈ R {1}, established through the application of center manifold theory due to the non-hyperbolic nature of the critical point (0, 1, 0).
In this instance, considering x as the central variable and y and u as stable variables, the associated matrices are A = 0 and B = −2 0 0 −3( 5−4n 2−4n ) . Consequently, the resulting center manifold adopts the form of y = h 1 (x) and u = h 2 (x), while the approximation N encompasses two components.
For zeroth approximation: ). Therefore the reduced equation gives us The analysis of the negative linear portion concludes that the system of equations ( 50)-( 52) exhibits asymptotic stability at the equilibrium point, aligning with the principles established in the central manifold theory.
In conclusion, the assessment utilizing center manifold theory revealed the asymptotic stability of the point (1 − y, y, 0), indicated by eigenvalues 0, −2, − 3 Unlike the curve (1 − y, y, 0), the curve (1 − y, y, −2) lacks the capacity to characterize a late-time dark-energy dominated universe at the perturbation level.However, this point serves as an illustrative depiction of the inflationary epoch of the universe, owing to its saddle-like nature and the presence of a negative w ef f value.
2 ) : This particular point corresponds to a scenario of matter domination at the background level, characterized by Ω m = 1, Ω Q = 0, and w ef f = 0.The parameter u = − 3  2 results in a variation of matter overdensity δ proportional to a − 3 2 .However, this characterization does not adequately describe the mechanism behind structure formation at the perturbation level.Analyzing the eigenvalues of the Jacobian matrix, they are 5  2 , −3(n−2), and −3(n−1).This critical point displays instability for n < 2 and assumes a saddlelike behavior for n ≥ 2. These observations highlight the dynamic nature of the system, portraying instability for certain parameter ranges and a saddle-like characteristic for others concerning the perturbation levels in the context of structure formation.
* Critical Point (0, 0, 1) : The background parameters Ω m = 1, Ω Q = 0, and w ef f = 0 characterize this critical point as a matter-dominated solution.At the perturbation level, u = 1 signifies a variation of matter overdensity δ directly proportional to the scale factor a, escalating with the universe's expansion.In the context of the Hyperbolic tangent power model, adherence to the BBN constraints (34) remains valid under n ≤ 1.9.The critical point (0, 0, 1) consistently assumes a saddlelike nature for any value of n < 2, as discerned from its associated Jacobian matrix featuring eigenvalues − 5 2 , −3(n − 2), and −3(n − 1).Trajectories converge toward this point initially, subsequently diverging as they converge towards a late-time stable point.
This characteristic positions the point as a compelling choice to elucidate the formation of structures during the era of matter dominance, offering insights at both the background and perturbation levels within the cosmological framework.

VII. DISCUSSION AND CONCLUSION
In this investigation, we integrated the dynamical system analysis of both background and perturbed equations, aiming to corroborate our findings through an independent validation process.The flat FLRW metric stands as the most suitable depiction of the universe within our study.Our motivation stemmed from the proven effectiveness of f (Q) gravity-based cosmological models in accurately aligning with observable datasets across various background levels.By transforming both background and perturbed equations into an autonomous system, our focus centered on exploring models from existing literature related to higher-curvature gravity with boundary terms.Specifically, our inquiry delved into the Log-square-root model and the Hyperbolic tangent-power model, seeking to evaluate their implications within the cosmological context.
In the Log-square-root model, four critical points were identified, among which the (1 − y, y, 0) point emerged as asymptotically stable through center manifold theory.This critical point signifies an energy density scenario with Ω m = 0, denoting the absence of matter contribution to the universe's overall energy density.Notably, Ω Q = 1, indicating complete dominance of dark energy in the total energy density.This observation suggests that the observed accelerated expansion of the universe is exclusively propelled by the repulsive effects of dark energy, without substantial matter contribution.Additionally, the universe accelerates with w ef f = −1, akin to the behavior of a cosmological constant, while u = 0 implies a constant matter perturbation.Contrastingly, the curve at the (1 − y, y, −2) critical point portrays a scenario where the effective dark energy component dominates, inducing universe acceleration with w ef f = −1.Despite being dominated by the geometric component of the f (Q) model and featuring matter perturbation decay (u = −2), this critical point assumes a saddle nature.As a result, it fails to characterize a late-time dark-energy dominated universe at the perturbation level.However, owing to its saddle-like nature and negative w ef f value, it provides insights into the inflationary epoch of the universe.Another critical point, (0, 0, − 3 2 ), was identified as an Unstable Spiral for n > 1 and Unstable for n ≤ 1, corresponding to matter supremacy at the background level.Despite this characterization, it fails to account for structure development at the perturbation level due to the matter overdensity's variation as ρ ∝ a − 3 2 when u = − 3 2 .Similarly, the point (0, 0, 1) with eigenvalues exhibits saddle-like behavior across all n values.Trajectories traverse this point before diverging towards a late-time stable point, suggesting its suitability in explaining structure formation during the matter-dominant epoch at both background and perturbation levels.Additionally, phase portraits in Figure 2 and Figure 3 visually illustrate the behaviors of these critical points more comprehensively.
Our investigations reveal distinct manifestations of matter disturbances across various critical points.Notably, similar background critical points exhibit disparate behaviors at the perturbation level, a notable observation we highlight.For instance, among the critical points (0, 0, 1) and (0, 0, 3  2 ), only (0, 0, 1) signifies the appropriate expansion pattern for matter structure, marking a significant divergence in their roles within the context of the matter-dominated era at the background level.The critical point (0, 0, 1) presents an intriguing transition toward a late-time accelerated epoch due to its saddle-like nature.However, while characterizing the accelerated dark-energy dominant epoch later on, the curves of critical points (1 − y, y, 0) and (1 − y, y, −2) exhibit identical behaviors at the background level.It's important to note that only the curve (1−y, y, 0) maintains continuous matter perturbations and stability at the perturbation level, rendering it the sole curve of substantial physical and observational interest.
In the context of the Hyperbolic tangent power model, our analysis yielded four critical points, among which the point (1 − y, y, 0) emerged as asymptotically stable via center manifold theory.The associated eigenvalues 0, −2, − 3 across all y, assumes a saddle-like nature.Despite the decay of matter perturbations (u = −2), this critical point fails to describe a late-time dark-energy dominated universe at the perturbation level, distinguishing it from the curve (1 − y, y, 0).Another critical point, (0, 0, − 3 2 ), characterizes matter domination at the background level, exhibiting unstable behavior for n < 2 and assuming a saddle-like nature for n ≥ 2. Despite depicting a saddle-like nature and featuring a negative w ef f value, this point offers insights into the inflationary era of the universe.Furthermore, the point (0, 0, 1) corresponds to a matter-dominated critical solution with u = 1 at the perturbation level.This configuration illustrates the matter overdensity's variation directly proportional to the scale factor a, escalating with the universe's expansion.Notably, this critical point consistently assumes a saddle-like nature for all n < 2 due to the associated Jacobian matrix's eigenvalues − 5 2 , −3(n−2), and −3(n−1).Thus, this critical point appears to offer the most fitting explanation for the formation of structures during the matter-dominated period across both background and perturbation levels.Furthermore, Figure 4's phase portraits offer a clearer representation of behaviors at these critical points, while Figure 1's plot of δT f T f versus the model parameter n reveals that the BBN constraints (34) remain satisfied within the range n ≤ 1.9.
Perturbations play a decisive role in distinguishing critical points that exhibit similarity at the background level, as observed in Model I. Our analysis leads to the conclusion that only the curve (0, 0, 1) bears substantial physical significance in characterizing the matter-dominated epoch.This critical point is pivotal in delineating the formation of matter perturbations when considering both background and perturbation aspects.In contrast, the curve (1 − y, y, 0) meets the observational criterion by depicting a consistent evolution of matter perturbations.This critical point aligns with late-time dark-energy dominance, marking a significant aspect within the cosmological context.Lastly, our investigation reveals the absence of critical points at infinity within our model, underscoring the finite nature of the critical point configurations in our analysis.
In summary, the comprehensive dynamical analysis spanning both background and perturbation levels effectively reinforces the outcomes derived from observational assessments.This independent validation underscores the potential viability of f (Q) gravity as a compelling alternative to the ΛCDM concordance model.As a closing observation, it's noteworthy that stability and the acceleration phase are attained within the framework of this study without the inclusion of an "unobserved" component, specifically dark energy.
More higher-curvature gravity with boundary term models that adhere to the BBN (Big-Bang Nucleosynthesis) restrictions and satisfy all viability conditions could be included in this research to further it.In addition, other applications of dynamical system analysis as well as stability assessments for such systems may be made.These include taking into account other categories of modification, such as f (Q, T ) or Weyl-type theories.With changes to gravity, other linear and nonlinear interactions could also be investigated.

*
Critical Point (1 − y, y, −2): The critical point denoted by (1 − y, y, −2) shares similarities with the critical point (1 − y, y, 0) in the context of the f (Q) model, both corresponding to solutions characterized by Ω m = 0, Ω Q = 1, and an effective equation of state parameter w ef f = −1, indicative of a dominant geometric component.Specifically, the former point, distinguished by u = −2, signifies the decay of matter disturbances.This critical point exhibits a saddle behavior, manifested by eigenvalues (0, −3, 2).Unlike its counterpart (1 − y, y, 0), the critical point (1 − y, y, −2) does not, at the perturbation level, represent a late-time universe dominated by dark energy.Nonetheless, owing to its saddle nature and the negative w ef f value, this point portrays the inflationary epoch of the universe.* Critical Point (0, 0, − 3 2 ):

2 2+4n−4n 2 +3y−8ny+4n 2 y 1+2n−2n 2 −4ny+2n 2 y 2 2+4n−4n 2
for y ∈ [0, 1].This characterization signifies the prevalence of dark energy dominance in the late-stage universe alongside constant matter perturbation (u = 0).Notably, matter exhibits negligible impact on the repulsive forces responsible for the observable accelerated expansion of the universe.Consequently, this critical point encapsulates the epoch of a late-time-dark-energy dominated universe.*Critical Point(1−y, y, −2) : The trajectory associated with this point represents a solution where the dominance of the effective dark energy component, indicated Critical Points, Stability Conditions, Matter perturbation and EoS parameterby Ω Q = 1, induces cosmic acceleration with an effective equation of state parameter w ef f = −1, akin to behavior reminiscent of a cosmological constant.In the context of the f (Q) model, this critical point is predominantly governed by the geometric component.Furthermore, it is characterized by the decay of matter perturbations, denoted by u = −2.However, despite these traits, this critical point exhibits a saddle-like nature with eigenvalues 0, 2, − 3 +3y−8ny+4n 2 y 1+2n−2n 2 −4ny+2n 2 y for all y.