1 Introduction

Despite the undeniable success of the ‘standard model of cosmology’, governed by general relativity (GR), there are atleast two motivations behind the popularity of the modified theories of gravity in modern research. The first being to alleviate the requirement of dark sector to reason the present acceleration of the Universe. The other is purely theoretical, towards the renormalizability of GR and to extend it to an ultimate theories of quantum gravity. Teleparallel gravity has been extensively researched in this regard in the past few years. To be specific, the Levi-Civita connection, which serves as the foundation of GR, may be replaced by a general affine connection on spatially flat spacetime with vanishing torsion, allowing its non-metricity to assume complete responsibility for defining gravity. ‘Symmetric teleparallel gravity’ is the name given to this particular theory. Long back, in order to unify gravity with electromagnetism, Einstein [1] attributed gravity to the torsion of spacetime in a more developed ‘metric teleparallel theory’, which is based on an affine connection with vanishing curvature and non-metricity both. Under teleparallel theories, one can construct either the so-called torsion scalar \({\mathbb {T}}\) from this torsion tensor in the metric teleparallel branch or the non-metricity scalar Q in its symmetric counterpart. Thereafter, in the Einstein-Hilbert action term of GR, one can replace the Ricci scalar by \({\mathbb {T}}\) in the metric teleparallel theory and by Q in the symmetric teleparallel theory to obtain the respective field equations. However, other than a total divergence term, the theories are actually equivalent to GR and thus also rely on dark components of the Universe as GR does. To address this dark sector issue of the Universe, without inviting a scalar field, natural extensions in terms of \(f({\mathbb {T}} )\) and f(Q)-theories were motivated. While these three theories of gravity (often referred to as the geometrical trinity of gravity) are regarded to be entirely equivalent, these two counterparts of GR nearly went ignored until recently, when the extension in these two theories, respectively \(f({\mathbb {T}})\) and f(Q) theories of gravity in line with the f(R) extension of GR were studied as a possible alternate source of dark energy to beautifully explain the present accelerating Universe [2,3,4,5,6,7,8,9,10,11]. The foremost advantage of these two theories are their second order field equations, as in GR, and unlike in the modified f(R) gravity theory where the field equation is of the fourth order [12]. However, the older and much matured \(f({\mathbb {T}})\) formulation displays certain issues as addressed in an array of publications [13,14,15,16] which the f(Q) theory is free from. On generic FLRW backgrounds, the significant coupling difficulties that may be observed in \(f({\mathbb {T}})\) theories are absent in f(Q) models. Also, the predictions of the f(Q) and \(f({\mathbb {T}})\) models correspond in the small-scale quasi-static limit, but that at higher scales the f(Q) models generically transmit 2 scalar degrees of freedom that are absent in the case of \(f({\mathbb {T}})\). These two degrees of freedom vanish around maximally symmetric backgrounds, resulting in the strong coupling problem that has been explored [17]. In addition, one can investigate the fascinating work [18], in which the formalism of \(f({\mathbb {T}})\) gravity was created in agreement with GR. In a further serious scrutiny of the field equations of \(f({\mathbb {T}})\) theories we can observe an unusual presence of anti-symmetric elements, which is not the case in f(Q) or f(R) or most other theories of gravity. This anti-symmetric part of the field equation of \(f({\mathbb {T}})\) obtained from the variation of the action term with respect to the metric tensor actually equals the field equation obtained from the variation of the same action by the affine connection [19]. As a result, \(f({\mathbb {T}})\) theories (except \(f({\mathbb {T}})={\mathbb {T}}\)) are also not locally Lorentz invariant and possess extra degrees of freedom which remain absent from GR. This is a devastating blow, as the absence of Lorentz invariance necessitates the implementation of a system of 16 equations in \(f({\mathbb {T}})\) theories instead of 10 equations in GR. On a side by side comparison, the \(f({\mathbb {T}})\) connections must comply with \(Q_{\mu \nu \gamma }=0\), which are 40 independent equations due to the symmetry of the non-metricity tensor \(Q_{\mu \nu \gamma }\) in the second and third indices. Whereas, the f(Q) connections must obey the vanishing torsion tensor criteria, \({\mathbb {T}}^\gamma _{\,\,\,\mu \nu }=0\); only 24 independent equations since the \(\mu =\nu \) case is trivially satisfied, making \(f({\mathbb {T}})\) connections much more restrictive than f(Q). Consequently, in a spatially flat FLRW spacetime, there are in total 3 classes of f(Q) formulations possible, yielding 3 different aspects of cosmology under the umbrella of the same f(Q) theory [20]. Whereas, there is only a single class of \(f({\mathbb {T}})\) formulation in that spacetime, and this class is predominantly equivalent to the one obtained using the diagonal tetrad in the Cartesian coordinates [21]. For an in depth comparison of these two theories, one can further look into [22,23,24,25,26,27]. However, as mentioned earlier, the modified f(Q) theory under the symmetric teleparallelism is merely at its infancy; a lot of theoretical investigation is still due before we can vouch for its robustness. Specially, the newly discovered f(Q) constructions arising from the non-vanishing classes of affine connections in both isotropic and homogeneous spacetimes [20] and static spherically symmetric spacetimes [28] are yet to be tested.

On the other hand, the ability to formulate practical guidelines for cosmological applications gives birth to the relationship between theory and observational datasets. For a certain class of f(Q) gravity theories under symmetric teleparallelism, the Hubble parameter may be computed analytically by simply adhering to some partial extensions of generally accepted conventional approaches. In this regard, we should mention some of the introductory data analysis works conducted in f(Q) theory. In [29], f(Q) theory challenged the \(\Lambda \)CDM model for the first time, in the sense that the new gravity theory despite having the same number of free parameters as in \(\Lambda \)CDM, at a cosmological framework it possibly can avoid \(\Lambda \)CDM as a limit and thus can alleviate the cosmological constant problem. Moreover, contact with observations at both background and perturbation levels reveals that the model, in some datasets is slightly preferred than \(\Lambda \)CDM cosmology, although in all cases the two models are statistically indiscriminate. The f(Q) model also does not exhibit early dark energy features, and thus it immediately passes BBN constraints. Within few weeks of this work, support in the form of [30] appeared. Modern cosmology probes are used in the statistical analysis, and a certain matter-energy composition is assumed, such that the candidate models provide promising modified gravity candidates to represent the cosmic backdrop. Starting with an expression for the pressure-energy density ratio, the so-called equation of state (EoS) parameter \(\omega \), the technique moves on to the derivation and analysis of other cosmic parameters corresponding to two particular f(Q) models, namely, \(f(Q)=-Q+ \frac{\alpha }{Q}\) and \(f(Q)= -\alpha Q^n\), both of which demonstrate excellent fit with the cosmological data. The model parameters \(\alpha \) and n can be easily adjusted to retrieve GR, providing the opportunity for a clean comparison.

The present paper is organised as follows:

After the brief introductory discussion presented above, in Sect. 2 the fundamental mathematical formulation of the modified f(Q) theories is presented, followed by the (effective) pressure and energy equations in the spatially flat Friedmann–Lemaître–Robertson–Walker (FLRW) geometric background in Sect. 3. We consider a form of the EoS parameter \(\omega (z)\) as a function of red-shift parameter z in Sect. 4. Under this section, two separate f(Q) models are analysed in two separate subsections and the required cosmological parameters are derived. Section 4.1 assumes \(f(Q)=-Q+\frac{\alpha }{Q}\) and Sect. 4.2 assumes the model \(f(Q)=-\alpha Q^n\). Section 5 delivers all the relevant observational data analysis of the above-mentioned models. A detailed concluding remarks on all the acquired findings in Sect. 7 is preceded by Sect. 6 discussing briefly the behavior of the cosmological parameters.

2 Basic concepts of f(Q) gravity

In this section We discuss the detailed formulation of the symmetric teleparallelism, specially its extension the modified f(Q) theory. We begin with a 4-dimensional Lorentzian manifold \(M^{4}\), a line element governed by the metric tensor \(g_{\mu \nu }\) in certain coordinate system \( \{x^{0},x^{1},x^{2},x^{3}\}\) and a non-tensorial term, the affine connection \(\Gamma _{\,\,\,\mu \nu }^{\alpha }\), defining the covariant derivative \( \nabla \) and also taking care of the three main aspects of the spacetime geometry corresponding to this connection, the curvature, torsion and non-metricity. However, once we restrict ourselves to specifically consider vanishing of both the non-metricity and the torsion tensors corresponsing to the connection, we can assert that there is only a unique connection available, the Levi-Civita connection \(\mathring{\Gamma }\) and it has a well-known relation with the metric g given by

$$\begin{aligned} \mathring{\Gamma }_{\,\,\,\mu \nu }^{\alpha }=\frac{1}{2}g^{\alpha \beta }\left( \partial _{\nu }g_{\beta \mu }+\partial _{\mu }g_{\beta \nu }-\partial _{\beta }g_{\mu \nu }\right) . \end{aligned}$$
(1)

So, the Levi-Civita connection is basically a function of the metric g and not an independent player in the spacetime geometry. Here, instead we consider a torsion-free affine connection \(\Gamma \) on a flat spacetime which is not metric-compatible, the incompatibility is characterised by the non-metricity tensor

$$\begin{aligned} Q_{\lambda \mu \nu }=\nabla _{\lambda }g_{\mu \nu }=\partial _{\lambda }g_{\mu \nu }-\Gamma _{\,\,\,\lambda \mu }^{\beta }g_{\beta \nu }-\Gamma _{\,\,\,\lambda \nu }^{\beta }g_{\beta \mu }\ne 0, \end{aligned}$$
(2)

We can always express

$$\begin{aligned} \Gamma ^{\lambda }{}_{\mu \nu }=\mathring{\Gamma }^{\lambda }{}_{\mu \nu }+L^{\lambda }{}_{\mu \nu } \end{aligned}$$
(3)

where \(L^{\lambda }{}_{\mu \nu }\) is the disformation tensor. It follows that

$$\begin{aligned} L^{\lambda }{}_{\mu \nu }=\frac{1}{2}(Q^{\lambda }{}_{\mu \nu }-Q_{\mu }{}^{\lambda }{}_{\nu }-Q_{\nu }{}^{\lambda }{}_{\mu }). \end{aligned}$$
(4)

We can construct two different types of non-metricity vectors,

$$\begin{aligned} Q_{\mu }=g^{\nu \lambda }Q_{\mu \nu \lambda }=Q_{\mu }{}^{\nu }{}_{\nu },\qquad {\tilde{Q}}_{\mu }=g^{\nu \lambda }Q_{\nu \mu \lambda }=Q_{\nu \mu }{}^{\nu }.\nonumber \\ \end{aligned}$$
(5)

The non-metricity conjugate tensor \(P^{\lambda }{}_{\mu \nu }\) is given by

$$\begin{aligned}{} & {} P^{\lambda }{}_{\mu \nu }=\frac{1}{4}\nonumber \\{} & {} \quad \times \left( -2L^{\lambda }{}_{\mu \nu }+Q^{\lambda }g_{\mu \nu }-{\tilde{Q}}^{\lambda }g_{\mu \nu }-\frac{1}{2} \delta _{\mu }^{\lambda }Q_{\nu }-\frac{1}{2}\delta _{\nu }^{\lambda }Q_{\mu }\right) . \nonumber \\ \end{aligned}$$
(6)

Finally, one can define the non-metricity scalar

$$\begin{aligned} Q=-Q_{\lambda \mu \nu }P^{\lambda \mu \nu }=\frac{1}{2}Q_{\lambda \mu \nu }L^{\lambda \mu \nu }+\frac{1}{2}Q_{\lambda }{\tilde{Q}}^{\lambda }-\frac{1}{4} Q_{\lambda }Q^{\lambda }. \nonumber \\ \end{aligned}$$
(7)

By varying the action term

$$\begin{aligned} S=\int \left[ \frac{1}{2\kappa }f(Q)+{\mathcal {L}}_{M}\right] \sqrt{-g}\,d^{4}x \end{aligned}$$
(8)

with respect to the metric tensor \(g^{\mu \nu }\), we obtain the field equation [31]

$$\begin{aligned}{} & {} \frac{2}{\sqrt{-g}}\nabla _{\lambda }(\sqrt{-g}f_{Q}P^{\lambda }{}_{\mu \nu })+\frac{1}{2}fg_{\mu \nu }\nonumber \\{} & {} \quad +f_{Q}(P_{\nu \rho \sigma }Q_{\mu }{}^{\rho \sigma }-2P_{\rho \sigma \mu }Q^{\rho \sigma }{}_{\nu })=-\kappa T_{\mu \nu }. \end{aligned}$$
(9)

\(T_{\mu \nu }\) is the energy-momentum tensor generated from the matter Lagrangian \({\mathcal {L}}_{M}\). We assume a barotropic perfect fluid given by

$$\begin{aligned} T_{\mu \nu }=(p+\rho )u_{\mu }u_{\nu }+pg_{\mu \nu }, \end{aligned}$$

with isotropic pressure p, energy-density \(\rho \) and the four-velocity vector \(u^{\mu }\). Recently, the covariant representation of the field equation (9) was derived [32]

$$\begin{aligned} f_{Q}\mathring{G}_{\mu \nu }+\frac{1}{2}g_{\mu \nu }(f-Qf_{Q})+2f_{QQ} \mathring{\nabla }_{\lambda }QP^{\lambda }{}_{\mu \nu }=-\kappa T_{\mu \nu }, \nonumber \\ \end{aligned}$$
(10)

where \(\mathring{G}_{\mu \nu }=\mathring{R}_{\mu \nu }-\frac{1}{2}g_{\mu \nu }\mathring{R}\), is the Einstein tensor corresponding to the Levi-Civita connection. The field equations (10) can be equivalently written in the effective form

$$\begin{aligned} \mathring{G}_{\mu \nu }=-\frac{\kappa }{f_{Q}}T_{\mu \nu }^{eff}=-\frac{ \kappa }{f_{Q}}T_{\mu \nu }+T_{\mu \nu }^{DE}, \end{aligned}$$
(11)

where \(\frac{-\kappa }{f_{Q}}\) is the effective gravitational constant and for its positivity in our construction, we assume \(f_{Q}<0\). The dark energy component emerged from this modification of STEGR is given by

$$\begin{aligned} T_{\mu \nu }^{DE}=\frac{-1}{f_{Q}}\left[ \frac{1}{2}g_{\mu \nu }(f-Qf_{Q})+2f_{QQ} \mathring{\nabla }_{\lambda }QP^{\lambda }{}_{\mu \nu }\right] . \end{aligned}$$
(12)

On the other hand, since the affine connection is an independent entity in the symmetric teleparallel theory, we vary the action (8) with regard to the affine connection \(\Gamma \) to obtain the connection field equation as [32]

$$\begin{aligned} \nabla _{\mu }\nabla _{\nu }(\sqrt{-g}f_{Q}P^{\mu \nu }{}_{\gamma })=0, \end{aligned}$$
(13)

on the basis of the assumption that the matter Lagrangian \({\mathcal {L}}_{M}\) is not a function of the affine connection. Moreover, it has been showed that [19] this second field equation (13) is trivially satisfied in a model-independent manner in the spacetime geometry we are going to consider in the Sect. 3. So in the present article our sole attention is devoted to the metric field equation (10). It should be mentioned that \(R_{\sigma \mu \nu }^{\rho }=0\) is one of the restrictions that was utilised when developing the f(Q)-theory. This indicates that there is a unique coordinate system that can be chosen to make the affine connection disappears, denoted by the expression \(\Gamma _{\mu \nu }^{\lambda }=0\). The term “coincident gauge” refers to this kind of circumstance. For more details of the f(Q) theory and its cosmological applications, one can see [29, 33,34,35,36,37,38,39,40,41,42] and the references therein.

3 f(Q) cosmology in isotropic and homogeneous Universe

The “cosmological principle” states that on a large enough scale our Universe is homogenous and isotropic, that is, it is the same at every point and in every direction. Based on this, the most reasonable and theoretically and observationally supported model of the present Universe is the spatially flat Friedmann–Lemaitree–Robertson–Walker (FLRW) spacetime given by the line element in Cartesian coordinates

$$\begin{aligned} ds^{2}=-dt^{2}+a^{2}(t)[dx^{2}+dy^{2}+dz^{2}], \end{aligned}$$
(14)

where a(t) is the scale factor of the Universe. We proceed with the coincident gauge choice as discussed above and obtain the non-metricity scalar as \(Q=6H^{2}\), where \(H=\frac{\overset{.}{a}}{a}\) is the Hubble parameter, which measures the expansion rate of the Universe, and \(\dot{()}\) indicates a derivative with regard to cosmic time t.

In this context, the field equations (10) give the following expressions of the energy density \(\rho \) and the isotropic pressure p [32, 37]

$$\begin{aligned}{} & {} \rho =\frac{f}{2}-6H^{2}f_{Q}, \end{aligned}$$
(15)
$$\begin{aligned}{} & {} p=\left( \overset{.}{H}+\frac{\overset{.}{f_{Q}}}{f_{Q}}H\right) \left( 2f_{Q}\right) - \left( \frac{f}{2}-6H^{2}f_{Q}\right) , \end{aligned}$$
(16)

Furthermore, to explain cosmic history and the possible transition to an accelerated phase, we use the equation of state (EoS) parameter \(\omega \), given by

$$\begin{aligned} \omega =\frac{p}{\rho }=-1+\frac{\left( \overset{.}{H}+\frac{\overset{.}{ f_{Q}}}{f_{Q}}H\right) \left( 2f_{Q}\right) }{\left( \frac{f}{2} -6H^{2}f_{Q}\right) }. \end{aligned}$$
(17)

On the other hand, using (11) the effective energy density \( \rho _{eff}\) and effective pressure \(p_{eff}\) of the cosmic fluid can be written as [37]

$$\begin{aligned} 3H^{2}&=-\frac{1}{f_{Q}}\rho _{eff}\,=-\frac{1}{f_{Q}}\left( \rho +\frac{ 6H^{2}f_{Q}-f}{2}\right) , \end{aligned}$$
(18)
$$\begin{aligned} -(2{\dot{H}}+3H^{2})&=-\frac{1}{f_{Q}}p_{eff}\,\nonumber \\&=-\frac{1}{f_{Q}}\left( p+\frac{ f-6H^{2}f_{Q}-4\overset{.}{f}_{Q}H}{2}\right) . \end{aligned}$$
(19)

Further, in the limiting situation \(f\left( Q\right) =-Q=-6H^{2}\), the gravitational action (8) is reduced to the standard Hilbert-Einstein form and Eqs. (18) and (19) reduce to the standard Friedmann equations of GR, \(3H^{2}=\rho \), and \(2{\dot{H}} +3H^{2}=-p\), respectively. Thus, the effective EoS parameter \(\omega _{eff}\) is

$$\begin{aligned} \omega _{eff}=\frac{p_{eff}}{\rho _{eff}}=\frac{2p+f-6H^{2}f_{Q}-4\overset{.}{f}_{Q}H}{2\rho +6H^{2}f_{Q}-f}. \end{aligned}$$
(20)

To describe the accelerated/decelerated aspect of the cosmic expansion, we consider the deceleration parameter q, expressed as,

$$\begin{aligned} q=-\frac{\overset{.}{H}}{H^{2}}-1, \end{aligned}$$
(21)

4 Cosmological models with specific form of EoS

In this section, we examine dark energy parametrization, which displays quintessence behavior with the evolution of the cosmos. Our principal objective is to study this parametrization using existing cosmological data. For simplicity, we will adopt the red-shift as the independent variable, expressed as \(z=\frac{a_{0}}{a}-1\), with the current scale factor \(a_{0}\) fixed to 1. In general, there is no theoretical method for selecting the optimal \(\omega \left( z\right) \), but by utilizing observational data, suitable parametrizations can be found. In literature, several EoS dark energy parametrization models were proposed and fitted with observational data. Ref. [43] proposed an one-parameter family of EoS dark energy model. Two-parameters family of EoS dark energy parametrizations, especially the Chevallier-Polarski-Linder parametrization [44, 45], the Linear parametrization [45,46,47,48], the Logarithmic parametrization [49], the Jassal-Bagla-Padmanabhan parametrization [50], and the Barboza-Alcaniz parametrization [51], were also investigated. Further, in [52,53,54] three and four parameters family of EoS dark energy parametrizations are examined. Here, we assume that the EoS parameter is parametrized as a function of red-shift z,

$$\begin{aligned} \omega \left( z\right) =-\frac{1}{1+3\beta \left( 1+z\right) ^{3}}, \end{aligned}$$
(22)

where \(\beta \) is the only free parameter. The purpose of choosing this parametrization for \(\omega \left( z\right) \) is that for very large red-shift \(z\gg 1\), i.e., at the early phases of cosmological evolution, \( \omega \) is approximately zero, representing the behavior of the EoS parameter for a pressureless fluid (ordinary matter), but it gradually decreases to negative values at the present i.e. \(z=0\), leads to negative pressure and \(\omega =-\frac{1}{1+3\beta }\). This last equation clearly shows that for an accelerated Universe scenario: \(\omega \le -\frac{1}{3}\), resulting in a constraint of the value of \(\beta \) as \(\beta \le \frac{2}{3} \). In addition, we can see that at later times \(z\rightarrow -1\), \(\omega \) tends to \(-1\), which is similar to the behavior of the cosmological constant \(\Lambda \). Also, for \(\beta =0\), \(\omega \) reduces to \(-1\). These situations are summarized below,

  • \(\omega \longrightarrow -1\), as \(z\longrightarrow -1\),

  • \(\omega \longrightarrow 0\), for \(z\gg 1\),

  • \(\omega =-\frac{1}{1+3\beta }\), for \(z=0\).

To make comparisons of theoretical results with cosmological data simple, we utilize the red-shift z instead of the usual time variable t. As a result, we can change the derivatives with regard to time with the derivatives with regard to red-shift using the relationship,

$$\begin{aligned} \frac{d}{dt}=\frac{dz}{dt}\frac{d}{dz}=-(1+z)H(z)\frac{d}{dz}. \end{aligned}$$
(23)

The deceleration parameter q can be calculated as a function of cosmic red-shift,

$$\begin{aligned} q\left( z\right) =\left( 1+z\right) \frac{1}{H\left( z\right) }\frac{ dH\left( z\right) }{dz}-1. \end{aligned}$$
(24)

Also, the derivative of the Hubble parameter can be expressed as,

$$\begin{aligned} {\dot{H}}=-(1+z)H(z)\frac{dH}{dz}. \end{aligned}$$
(25)

Now, we examine several specific cosmological models in f(Q) gravity theory, models that correspond to different classes of the function f(Q). We also analyze the behavior of geometric and physical cosmological parameters in f(Q) gravity such as energy density, pressure, and deceleration parameter.

4.1 \(f\left( Q\right) =-Q+\frac{\alpha }{Q}\)

For a first case of a cosmological model in f(Q) gravity, consider the scenario where the function f(Q) can be expressed as, \(f(Q)=-Q+\frac{ \alpha }{Q}\), where \(\alpha \) is a constant. So, we get \(f_{Q}=-1-\frac{ \alpha }{Q^{2}}\) and \(f_{QQ}=\frac{2\alpha }{Q^{3}}\). The Friedmann equations (15) and (16) for this particular f(Q) model, reduce to

$$\begin{aligned} \rho =\frac{\alpha }{4H^{2}}+3H^{2}, \end{aligned}$$
(26)

and

$$\begin{aligned} p=\frac{\alpha \overset{.}{H}}{6H^{4}}-2\overset{.}{H}-\frac{\alpha }{4H^{2}} -3H^{2}. \end{aligned}$$
(27)

Using (17), we obtain the EoS parameter in terms of Hubble parameter as

$$\begin{aligned} \omega =\frac{2\overset{.}{H}}{3H^{2}}-\frac{16\overset{.}{H}H^{2}}{\alpha +12H^{4}}-1. \end{aligned}$$
(28)

The differential equation for \(H\left( z\right) \) is obtained by (28) the presumed ansatz of \(\omega \) as indicated in (22),

$$\begin{aligned} \frac{2\overset{.}{H}\left( \alpha -12H^{4}\right) }{3H^{2}\left( \alpha +12H^{4}\right) }-\frac{3\beta (z+1)^{3}}{3\beta (z+1)^{3}+1}=0. \end{aligned}$$
(29)

Therefore, the solution obtained for the Hubble parameter \(H\left( z\right) \) as a function of red-shift z is

$$\begin{aligned}{} & {} \!\!\! H\left( z\right) \nonumber \\{} & {} \ =\frac{\sqrt{\sqrt{\frac{\left( \alpha {+}12H_{0}^{4}\right) ^{2}\left( 3\beta (z+1)^{3}+1\right) ^{2}}{(3\beta {+}1)^{2}H_{0}^{4}} -48\alpha }{+}\frac{\left( \alpha {+}12H_{0}^{4}\right) \left( 3\beta (z+1)^{3}+1\right) }{(3\beta {+}1)H_{0}^{2}}}}{2\sqrt{6}}, \nonumber \\ \end{aligned}$$
(30)

where \(H_{0}\) represents the current value (i.e. at \(z=0\)) of the Hubble parameter.

For this specific case, the expression of the deceleration parameter obtained by including (30) into (24) as

$$\begin{aligned} q\left( z\right) =\frac{9\beta (z+1)^{3}\left( \alpha +12H_{0}^{4}\right) }{ 2(3\beta {+}1)H_{0}^{2}\sqrt{\frac{\left( \alpha {+}12H_{0}^{4}\right) ^{2}\left( 3\beta (z+1)^{3}{+}1\right) ^{2}}{(3\beta {+}1)^{2}H_{0}^{4}} {-}48\alpha }}-1. \nonumber \\ \end{aligned}$$
(31)

4.2 \(f\left( Q\right) =-\alpha Q^{n}\)

For a second case of a cosmological model in f(Q) gravity, we assume the scenario where the function f(Q) can be expressed as a power-law form, \( f(Q)=-\alpha Q^{n}\), where \(\alpha \) and n are model parameters. So, we get \(f_{Q}=-\alpha nQ^{n-1}\) and \(f_{QQ}=-\left( n-1\right) \alpha nQ^{n-2}\). The Friedmann equations (15) and (16) for this particular f(Q) model, reduce to

$$\begin{aligned} \rho =\alpha 2^{n-1} 3^n (2 n-1) \left( H^2\right) ^n, \end{aligned}$$
(32)

and

$$\begin{aligned} p=\alpha \left( -6^{n-1}\right) (2 n-1) \left( H^2\right) ^{n-1} \left( 2\overset{.}{H} n+3 H^2\right) . \end{aligned}$$
(33)

Using (17), we obtain the EoS parameter in terms of Hubble parameter as

$$\begin{aligned} \omega =-\frac{2\overset{.}{H}n}{3H^{2}}-1. \end{aligned}$$
(34)

The differential equation for \(H\left( z\right) \) is obtained by Eq. (34) the presumed ansatz of \(\omega \) as indicated in Eq. (22),

$$\begin{aligned} \frac{2\overset{.}{H}n}{3H^{2}}+\frac{3\beta (1+z)^{3}}{3\beta (1+z)^{3}+1}=0. \end{aligned}$$
(35)

Therefore, the solution obtained for the Hubble parameter \(H\left( z\right) \) as a function of red-shift z is,

$$\begin{aligned} H\left( z\right) =H_{0}\left[ \frac{3\beta (1+z)^{3}+1}{3\beta +1}\right] ^{ \frac{1}{2n}}, \end{aligned}$$
(36)

where \(H_{0}\) represents the current value (i.e. at \(z=0\)) of the Hubble parameter. The previous equation can be rewritten as,

$$\begin{aligned} H\left( z\right) =H_{0}\left[ \left( 1-\gamma \right) (1+z)^{3}+\gamma \right] ^{\frac{1}{2n}}, \end{aligned}$$
(37)

For reasons of simplicity, we have introduced \(\gamma =\frac{1}{1+3\beta }\). It should be observed that the standard \(\Lambda \)CDM model is equivalent to the scenario \(n=1\), with the current cold dark matter density parameter \( \Omega _{m}^{0}=(1-\gamma ) \). As a result, the model parameter n is an excellent indicator of the current model’s deviation from the \(\Lambda \)CDM model.

For this specific case, the expression of the deceleration parameter obtained by including Eq. (36) into Eq. (24) as,

$$\begin{aligned} q\left( z\right) =\frac{9\beta (1+z)^{3}}{6\beta n(1+z)^{3}+2n}-1. \end{aligned}$$
(38)

In the next section, we will attempt to estimate the values of \(H_{0}\), n, \(\alpha \), and \(\beta \) using Hubble, SNe Ia, and BAO datasets. The behavior of cosmological parameters such as density, pressure, and deceleration parameter can be examined for each model (Models 1 and 2) using the values of \(H_{0}\), n, \(\alpha \), and \(\beta \).

5 Observational constraints from Hubble, SNe Ia and BAO

This section is involved with the different observational datasets used to restrict the parameters \(\alpha \), \(\beta \), n and \(H_{0}\). To get the posterior distributions of the parameters, we use the usual Bayesian algorithms and a Markov Chain Monte Carlo (MCMC) approach with the emcee python package [55]. This stimulation is accomplished by the use of Hubble measurements (i.e., Hubble datasets), Type Ia supernovae (SNe Ia) datasets, and BAO datasets. The probability function is used to maximize the best fits of the parameters

$$\begin{aligned} {\mathcal {L}}\propto exp(-\chi ^{2}/2), \end{aligned}$$
(39)

where \(\chi ^{2}\) is the chi-square function. The \(\chi ^{2}\) functions are explained below for different datasets.

5.1 Hubble datasets

The Hubble results are the first observational data sample used during our research. The Hubble parameter is written as \(H(z)=-dz/[dt(1+z)]\). Given that dz is determined via a spectroscopic survey, the model-independent value of the H(z) can be estimated by measuring dt. It is commonly known that the Hubble parameter can directly estimate the rate of cosmological expansion. In principle, two methods for calculating the Hubble parameter at various red-shifts are widely used: the differential ages \(\Delta t\) of galaxies and the line of sight BAO technique. In this work, we restrict the model using a collection of 57 Hubble parameter data points in the red-shift range \(0.07\le z\le 2.41\) published by Sharov and Vasiliev [56]. For Hubble datasets, the \(\chi _{Hubble}^{2}\) function is given respectively for the two models

$$\begin{aligned}{} & {} \chi _{Hubble\_1}^{2}\left( H_{0},\alpha ,\beta \right) \nonumber \\{} & {} \quad =\sum \limits _{k=1}^{57}\frac{[H_{th}(z_{k},H_{0},\alpha ,\beta )-H_{obs}(z_{k})]^{2}}{\sigma _{H(z_{k})}^{2}}, \end{aligned}$$
(40)

and

$$\begin{aligned}{} & {} \chi _{Hubble\_2}^{2}\left( H_{0},n,\beta \right) \nonumber \\{} & {} \quad =\sum \limits _{k=1}^{57} \frac{[H_{th}(z_{k},H_{0},n,\beta )-H_{obs}(z_{k})]^{2}}{\sigma _{H(z_{k})}^{2}}. \end{aligned}$$
(41)

Here, \(H_{obs}\) is the Hubble parameter value recovered from cosmic observations, \(H_{th}\) is its theoretical value estimated, and \(\sigma _{H}\) is the standard deviation in the observed value of \(H\left( z\right) \). The \( 1-\sigma \) and \(2-\sigma \) contour graphs in Figs. 1 and 2 show the best fit values for the parameters of both models obtained from the Hubble datasets. The likelihoods are extremely closely adapted to Gaussian distributions. The best-fit values of the model parameters found are: \( H_{0}=69.3_{-2.2}^{+2.4}\), \(\alpha =3.9_{-4.0}^{+5.2}\), and \(\beta =0.129_{-0.024}^{+0.027}\) for the first model, and \(H_{0}=65.5_{-4.6}^{+4.5}\), \(n=1.16_{-0.15}^{+0.14}\), and \(\beta =0.33_{-0.24}^{+0.39}\) for the second model. Also, Figs. 3 and 4 show the error bar plots for both models and the \(\Lambda \)CDM, with the cosmological constant density parameter \(\Omega _{\Lambda }^{0}=0.7\), the matter density parameter \(\Omega _{m}^{0}=0.3\) and \(H_{0}=69\) km/s/Mpc. It is shown that the f(Q) models completely suit the observational data while deviating somewhat from the \(\Lambda \)CDM.

Fig. 1
figure 1

Constraints on the model parameters at \(1-\sigma \) and \(2- \sigma \) confidence interval using the Hubble datasets (Model 1)

Fig. 2
figure 2

Constraints on the model parameters at \(1-\sigma \) and \(2- \sigma \) confidence interval using the Hubble datasets (Model 2)

Fig. 3
figure 3

The plot of H(z) vs the red-shift z for our f(Q) model 1, which is shown in red, and \(\Lambda \)CDM, which is shown in black dashed lines, shows an excellent match to the 57 points of the Hubble datasets

Fig. 4
figure 4

The plot of H(z) vs the red-shift z for our f(Q) model 2, which is shown in red, and \(\Lambda \)CDM, which is shown in black dashed lines, shows an excellent match to the 57 points of the Hubble datasets

Fig. 5
figure 5

Constraints on the model parameters at \(1-\sigma \) and \(2- \sigma \) confidence interval using the SNe Ia datasets (Model 1)

Fig. 6
figure 6

Constraints on the model parameters at \(1-\sigma \) and \(2- \sigma \) confidence interval using the BAO datasets (Model 1)

5.2 SNe Ia datasets

Type Ia Supernovae (SNe Ia) are a strong distance indicator that can be employed to investigate the background expansion of the cosmos. To restrict the aforementioned parameters, we use the latest Pantheon SNe Ia collection, which consists of 1048 SNe Ia data points collected from several SNe Ia samples in the red-shift range \(z\in [0.01,2.3]\) such as SDSS, SNLS, Pan-STARRS1, low-red-shift survey, and HST surveys [57]. For SNe Ia datasets, the \(\chi _{SNe}^{2}\) function is given as

$$\begin{aligned} \chi _{SNe\_1}^{2}=\sum _{i,j=1}^{1048}\Delta \mu _{i}\left( C_{SNe}^{-1}\right) _{ij}\Delta \mu _{j}, \end{aligned}$$
(42)

where \(C_{SNe}\) represents the covariance matrix [57], and

$$\begin{aligned} \quad \Delta \mu _{i}=\mu ^{th}(z_{i},H_{0},\alpha ,\beta )-\mu _{i}^{obs}, \end{aligned}$$
(43)

is the difference between the observable distance modulus value from astronomical data and the theoretical values calculated from the model with the given parameter \(\alpha \), \(\beta \), n and \(H_{0}\). Furthermore, the distance modulus is derived as, \(\mu =m_{B}-M_{B}\), where \(m_{B}\) and \(M_{B}\) signify the measured apparent magnitude and absolute magnitude at a given red-shift z (Trying to retrieve the nuisance parameter using the new BEAMS with Bias Correction technique (BBC) [58]). Its theoretical value is also given by

$$\begin{aligned} \mu (z)=5log_{10}\left[ \frac{D_{L}(z)}{1Mpc}\right] +25, \end{aligned}$$
(44)

where

$$\begin{aligned} D_{L}(z)=c(1+z)\int _{0}^{z}\frac{dz^{\prime }}{H(z^{\prime },H_{0},\alpha ,\beta )}. \end{aligned}$$
(45)

The \(1-\sigma \) and \(2-\sigma \) contour graphs in Fig. 5 show the best-fit values for the parameters of model 1 obtained from the SNe Ia datasets. The best-fit values of the model parameters found are: \( H_{0}=67.9_{-2.6}^{+2.5}\), \(\alpha =4.0_{-4.1}^{+5.3}\), and \(\beta =0.160_{-0.074}^{+0.095}\). Furthermore, Fig. 7 shows the error bar plots for the model and the \(\Lambda \)CDM, with the cosmological constant density parameter \(\Omega _{\Lambda }^{0}=0.7\), the matter density parameter \(\Omega _{m}^{0}=0.3\) and \(H_{0}=69\) km/s/Mpc. The graphic also displays the SNe Ia findings, 1048 data points, with their errors, allowing for a direct comparison of the two models.

5.3 BAO datasets

The last restrictions in this study are obtained by BAO observation. BAO studies oscillations in the early cosmos generated by cosmic perturbations in a fluid composed of photons, baryons, and dark matter and connected via Thompson scattering. BAO observations include the Sloan Digital Sky Survey (SDSS), the Six Degree Field Galaxy Survey (6dFGS), and the Baryon Oscillation Spectroscopy Survey (BOSS) [59, 60]. The equations employed in BAO analysis are (Table 1)

$$\begin{aligned} d_{A}(z)=c\int _{0}^{z}\frac{dz^{\prime }}{H(z^{\prime })}, \end{aligned}$$
(46)
$$\begin{aligned} D_{V}(z)=\left[ \frac{d_{A}^2(z) c z}{H(z)} \right] ^\frac{1}{3}, \end{aligned}$$
(47)

and

$$\begin{aligned} \chi _{BAO\_1}^{2}=X^{T}C_{BAO}^{-1}X, \end{aligned}$$
(48)
$$\begin{aligned} X=\left( \begin{array}{c} \frac{d_{A}(z_{\star })}{D_{V}(0.106)}-30.95 \\ \frac{d_{A}(z_{\star })}{D_{V}(0.2)}-17.55 \\ \frac{d_{A}(z_{\star })}{D_{V}(0.35)}-10.11 \\ \frac{d_{A}(z_{\star })}{D_{V}(0.44)}-8.44 \\ \frac{d_{A}(z_{\star })}{D_{V}(0.6)}-6.69 \\ \frac{d_{A}(z_{\star })}{D_{V}(0.73)}-5.45\\ \end{array} \right) , \end{aligned}$$
Table 1 Values of \(d_{A}(z_{*})/D_{V}(z_{BAO})\) for various values of \(z_{BAO}\) [61]

where \(d_{A}(z)\) represents the angular diameter distance, \( D_{V}(z)\) represents the dilation scale, and \(C_{BAO}\) represents the covariance matrix defined as [61],

$$\begin{aligned} C_{BAO}^{-1}=\left( \begin{array}{cccccc} 0.48435 &{}\quad -0.101383 &{}\quad -0.164945 &{}\quad -0.0305703 &{}\quad -0.097874 &{}\quad -0.106738 \\ -0.101383 &{}\quad 3.2882 &{}\quad -2.45497 &{}\quad -0.0787898 &{}\quad -0.252254 &{}\quad -0.2751 \\ -0.164945 &{}\quad -2.454987 &{}\quad 9.55916 &{}\quad -0.128187 &{}\quad -0.410404 &{}\quad -0.447574 \\ -0.0305703 &{}\quad -0.0787898 &{}\quad -0.128187 &{}\quad 2.78728 &{}\quad -2.75632 &{}\quad 1.16437 \\ -0.097874 &{}\quad -0.252254 &{}\quad -0.410404 &{}\quad -2.75632 &{}\quad 14.9245 &{}\quad -7.32441 \\ -0.106738 &{}\quad -0.2751 &{}\quad -0.447574 &{}\quad 1.16437 &{}\quad -7.32441 &{}\quad 14.5022 \end{array} \right) . \end{aligned}$$

The \(1-\sigma \) and \(2-\sigma \) contour graphs in Fig. 6 shows the best-fit values for the parameters of model 2 obtained from the BAO datasets. The model parameter constraints are derived by minimizing the associated \(\chi ^{2}\) using MCMC and the emcee library. The best-fit values of the model parameters found are: \(H_{0}=69_{-9}^{+10}\), \(\alpha =3.7_{-3.9}^{+5.2}\), and \(\beta =0.148_{-0.035}^{+0.037}\) (Fig. 7).

Fig. 7
figure 7

The plot of \(\mu (z)\) vs the red-shift z for our f(Q) model 1, which is shown in red, and \(\Lambda \)CDM, which is shown in black dashed lines, shows an excellent match to the 1048 points of the Pantheon datasets

6 Behavior of cosmological parameters

The behavior of cosmological parameters is determined by the underlying cosmological model. Several parameters are used to characterize the current state and evolution of the Universe in the classic \(\Lambda \)CDM model, which defines the Universe as spatially flat, homogeneous and isotropic, and composed of baryonic matter, dark matter, and dark energy such as the deceleration parameter, energy density, and EoS parameter. In this work, we study two cosmological models in \(f\left( Q\right) \) gravity. Now, we discuss the behavior of some of the previously mentioned cosmological parameters.

The deceleration parameter is a measure that estimates the rate of cosmic expansion. It is obtained as shown in (21). The deceleration parameter’s value can be positive, zero, or negative, and it is determined by the density of matter and the cosmological constant in the Universe. A positive value of q suggests that the expansion of the Universe is decelerating, whereas a negative value of q indicates that the expansion is accelerating. From Figs. 8 and 12, we can observe that the deceleration parameter of our models is positive (\(q>0\)) in the early Universe and negative (\(q<0\)) in the late cosmos. As a result, it shows that the cosmos is transitioning from deceleration to acceleration after a transition red-shift \(z_{t}\). Also, the Universe achieves the exponentially accelerating de Sitter stage with \(q=-1\) in the enormous time limit, a finding that is independent of the model parameters, q decreases as cosmic time increases or vice versa in terms of red-shift. This evolution is compatible with the recent Universe’s behavior, which passed through three phases: decelerating dominated, accelerating expansion, and late-time acceleration. Finally, we observe that the current values of the deceleration parameter \(q_{0}\left( z=0\right) \) and \(z_{t}\left( q=0\right) \) agree with the Hubble, SNe Ia, and BAO datasets.

Figures 9 and 13 show that the energy density of the cosmos remains positive all through the Universe’s history and decreases as cosmic time t increases in both models. It starts as positive and decreases to a small value as \(t\rightarrow \infty \) (or \(z\rightarrow -1\)). Another attempt to understand the presence of dark energy is to determine the equation of state (EoS) value and its evolution. Figures 10 and 14 show the behavior of the EoS parameter for two models. So, the EoS parameter is a dimensionless quantity that represents the pressure-to-energy density ratio in a cosmic fluid i.e. \(\omega =\frac{p}{\rho }\). It is frequently used in cosmology to explain the behavior of dark energy and dark matter, which are considered to make up the majority of the Universe. The EoS parameter can also have values varying from \(-1\) to 1. A fluid with a value of \(\omega =-1\) behaves like a cosmological constant, such as dark energy, whereas a fluid with a value of \(\omega =0\) behaves like non-relativistic matter, such as dark matter. The effective EoS parameter in Figs. 11 and 15 appears to be similar, takes negative values for all values of z. As a result, both EoS are in the quintessence region ( \(-1<\omega <-\frac{1}{3}\)), approaching the cosmological constant at high red-shifts. Finally, we observe that the current values of the EoS parameter \( \omega _{0}\left( z=0\right) \) agree with the Hubble, SNe Ia, and BAO datasets (Figs. 12, 13, 14, 15).

Fig. 8
figure 8

Evolution of the deceleration parameter vs red-shift z (Model 1)

Fig. 9
figure 9

Evolution of the energy density vs red-shift z (Model 1)

Fig. 10
figure 10

Evolution of the EoS parameter vs red-shift z (Model 1)

Fig. 11
figure 11

Evolution of the effective EoS parameter vs red-shift z (Model 1)

Fig. 12
figure 12

Evolution of the deceleration parameter vs red-shift z (Model 2)

Fig. 13
figure 13

Evolution of the energy density vs red-shift z (Model 2)

Fig. 14
figure 14

Evolution of the EoS parameter vs red-shift z (Model 2)

Fig. 15
figure 15

Evolution of the effective EoS parameter vs red-shift z (Model 2)

7 Concluding remarks

The present era of accelerating Universe has gotten increasingly intriguing over time. To develop a proper description of the accelerating Universe, a variety of dynamical DE models and modified gravity models have been used in different ways. We have investigated the accelerated expansion of the Universe in this paper by using the parametric form of the equation of state parameter in the context of f(Q) gravity, where the gravitational interaction is represented by the non-metricity scalar Q. We have studied two functional forms of f(Q), specifically \(f(Q)=-Q+\frac{\alpha }{Q}\) and \(f(Q)=-\alpha Q^{n}\), where \(\alpha \) and n are model parameters, and the parametrization form of the equation of state parameter as \(\omega \left( z\right) =-\frac{1}{1+3\beta \left( 1+z\right) ^{3}}\), where \(\beta \) is a constant parameter.

Using the aforementioned parametric form, we can derive the Hubble parameter as shown in (30) and (36). In addition, we have utilized Hubble datasets with 57 data points, SNe Ia datasets with 1048 data points, and BAO datasets with six data points to find the best-fit values for the model parameter. The values for the first model are as: \(H_{0}=69.3_{-2.2}^{+2.4}\), \(\alpha =3.9_{-4.0}^{+5.2}\), and \(\beta =0.129_{-0.024}^{+0.027}\) for Hubble datasets, \(H_{0}=67.9_{-2.6}^{+2.5}\), \(\alpha =4.0_{-4.1}^{+5.3}\), and \(\beta =0.160_{-0.074}^{+0.095}\) for SNe Ia datasets, and \( H_{0}=69_{-9}^{+10}\), \(\alpha =3.7_{-3.9}^{+5.2}\), and \(\beta =0.148_{-0.035}^{+0.037}\) for BAO datasets. In the case of the second model: \(H_{0}=65.5_{-4.6}^{+4.5}\), \(n=1.16_{-0.15}^{+0.14}\), and \(\beta =0.33_{-0.24}^{+0.39}\) for Hubble datasets. We have examined the behavior of several cosmological parameters for both models according to these best-fit model parameter values. The deceleration parameter behavior of both models has predicted the phase transition from deceleration (\(q>0\)) to acceleration (\(q<0\)). This shows that the current expansion of the Universe is accelerating. Furthermore, the effective EoS parameter in both models presently behaves in the same method as the quintessence behavior, i.e.\( -1<\omega <-\frac{1}{3}\), with the current values in the negative range and close to the observed value. It shows that the current Universe is accelerating. Finally, we have concluded that the proposed parameterization form of the EoS parameter in the context of f(Q) gravity theory plays a significant role in demonstrating the late-time accelerated expansion of the Universe.