1 Introduction

Alternative theories of gravity [1], where the Lorentz symmetry is violated, have drawn the attention of gravitation physicists in the last decades. Hořava–Lifshitz gravity and Einstein-aether gravity are two theories which have been widely studied because they provide Lorentz violation.

Hořava–Lifshitz gravity is a power-counting renormalization theory with consistent ultra-violet behavior exhibiting an anisotropic Lifshitz scaling between time and space at the ultra-violet limit, while general relativity is provided as a limit [2]. There are various physical applications of Hořava–Lifshitz gravity, some results on compact stars, black holes, universal horizons, non-relativistic gravity duality and other subjects are discussed in the review [3]. Recently, it has been found that Hořava–Lifshitz gravity is in agreement with the observations of the gravitational-wave event GW170817 [4]. Cosmological applications of Hořava–Lifshitz theory are discussed in [5,6,7,8,9] , while cosmological constraints on Hořava–Lifshitz theory using some of the recent cosmological data can be found in [10].

On the other hand, in Einstein-aether theory, the Lorentz symmetry is violated by the introduction of a unitary time-like vector field, known as the ‘æther’, in the Einstein–Hilbert action [11,12,13,14,15] . In particular, the terms quadratic in the covariant derivatives of the kinetic part of the aether lagrangian modify the gravitational action of general relativity. The Einstein-aether action is the most general second-order theory which is defined by the spacetime metric \(g_{ab}\) and the aether field \(u^{a}\) involving no more than two derivatives [14, 15]. A main property of the Einstein-aether theory is that it can be seen as the classical limit of Hořava–Lifshitz gravity. This means that exact solutions of Einstein-aether theory are solutions of Hořava–Lifshitz theory; however, the inverse is not true. While the two theories are equivalent in terms of exact solutions, in general, the equivalence is not true when the full field equations are considered [16].

Black hole solutions in Einstein-aether gravity are discussed in [17], where it was found that the exterior solution is close to the Schwarzschild solution of general relativity. On the other hand, the interior solution differs from that of general relativity. In [18], the authors studied spherically symmetric spacetimes with fluid source which describe neutron stars. A detailed study of spherically symmetric spacetimes with perfect fluid models in Einstein-æther theory was performed in Ref. [19]. The authors studied the local stability of the equilibrium points for the gravitational field equations for various perfect fluid models. The results of [19] describe inhomogeneous cosmological models and astrophysical objects. Kantowski–Sachs Einstein-aether perfect fluid models were studied in [20]. In the presence of a scalar field, static spherically symmetric solutions were determined in [21, 22]. Moreover, exact solutions of homogeneous and anisotropic spacetimes were found recently in [23, 24]. In addition, exact solutions in the presence of a modified Chaplygin gas or in the presence of a Maxwell field can be found in [25, 26] while the Einstein-aether theory has been studied as dark-energy candidate to explain the late-time acceleration phase of the universe in [27, 28].

A Lorentz violating inflationary model has been proposed by Kanno and Soda [29]. More precisely, it has been proposed a nonminimally coupling of a scalar field with the aether field, where the Einstein-aether coefficients become functions of the scalar field. In this model, the inflationary stage is divided into two parts; the Lorentz-violating stage and the standard slow-roll stage. In the first stage, the universe expands as an exact de Sitter spacetime, although the inflaton field is rolling down the potential. Another Einstein-aether scalar-field inflaton model coupled bilinearly to the expansion of the aether was proposed by Donnelly and Jacobson in [30]. In Friedmann–Lemaître–Robertson–Walker (FLRW) spacetime the inflaton-aether expansion coupling leads to a driving force on the inflaton that is proportional to the Hubble parameter. This force affects the slow-roll dynamics, but still allows a graceful exit of inflation. In [32], several families of inflationary exact solutions were derived for the model of Donnelly and Jacobson, while the effects of the Lorentz violation during the slow-roll period were studied in [33]. Recently, homogeneous and isotropic exact solutions for an Einstein-aether scalar field model were determined in [34,35,36,37], while studies of homogeneous and anisotropic models in the Einstein-aether scalar field model can be found in [38,39,40].

In this paper, we study the extension of Einstein-aether theory in Weyl geometry, specifically in Weyl integrable theory [41]. Weyl geometry is a torsion-free manifold equipped with a connection which preserves the conformal structure. In the case, where the conformal structure is analogous to a scalar field, we have the Weyl Integrable theory where the connection structure of the geometry differs from the Levi–Civita connection by a scalar field of the conformal metric, which defines the Levi–Civita connection. Hence, a scalar field can be introduced in the gravitational action integral by geometric quantities of the underlying manifold. For various applications of Weyl integrable theory in gravitational physics, we refer the reader to [42,43,44,45,46,47] and references therein. As we shall see from the following analysis, in our approach we are able to introduce a geometric scalar field coupled to the aether field, by considering the Einstein-æther action in Weyl geometry. The scalar field is introduced by the symmetric parts of the covariant derivatives for the aether field. We consider the simplest scenario of a homogeneous and isotropic spacetime, where we observe that the scalar field is introduced by the expansion rate of the aether field. The presence of an additional matter source is also discussed. The plan of this paper is as follows.

In Sect. 2, we briefly discuss the main properties of Weyl geometry and present the Einstein–Hilbert action in Weyl integrable theory. In Sect. 3, we present the model of our study, which is Einstein-aether theory in Weyl geometry. We produce the gravitational field equations and we focus on the case of a FLRW background space. We find that the field equations can be solved explicitly, where the scale factor can describe an accelerated universe for a specific range for the values of the free parameters of the model. In Sects. 4 and 5 we consider cosmological models with a dust fluid source, or with a non-zero scalar field potential respectively. Finally, in Sect. 7 we discuss our results and we draw our conclusions.

2 Weyl integrable gravity

Consider a four-dimensional manifold M described by the metric tensor, \(g_{\mu \nu }\), and the covariant derivative, \({\tilde{\nabla }}_{\mu }\), defined by the affine connection \({\tilde{\Gamma }}_{\mu \nu }^{\kappa }\) with the property,

$$\begin{aligned} {\tilde{\nabla }}_{\kappa }g_{\mu \nu }=\omega _{\kappa }g_{\mu \nu }, \end{aligned}$$
(1)

where \(\omega _{\mu }~\)is a gauge field which defines the geometry.

By definition (1), it follows that the affine connection \({\tilde{\Gamma }}_{\mu \nu }^{\kappa }\) is related to the Christoffel symbols \(\Gamma _{\mu \nu }^{\kappa }\left( g\right) \) of the metric tensor \(g_{\mu \nu }\) by the relation

$$\begin{aligned} {\tilde{\Gamma }}_{\mu \nu }^{\kappa }=\Gamma _{\mu \nu }^{\kappa }-\omega _{(\mu } \delta _{\nu )}^{\kappa }+\frac{1}{2}\omega ^{\kappa }g_{\mu \nu }, \end{aligned}$$
(2)

from which we observe the importance of the gauge vector, \(\omega _{\mu }\). When \(\omega _{\mu }\) is a gradient vector field, that is, when there exists a scalar field, \(\phi \), such that \(\omega _{\mu }=\phi _{,\mu }\), then the Weyl geometry reduces to the so-called Weyl integrable geometry.

In Weyl integrable geometry, the Ricci tensor \({\tilde{R}}_{\mu \nu }\) is related to the Riemannian Ricci tensor \(R_{\mu \nu }\) by [48],

$$\begin{aligned} {\tilde{R}}_{\mu \nu }= & {} R_{\mu \nu }-{\tilde{\nabla }}_{\nu }\left( {\tilde{\nabla }}_{\mu }\phi \right) -\frac{1}{2}\left( {\tilde{\nabla }}_{\mu }\phi \right) \left( {\tilde{\nabla }}_{\nu }\phi \right) \nonumber \\&-\frac{1}{2}g_{\mu \nu }\left( \frac{1}{\sqrt{-g}}\left( g^{\mu \nu }\sqrt{-g}\phi \right) _{,\mu \nu }\right. \nonumber \\&\left. -g^{\mu \nu }\left( {\tilde{\nabla }}_{\mu }\phi \right) \left( {\tilde{\nabla }}_{\nu } \phi \right) \right) , \end{aligned}$$
(3)

from which the Ricci scalar follows:

$$\begin{aligned} {\tilde{R}}=R-\frac{3}{\sqrt{-g}}\left( g^{\mu \nu }\sqrt{-g}\phi \right) _{,\mu \nu }+\frac{3}{2}\left( {\tilde{\nabla }}_{\mu }\phi \right) \left( {\tilde{\nabla }}_{\nu }\phi \right) . \end{aligned}$$
(4)

The simplest gravitational action Integral which can be defined in Weyl integrable theory as an extension of the Einstein–Hilbert action has been proposed to be,

$$\begin{aligned} S_{W}=\int dx^{4}\sqrt{-g}\left( {\tilde{R}}+\xi \left( {\tilde{\nabla }}_{\nu }\left( {\tilde{\nabla }}_{\mu }\phi \right) \right) g^{\mu \nu }\right) , \end{aligned}$$
(5)

where \(\xi \) is an arbitrary coupling constant,

Variation with respect to the metric tensor of (5) provides the field equations [48],

$$\begin{aligned}&{\tilde{G}}_{\mu \nu }+{\tilde{\nabla }}_{\nu }\left( {\tilde{\nabla }}_{\mu } \phi \right) -\left( 2\xi -1\right) \left( {\tilde{\nabla }}_{\mu }\phi \right) \left( {\tilde{\nabla }}_{\nu }\phi \right) \nonumber \\&\quad +\xi g_{\mu \nu }g^{\kappa \lambda }\left( {\tilde{\nabla }}_{\kappa }\phi \right) \left( {\tilde{\nabla }}_{\lambda }\phi \right) =0, \end{aligned}$$
(6)

while variation with respect to the field \(\phi \) gives

$$\begin{aligned} \left( {\tilde{\nabla }}_{\nu }\left( {\tilde{\nabla }}_{\mu }\phi \right) \right) g^{\mu \nu }+2g^{\mu \nu }\left( {\tilde{\nabla }}_{\mu }\phi \right) \left( {\tilde{\nabla }}_{\nu }\phi \right) =0 \end{aligned}$$
(7)

where \({\tilde{G}}_{\mu \nu }\) is the Einstein tensor in Weyl theory, that is, \({\tilde{G}}_{\mu \nu }={\tilde{R}}_{\mu \nu }-\frac{1}{2}{\tilde{R}}g_{\mu \nu }.\)

However, we always write the field equations by using the geometric definitions of Riemannian geometry, which means that Eq. (6) are written in a similar form

$$\begin{aligned} G_{\mu \nu }-\zeta \left( \phi _{,\mu }\phi _{,\nu }-\frac{1}{2}g_{\mu \nu } \phi ^{,\kappa }\phi _{,\kappa }\right) =0, \end{aligned}$$
(8)

where \(2\zeta \equiv 4\xi -3\), while for the scalar field \(\phi \), the second-order differential equation (7) is the usual Klein–Gordon equation \(g^{\mu \nu }\nabla _{\nu }\nabla _{\mu }\phi =0.\)

The origin of the scalar field \(\phi \) is geometrical and it is related to the nature of the affine connection, \({\tilde{\Gamma }}_{\mu \nu }^{\kappa }\).

3 Einstein-aether theory in Weyl integrable gravity

We consider the contribution of a timelike vector field, the so called aether field \(u^{\mu }\), in the gravitational action integral. In particular, we assume the gravitational action Integral

$$\begin{aligned} S_{AE}=S_{W}+S_{AE}, \end{aligned}$$
(9)

where \(S_{W}\) is the action integral (5), and \(S_{AE}\) is the action integral which is defined by the aether field defined in the Weyl manifold, that is, [12]

$$\begin{aligned} S_{AE}=\int d^{4}x\sqrt{-{\tilde{g}}}\left( {\tilde{K}}^{\alpha \beta \mu \nu } {\tilde{\nabla }}_{\alpha }u_{\mu }{\tilde{\nabla }}_{\beta }u_{\nu }+\lambda \left( {\tilde{g}}_{\mu \nu }u^{\mu }u^{\nu }+1\right) \right) , \end{aligned}$$
(10)

where the Lagrange multiplier \(\lambda ~\)ensures the unitarity of the aether field, i.e. \({\tilde{g}}_{\mu \nu }u^{\mu }u^{\nu }=-1,\) and the fourth-rank tensor \(K^{\alpha \beta \mu \nu }\) is expressed as follows:

$$\begin{aligned} {\tilde{K}}^{\alpha \beta \mu \nu }\equiv c_{1}{\tilde{g}}^{\alpha \beta }{\tilde{g}} ^{\mu \nu }+c_{2}{\tilde{g}}^{\alpha \mu }{\tilde{g}}^{\beta \nu }+c_{3}\tilde{g}^{\alpha \nu }{\tilde{g}}^{\beta \mu }+c_{4}{\tilde{g}}^{\mu \nu }u^{\alpha }u^{\beta }. \end{aligned}$$
(11)

where \({\tilde{g}}_{\mu \nu }=e^{\phi }g_{\mu \nu }\) is the conformal related metric associated with the covariant derivative \({\tilde{\nabla }}_{\mu }\). Moreover, parameters \(c_{1},~c_{2},~c_{3}\) and \(c_{4}\) are dimensionless constants and define the coupling between the æther field and gravity.

The gravitational field equations are

$$\begin{aligned}&{\tilde{G}}_{\mu \nu }+{\tilde{\nabla }}_{\nu }\left( {\tilde{\nabla }}_{\mu } \phi \right) -\left( 2\xi -1\right) \left( {\tilde{\nabla }}_{\mu }\phi \right) \left( {\tilde{\nabla }}_{\nu }\phi \right) \nonumber \\&\quad +\xi g_{\mu \nu }g^{\kappa \lambda }\left( {\tilde{\nabla }}_{\kappa }\phi \right) \left( {\tilde{\nabla }}_{\lambda }\phi \right) ={{\tilde{T}}_{\mu \nu }^{\ae },} \end{aligned}$$
(12)

where \({{\tilde{T}}_{\mu \nu }^{\ae }}\) is the energy momentum tensor which correspond to the Aether field defined in the Weyl manifold.

The energy momentum tensor \({{\tilde{T}}_{\mu \nu }^{\ae }}\) is derived by variation with respect to to the coformal metric \({\tilde{g}}_{\mu \nu }\), because the energy momentum tensor to be defined in the Weyl integrable manifold, see also the discussion in [41, 49],

$$\begin{aligned} {{\tilde{T}}_{\mu \nu }^{\ae }}&=2c_{1}({\tilde{\nabla }}_{\mu }u^{\alpha }\tilde{\nabla }_{\nu }u_{\alpha }-{\tilde{\nabla }}_{\alpha }u_{\mu }{\tilde{\nabla }}_{\beta }u_{\nu }{\tilde{g}}^{\alpha \beta })\nonumber \\&\quad +2\lambda u_{\mu }u_{\nu }+{\tilde{g}}_{\mu \nu }{\tilde{\Phi }} _{u}\nonumber \\&\quad -2[{\tilde{\nabla }}_{\alpha }(u_{(\mu }{\tilde{J}}^{\alpha }{}_{\nu )}) +{\tilde{\nabla }}_{\alpha }(u^{\alpha }{\tilde{J}}_{(a\nu )})-{\tilde{\nabla }}_{\alpha } (u_{(\mu }{\tilde{J}}_{\nu )}{}^{\alpha })]\nonumber \\&\quad -2c_{4}\left( {\tilde{\nabla }}_{\alpha }u_{\mu }u^{\alpha }\right) \left( {\tilde{\nabla }}_{\beta }u_{\nu }u^{\beta }\right) , \end{aligned}$$
(13)

where now \({{{\tilde{J}}}}{{^{a}}_{m}}\) and \(\Phi _{u}\) are defined as

$$\begin{aligned} {{{\tilde{J}}}}{{^{\mu }}_{\nu }}=-{\tilde{K}}{{^{\mu \beta }}_{\nu \alpha }{\tilde{\nabla }} }_{\beta }{u}^{\alpha }~,~{\tilde{\Phi }}_{u}=-{\tilde{K}}^{\mu \beta }{}_{\alpha \mu } {\tilde{\nabla }}_{\mu }u^{\alpha }{\tilde{\nabla }}_{\beta }u^{\mu }\,. \end{aligned}$$

Moreover, the Lagrange multiplier \(\lambda \) is defined by the constraint equation

$$\begin{aligned}&c_{4}{\tilde{g}}^{\mu \nu }u^{\alpha }{{\tilde{\nabla }}}_{\beta }u_{\nu }{\tilde{\nabla }}_{\alpha }u_{\mu }{\tilde{g}}^{\kappa \beta }-c_{4}{\tilde{g}}^{\mu \kappa }{\tilde{g}}^{\alpha \lambda }{{\tilde{\nabla }}}_{\beta }u_{\lambda }u^{\beta } {{\tilde{\nabla }}}_{\alpha }u_{\mu }\nonumber \\&\quad -c_{4}{\tilde{g}}^{\mu \kappa }u^{\alpha } {{\tilde{\nabla }}}_{\beta }u^{\beta }{{\tilde{\nabla }}}_{\alpha }u_{\mu }-\tilde{K}^{\alpha \beta \mu \kappa }{{\tilde{\nabla }}}_{\beta }{{\tilde{\nabla }}}_{\alpha }u_{\mu }-\lambda u^{\kappa }=0,\nonumber \\ \end{aligned}$$
(14)

and the field \(\phi \) satisfies Eq. (7). The energy-momentum tensor (13) has similar covariant form with the Einstein-aether theory, where now it is defined in Weyl manifold.

Using the kinematic quantities for the timelike vector field defined in they Weyl manifold, that is, the expansion rate \({\tilde{\theta }},\) the shear \({\tilde{\sigma }}_{\mu \nu }\), the vorticity \({\tilde{\omega }}_{\mu \nu }\), and the acceleration \({\tilde{\alpha }}^{\mu }\), the action integral (10) becomes [31],

$$\begin{aligned} S_{EA}=\int \sqrt{-{\tilde{g}}}dx^{4}\left( \frac{c_{\theta }}{3}\tilde{\theta }^{2}+c_{\sigma }{\tilde{\sigma }}^{2}+c_{\omega }{\tilde{\omega }}^{2}+c_{\alpha }{\tilde{\alpha }}^{2}\right) , \end{aligned}$$
(15)

where the new parameters \(c_{\theta },~c_{\sigma },~c_{\omega },~c_{a}\) are functions of \(c_{1},~c_{2},~c_{3}\) and \(c_{4}\), i.e.

$$\begin{aligned} c_{\theta }= & {} \left( c_{1}+3c_{2}+c_{3}\right) ~,\ c_{\sigma }=c_{1} +c_{3}~,\nonumber \\ c_{\omega }= & {} c_{1}-c_{3}~,\ c_{a}=c_{4}-c_{1}, \end{aligned}$$
(16)

and \({\tilde{\sigma }}^{2}={\tilde{\sigma }}_{\mu \nu }{\tilde{\sigma }}^{\mu \nu } ,~{\tilde{\omega }}^{2}={\tilde{\omega }}_{\mu \nu }{\tilde{\omega }}^{\mu \nu }\), and \({\tilde{\alpha }}^{2}={\tilde{\alpha }}_{\mu }{\tilde{\alpha }}^{\mu }\).

3.1 FLRW background spacetime

For the homogeneous and isotropic FLRW spacetime with zero spatial curvature

$$\begin{aligned} ds^{2}=-dt^{2}+a^{2}\left( t\right) \left( dx^{2}+dy^{2}+dz^{2}\right) , \end{aligned}$$
(17)

we calculate for the comoving aether field

$$\begin{aligned} {\tilde{\sigma }}^{2}=0~,~{\tilde{\omega }}^{2}=0~\text {and }{\tilde{\alpha }} ^{2}=0 \end{aligned}$$
(18)

while the energy-momentum tensor \({{\tilde{T}}_{\mu \nu }^{\ae }}\) can be written in the simplest form

$$\begin{aligned} {{\tilde{T}}_{\mu \nu }^{\ae }}={\rho }^{{\ae \ }}u_{\mu }u_{\nu }+{p}^{{\ae }}{\tilde{h}}_{\mu \nu } \end{aligned}$$
(19)

where \(\rho ^{{\ae }}=-\frac{c_{\theta }}{3}{\tilde{\theta }}^{2}\), \(~ p^{{\ae \ }}=\frac{c_{\theta }}{3}\left( 2{\tilde{\theta }}_{,t} +{\tilde{\theta }}^{2}\right) \) and \({\tilde{h}}_{\mu \nu }\) is the projection tensor \(h_{\mu \nu }={\tilde{g}}_{\mu \nu }+u_{\mu }u_{\nu }\).

At this point we remark that \({\tilde{\theta }}\) is the expansion rate defined in the Weyl manifold, and it is expanded as \({\tilde{\theta }}=\theta -{\dot{\phi }}\), where \(\theta \) is the expansion rate of General Relativity. From the latter it is clear that an interaction between the scalar field and the aether field exists. That interaction has geometric origin and it is direct related with the affine connection \({\tilde{\Gamma }}_{\mu \nu }^{\kappa }\).

For the line element (17), the gravitational field equations are

$$\begin{aligned}&\frac{\theta ^{2}}{3}-\rho _{\phi }-\rho ^{{\ae \ }}=0 \end{aligned}$$
(20)
$$\begin{aligned}&{\dot{\theta }}+\frac{\theta ^{2}}{3}+\frac{1}{2}\left( \rho _{\phi }+3p_{\phi }\right) +\frac{1}{2}\left( \rho ^{{\ae \ }}+3p^{{\ae }}\right) =0 \end{aligned}$$
(21)

where \(\theta \) is the expansion rate of general relativity, and \(\rho _{\phi }=\frac{\zeta }{2}{\dot{\phi }}^{2}\)\(p_{\phi }=\frac{\zeta }{2} {\dot{\phi }}^{2}\) are the energy density and pressure for the field \(\phi \).

The field equations (20)–(21) can be written in an equivalent way, as follows:

$$\begin{aligned}&\left( 1+c_{\theta }\right) \frac{\theta ^{2}}{3}-\frac{2}{3}c_{\theta } \theta {\dot{\phi }}-\left( \frac{\zeta }{2}-\frac{c_{\theta }}{3}\right) {\dot{\phi }}^{2}=0, \end{aligned}$$
(22)
$$\begin{aligned}&\frac{\left( 1+c_{\theta }\right) }{3}{\dot{\theta }}+\left( 1+c_{\theta }\right) \theta ^{2}-\frac{2}{3}c_{\theta }\theta {\dot{\phi }}\nonumber \\&\quad +\left( \frac{c_{\theta }}{3}+\zeta \right) {\dot{\phi }}^{2}-c_{\theta }\ddot{\phi }=0, \end{aligned}$$
(23)

while for the field \(\phi \), we have

$$\begin{aligned} \left( 2c_{\theta }-3\left( 1+c_{\theta }\right) \zeta \right) \ddot{\phi }+3\zeta c_{\theta }{\dot{\phi }}^{2}-3\left( 1+c_{\theta }\right) \zeta \theta {\dot{\phi }}=0. \end{aligned}$$
(24)

From (23) and (24), we find

$$\begin{aligned} \frac{{\dot{\theta }}}{\theta ^{2}}= & {} -\frac{1}{3}-\frac{c_{\theta }}{3\left( 1+c_{\theta }\right) \left( 3\zeta \left( 1+c_{\theta }\right) -2c_{\theta }\right) }\left( \frac{{\dot{\phi }}}{\theta }\right) \nonumber \\&\quad +\frac{\left( 2c_{\theta }-3\zeta \right) \left( c_{\theta }+3\zeta \left( 1+c_{\theta }\right) \right) }{3\left( 1+c_{\theta }\right) \left( 3\zeta \left( 1+c_{\theta }\right) -2c_{\theta }\right) }\left( \frac{{\dot{\phi }}}{\theta }\right) ^{2}, \end{aligned}$$
(25)

from which we infer that the parameter for the equation of state for the effective fluid is \(w_{eff}=-1-2\frac{{\dot{\theta }}}{\theta ^{2}}.\)

We have introduced a scalar field in Einstein-aether theory by using the Weyl geometry. Other attempts to introduce a scalar field in Einstein-aether theory have been proposed before by Kanno and Soda in [29] and latter by Donnelly and Jacobson in [30]. In both attempts, the scalar field interacts with the aether field. In [29], the interaction has been proposed to be in the coupling parameters \(c_{I}\) of the aether field. A more general consideration has been proposed in [30], where an arbitrary potential for the scalar field has been introduced in the action integral such that kinematic quantities of the aether field are included. Our approach is totally geometric and the interaction between the scalar field and the aether field is introduced by the covariant derivative which defines the kinematic quantities and the affine connection \({\tilde{\Gamma }}_{\mu \nu }^{\kappa }\) of the Weyl geometry.

3.1.1 Exact solution

We select the new dimensionless variable \(x\left( t\right) =\frac{2c_{\phi } }{1+c_{\phi }}\frac{{\dot{\phi }}}{\theta }\), to write the field equations as,

$$\begin{aligned}&-12c_{\theta }^{2}\left( 3\zeta +c_{\theta }\left( 3\zeta -2\right) \right) \frac{dx}{d\ln a} \nonumber \\&\quad =x\left( 8c_{\theta }^{2}\left( c_{\theta } +3\zeta \left( 1+c_{\theta }\right) \right) -8c_{\theta }^{2}\left( c_{\theta }+3\zeta \left( 1+c\theta \right) \right) x\right) \nonumber \\&\qquad +x^{3}\left( 1+c_{\theta }\right) \left( 3c_{\theta }\left( 1-3\zeta \right) \zeta -9\zeta ^{2}+\left( 2+6\zeta \right) c_{\theta } ^{2}\right) , \end{aligned}$$
(26)

while the constraint equation becomes,

$$\begin{aligned} 1+\left( \frac{\left( 1+c_{\theta }\right) \left( 2c_{\theta } -3\zeta \right) }{8c_{\theta }}x-1\right) x=0. \end{aligned}$$
(27)

From the latter algebraic equation we get,

$$\begin{aligned} \mathbf {x^{\pm }}=\frac{2c_{\theta }}{1+c_{\theta }}\frac{2c_{\theta }\pm \sqrt{6\zeta +c_{\theta }\left( 6\zeta -4\right) }}{2c_{\theta }-3\zeta },~2c_{\theta }-3\zeta \ne 0. \end{aligned}$$
(28)

The points are real when \(6\zeta +c_{\theta }\left( 6\zeta -4\right) \ge 0\,.\)These two points provide ideal-gas exact solutions with an equation of state parameter,

$$\begin{aligned}&w_{eff}\left( \mathbf {x^{\pm }}\left( c_{\theta },\zeta \right) \right) \nonumber \\&\quad =-1-2\frac{{\dot{\theta }}}{\theta ^{2}}\nonumber \\&\quad =\frac{9\zeta ^{2}+c_{\theta }^{2}\left( 6\zeta -4\right) +3\zeta c_{\theta }\left( 3\zeta +2\sqrt{6\left( 1+c_{\theta }\right) \zeta -4c_{\theta }}\right) }{\left( 2c_{\theta }-3\zeta \right) \left( 3\left( 1+c_{\theta }\right) \zeta -2c_{\theta }\right) }.\nonumber \\ \end{aligned}$$
(29)

In the special limit where \(\zeta =0\), the exact solution is that of a stiff fluid, that is, \(w_{eff}\left( \mathbf {x^{\pm }}\left( c_{\theta },0\right) \right) =1\).

In Fig. 1 we present the region in the space of parameters \(\left\{ c_{\theta },\zeta \right\} \) in which the exact solutions at the critical points\(~x^{\pm }\) describe accelerated universes, that is, \(w_{eff}\left( \mathbf {x^{\pm } }\left( c_{\theta },\zeta \right) \right) <-\frac{1}{3}\).

Fig. 1
figure 1

Region plots in the space of parameters \(\left\{ c_{\theta },\zeta \right\} \). The exact solutions at the critical points \(\mathbf {x^{+}}\) (left fig.) and \(\mathbf {x^{-}}\) (right fig.) describe accelerated universes

4 In the presence of matter

If we consider now the existence of a pressureless fluid source, then the field equations become,

$$\begin{aligned} G_{\mu \nu }-\zeta \left( \phi _{,\mu }\phi _{,\nu }-\frac{1}{2}g_{\mu \nu } \phi ^{,\kappa }\phi _{,\kappa }\right) ={T_{ab}^{\ae }+}T_{\mu \nu }^{\left( m\right) }, \end{aligned}$$
(30)

where \(T_{\mu \nu }^{\left( m\right) }=e^{-\frac{\phi }{2}}\rho _{m}u_{\mu }u_{\nu }\), while the conservation equation for the matter field reads \({\tilde{\nabla }}_{\nu }T^{\left( m\right) \mu \nu }=0\).

For the spatially-flat FLRW line element (17), the conservation equation for the dust fluid reads,

$$\begin{aligned} {\dot{\rho }}_{m}+\left( \theta -{\dot{\phi }}\right) \rho _{m}=0, \end{aligned}$$
(31)

while the rest of the field equations become

$$\begin{aligned}&\left( 1+c_{\theta }\right) \frac{\theta ^{2}}{3}-\frac{2}{3}c_{\theta } \theta {\dot{\phi }}-\left( \frac{\zeta }{2}-\frac{c_{\theta }}{3}\right) {\dot{\phi }}^{2}\nonumber \\&\quad -e^{-\frac{\phi }{2}}\rho _{m}=0, \end{aligned}$$
(32)
$$\begin{aligned}&\frac{\left( 1+c_{\theta }\right) }{3}{\dot{\theta }}+\left( 1+c_{\theta }\right) \theta ^{2}-\frac{2}{3}c_{\theta }\theta {\dot{\phi }}\nonumber \\&\quad +\left( \frac{c_{\theta }}{3}+\zeta \right) {\dot{\phi }}^{2}-c_{\theta }\ddot{\phi } +\frac{1}{2}e^{-\frac{\phi }{2}}\rho _{m}=0, \end{aligned}$$
(33)

while the Klein–Gordon equation is modified as

$$\begin{aligned}&\left( 2c_{\theta }-3\left( 1+c_{\theta }\right) \zeta \right) \ddot{\phi }+3\zeta c_{\theta }{\dot{\phi }}^{2}-3\left( 1+c_{\theta }\right) \zeta \theta {\dot{\phi }}\nonumber \\&\quad +3\left( 1-c_{\theta }\right) \rho _{m}e^{-\frac{\phi }{2}}=0. \end{aligned}$$
(34)

With the use of the dimensionless variables,

$$\begin{aligned} x\equiv \frac{2c_{\phi }}{1+c_{\phi }}\frac{{\dot{\phi }}}{\theta }\text { and } \Omega _{m}\equiv \frac{3e^{-\frac{\phi }{2}}\rho _{m}}{\left( 1+c_{\theta }\right) \theta ^{2}}, \end{aligned}$$
(35)

the field equations reduce to the following algebraic-differential system,

$$\begin{aligned}&\Omega _{m}=1+\left( \frac{\left( 1+c_{\theta }\right) \left( 2c_{\theta }-3\zeta \right) }{8c_{\theta }}x-1\right) x, \end{aligned}$$
(36)
$$\begin{aligned}&-16c_{\theta }^{2}\left( 3\left( 1+c_{\theta }\right) \zeta -2c_{\theta }\right) \frac{dx}{d\ln a}\nonumber \\&\quad =\left( \left( 1+c_{\theta }\right) \left( 2c_{\theta }-3\zeta \right) x^{2}+8c_{\theta }^{2}\left( 1-x\right) \right) \nonumber \\&\qquad \times \left( \left( 1+c_{\theta }\right) \left( c_{\theta }+3\zeta \right) x+2\left( 1-c_{\theta }\right) c_{\theta }\right) . \end{aligned}$$
(37)

The stationary points of Eq. (37) are,

$$\begin{aligned} \mathbf {x^{\pm }}= & {} \frac{2c_{\theta }}{1+c_{\theta }}\frac{2c_{\theta }\pm \sqrt{6\zeta +c_{\theta }\left( 6\zeta -4\right) }}{2c_{\theta }-3\zeta }\text { and }\\ \mathbf {x^{0}}= & {} \frac{2c_{\theta }\left( c_{\theta }-1\right) }{\left( 1+c_{\theta }\right) \left( c_{\theta }+3\zeta \right) }. \end{aligned}$$

Points \(\mathbf {x^{\pm }}\) are the critical points of the vacuum case, \(\Omega _{m}\left( \mathbf {x^{\pm }}\right) =0\), while point \(\mathbf {x^{0}}\) provides

$$\begin{aligned} \Omega _{m}\left( \mathbf {x^{0}}\right) =\frac{c_{\theta }+9\zeta ^{2}}{\left( c_{\theta }+3\zeta \right) ^{2}}-\frac{3\left( 1+c_{\theta }\left( c_{\theta }-10\right) \right) \zeta }{\left( 1+c_{\theta }\right) \left( c_{\theta }+3\zeta \right) ^{2}}. \end{aligned}$$
(38)

Point \(\mathbf {x^{0}}\) is physically accepted when \(0\le \Omega _{m}\left( \mathbf {x^{0}}\right) \le 1\), as presented in Fig. 2.

As far as the equation of state for the cosmological solution at point \(\mathbf {x^{0}}\) is concerned, it is calculated to be,

$$\begin{aligned} w_{eff}\left( \mathbf {x^{0}}\left( c_{\theta },\zeta \right) \right) =\frac{1-c_{\theta }}{2\left( c_{\theta }+6\zeta \right) }. \end{aligned}$$
(39)

The point \(\mathbf {x^{0}}~\)describes an accelerated universe, i.e. \(w\left( \mathbf {x^{0}}\left( c_{\theta },\zeta \right) \right) <-\frac{1}{3}\), when \(\zeta ,c_{\theta }\) take values in the range of Fig. 2.

Fig. 2
figure 2

Region plots in the space of parameters \(\left\{ c_{\theta },\zeta \right\} \), where the critical point \(\mathbf {x^{0}}\) is physically accepted (left fig.) and the exact solution at this point describes an accelerated universe (right fig.)

As far as the stability of the stationary points is concerned, we determine the eigenvalue \(e\left( {\mathbf {x}}\left( c_{\theta },\zeta \right) \right) \) of Eq. (37). A stationary point is an attractor when  \({\text {Re}}\left( e\left( {\mathbf {x}}\left( c_{\theta },\zeta \right) \right) \right) <0\,\). In particular, we find that \(e\left( \mathbf {x^{0} }\left( c_{\theta },\zeta \right) \right) =\frac{\left( 1+c_{\theta }\right) \left( c_{\theta }+3\zeta \right) }{6\left( 1+c_{\theta }\right) \zeta -4c_{\theta }}\Omega _{m}\left( \mathbf {x^{0}}\right) \), while, for the other two points, it follows

$$\begin{aligned} e\left( \mathbf {x^{\pm }}\left( c_{\theta },\zeta \right) \right)= & {} \frac{18\zeta ^{2}+6c_{\theta }\left( 3\zeta -1\right) +c_{\theta }^{2}\left( 6\zeta -4\right) \mp 3\left( \zeta -c_{\theta }\left( 2+9\zeta \right) \right) \sqrt{6\left( 1+c_{\theta }\right) \zeta -4c_{\theta }}}{2\left( 2c_{\theta }-3\zeta \right) \left( 3\left( 1+c_{\theta }\right) \zeta -2c_{\theta }\right) }. \end{aligned}$$
(40)

In Fig. 3, we present the regions in the space of parameters, \(\left( c_{\theta },\zeta \right) \), where the stationary points are stable.

Fig. 3
figure 3

Region plots in the space of parameters \(\left\{ c_{\theta },\zeta \right\} \), where the stationary points: \(\mathbf {x^{0}}\) (left fig.), \(\mathbf {x^{+}}\) (middle Fig.) and \(\mathbf {x^{-}}\) (right fig.) are attractors for the cosmological model with a dust fluid source

In the limit where \(\zeta =0\), for point \(\mathbf {x^{0}}\) it follows \(\Omega _{m}\left( \mathbf {x^{0}}\right) =\frac{1}{c_{\theta }}\) and \(w_{eff}\left( \mathbf {x^{0}}\left( c_{\theta },2\right) \right) =-\frac{1}{2}\left( 1-c_{\theta }\right) ,\) hence it follows that point \(\mathbf {x^{0}}\) is physically accepted when \(c_{\theta }\ge 1\) while for \(c_{\theta }>3\) the exact solution describes an accelerated universe. That is a very interesting result because \(\Omega _{m}\left( \mathbf {x^{0}}\right) \) becomes zero only for very large values of the coupling constant \(c_{\theta }\), while \(e\left( \mathbf {x^{0}}\left( c_{\theta },0\right) \right) \) is found to be always negative for \(c_{\theta }>1\) which means that \(\mathbf {x^{0}}\) is a future attractor.

5 Evolution of a scalar field under potential \(V\left( \phi \right) \) in vacuum

In the presence of a potential term \(V\left( \phi \right) \), the field equations (20)–(21) become,

$$\begin{aligned}&\left( 1+c_{\theta }\right) \frac{\theta ^{2}}{3}-\frac{2}{3}c_{\theta } \theta {\dot{\phi }}-\left( \frac{\zeta }{2}-\frac{c_{\theta }}{3}\right) {\dot{\phi }}^{2}-V\left( \phi \right) =0, \nonumber \\ \end{aligned}$$
(41)
$$\begin{aligned}&\frac{\left( 1+c_{\theta }\right) }{3}{\dot{\theta }}+\left( 1+c_{\theta }\right) \theta ^{2}-\frac{2}{3}c_{\theta }\theta {\dot{\phi }}\nonumber \\&\quad +\left( \frac{c_{\theta }}{3}+\zeta \right) {\dot{\phi }}^{2}-c_{\theta }\ddot{\phi }-V\left( \phi \right) =0, \end{aligned}$$
(42)

while the field \(\phi \) obeys,

$$\begin{aligned}&\left( 2c_{\theta }-3\left( 1+c_{\theta }\right) \zeta \right) \ddot{\phi }+3\zeta c_{\theta }{\dot{\phi }}^{2}-3\left( 1+c_{\theta }\right) \zeta \theta {\dot{\phi }}\nonumber \\&\quad -3\left( 1+c_{\theta }\right) V_{,\phi }=0. \end{aligned}$$
(43)

We define the new dimensionless variables, x,  from (35) and \(y^{2}\equiv \frac{3V\left( \phi \right) }{\left( 1+c_{\theta }\right) \theta ^{2}}\), and the field equations become,

$$\begin{aligned}&1+\left( \frac{\left( 1+c_{\theta }\right) \left( 2c_{\theta } -3\zeta \right) }{8c_{\theta }}x-1\right) x-y^{2}=0, \end{aligned}$$
(44)
$$\begin{aligned}&-24c_{\theta }^{2}\left( 3\left( 1+c_{\theta }\right) \zeta -2c_{\theta }\right) \frac{dx}{d\ln a} \nonumber \\&\quad =\left( \left( 1+c_{\theta }\right) \left( 2c_{\theta }-3\zeta \right) x^{2}+8c_{\theta }^{2}\left( 1-x\right) \right) \nonumber \\&\qquad \times \left( x\left( 2c_{\theta }+3\left( 1+c_{\theta }\right) \left( 2\zeta +\lambda c_{\theta }\right) \right) \right. \nonumber \\&\qquad \left. -6\lambda c_{\theta }\left( 1+c_{\theta }\right) \right) , \end{aligned}$$
(45)

where \(\lambda =-\frac{V_{\phi }}{V}\). For the scalar field, we assume the exponential potential \(V\left( \phi \right) =V_{0}e^{-\lambda \phi }\) where \(\lambda =\)constant. The exponential potential is of special interest, mathematically and physically. In terms of mathematics, it reduces the dimension of the dynamical system, while, in terms of physics ,when \(\lambda =2\), it includes the case of the cosmological constant term in Weyl integrable theory [48].

The latter dynamical system admits the stationary points \(\mathbf {x^{\pm }}\) and \(\mathbf {x^{1}}=\frac{2c_{\theta }\lambda }{3\zeta +c_{\theta }\lambda }\). The new stationary point \(\mathbf {x^{1}}\) describes a universe where the scalar field potential also contributes to the effective cosmological fluid, while the equation of state parameter \(w_{eff}\left( \mathbf {x^{1}}\right) \) for the effective fluid at the stationary point is found to be,

$$\begin{aligned} w_{eff}\left( \mathbf {x^{1}}\left( c_{\theta },\zeta ,\lambda \right) \right) =\frac{\lambda \left( c_{\theta }\left( \lambda -1\right) +\lambda \right) -3\zeta }{3\zeta +c_{\theta }\lambda }. \end{aligned}$$
(46)

The corresponding eigenvalues are derived

$$\begin{aligned}&2(3({c_{\theta }}+1)\zeta -2{c_{\theta }})({c_{\theta } }\lambda +3\zeta )^{2}e\left( \mathbf {x^{1}}\left( c_{\theta },\zeta ,\lambda \right) \right) \nonumber \\&\quad =-12\zeta ^{2}(3({c_{\theta }} +1)\zeta +{c_{\theta }})\nonumber \\&\qquad +{c_{\theta }}({c_{\theta } }+1)\lambda ^{3}(3({c_{\theta }}+1)\zeta -2{c_{\theta } })\nonumber \\&\qquad +2{c_{\theta }}\lambda ^{2}(3({c_{\theta }}+1)\zeta -2{c_{\theta }})\nonumber \\&\qquad -2{c_{\theta }}\zeta \lambda (15({c_{\theta } }+1)\zeta -4{c_{\theta }}) \end{aligned}$$
(47)
$$\begin{aligned}&E_{1}e\left( \mathbf {x^{+}}\left( c_{\theta },\zeta ,\lambda \right) \right) \nonumber \\&\quad =\left( 8{c_{\theta }}^{2}-4{c_{\theta }}\left( \sqrt{6({c_{\theta }}+1)\zeta -4{c_{\theta }}}+2\right) \right) \nonumber \\&\qquad \times \left( \begin{array}[c]{c} 3{c_{\theta }}({c_{\theta }}+1)^{2}\lambda (2{c_{\theta } }-3\zeta )-{c_{\theta }}\left( \sqrt{6({c_{\theta }} +1)\zeta -4{c_{\theta }}}+2\right) \times \\ (6({c_{\theta }}+1)\zeta +3({c_{\theta }}+1){c_{\theta } }\lambda +2{c_{\theta }}) \end{array} \right) \nonumber \\&\qquad -8({c_{\theta }}-1){c_{\theta }}^{2}\left( \sqrt{6({c_{\theta }}+1)\zeta -4{c_{\theta }}}-{c_{\theta } }+1\right) \nonumber \\&\qquad \times (6({c_{\theta }}+1)\zeta +3({c_{\theta }} +1){c_{\theta }}\lambda +2{c_{\theta }}). \end{aligned}$$
(48)
$$\begin{aligned}&E_{1}e\left( \mathbf {x^{-}}\left( c_{\theta },\zeta ,\lambda \right) \right) \nonumber \\&\quad =8({c_{\theta }}-1){c_{\theta }}^{2}\left( \sqrt{6({c_{\theta }}+1)\zeta -4{c_{\theta }}}+{c_{\theta } }-1\right) \nonumber \\&\qquad \times (6({c_{\theta }}+1)\zeta +3({c_{\theta }} +1){c_{\theta }}\lambda +2{c_{\theta }})\nonumber \\&\qquad +4{c_{\theta } }\left( \sqrt{6({c_{\theta }}+1)\zeta -4{c_{\theta }} }+2{c_{\theta }}-2\right) \nonumber \\&\qquad +\left( \begin{array}[c]{l} 3{c_{\theta }}({c_{\theta }}+1)^{2}\lambda (2{c_{\theta } }-3\zeta )+{c_{\theta }}\left( \sqrt{6({c_{\theta }} +1)\zeta -4{c_{\theta }}}-2\right) \\ \times (6({c_{\theta }}+1)\zeta +3({c_{\theta }}+1){c_{\theta }}\lambda +2{c_{\theta }}) \end{array} \right) . \end{aligned}$$
(49)

where \(E_{1}=12c_{\theta }^{2}(c_{\theta }+1)(2c_{\theta }-3\zeta )(2c_{\theta }-3(c_{\theta }+1)\zeta )\).

We focus on the special value where \(\lambda =2~\) which, as we described before, corresponds to the case of the cosmological constant. In Fig. 4, we present the regions in the space of the free parameters \(\left( c_{\theta },\zeta \right) ,\) where the three stationary points \(\mathbf {x^{1}},~\mathbf {x^{\pm }}\) are attractors, while in Fig. 5 we present the region of the free parameters \(\left( c_{\theta },\zeta \right) \) where \(w_{eff}\left( \mathbf {x^{1}}\left( c_{\theta },\zeta ,\lambda \right) \right) <-\frac{1}{3}\) for \(\lambda =2\).

Fig. 4
figure 4

Region plots in the space of parameters \(\left\{ c_{\theta },\zeta \right\} \), where the stationary points: \(\mathbf {x^{1}}\) (left fig.), \(\mathbf {x^{+}}\) (middle fig.) and \(\mathbf {x^{-}}\) (right fig.) are attractors for the dynamical system with an exponential potential and \(\lambda =2\)

Recall that for point \(\mathbf {x^{1}}\) is physically acceptable as long as \(y^{2}\left( \mathbf {x^{1}}\right) \ge 0,\) where for \(\lambda =2\) we find the range \(\left\{ c_{\theta }\le -1,~\zeta \ne -\frac{2}{3}c_{\theta }\right\} \)\(\left\{ -1<c_{\theta }<0,~\zeta <\frac{2}{3}c_{\theta }~,~\zeta >\frac{2}{3}\right\} \), \(\left\{ 0<c_{\theta }\le 1~,~\zeta \ne -\frac{2}{3}c_{\theta }\right\} \)\(\left\{ c_{\theta }>1~,~\zeta \ne -\frac{2}{3}c_{\theta } ,~\zeta <\frac{2}{3}~,~\zeta >\frac{2}{3}c_{\theta }~\right\} \).

In the case where \(\zeta =0\), point \(\mathbf {x^{1}}\) is physically accepted when \(c_{\theta }\ge 1\), while \(w_{eff}\left( \mathbf {x^{1}}\left( c_{\theta },\zeta ,\lambda \right) \right) =-1+\lambda \left( 1+\frac{1}{c_{\theta } }\right) \), and the point is an attractor when \(\lambda \le -1\) or \(\left\{ -1<\lambda <0,~c_{\theta }>-\left( 1+\frac{2}{\lambda }\right) \right\} \).

6 Evolution of a scalar field with potential \(V\left( \phi \right) \) and matter

In this example the Friedmann equation, Raychaudhuri, and Klein–Gordon equations are, respectively,

$$\begin{aligned}&(2c_{\theta }-3\zeta ){\dot{\phi }}^{2}-4c_{\theta }\theta {\dot{\phi }}+2(c_{\theta }+1)\theta ^{2}\nonumber \\&\quad -6\rho _{m}e^{-\frac{\phi }{2}}-6V(\phi )=0, \end{aligned}$$
(50)
$$\begin{aligned}&6c_{\theta }{\dot{\theta }}-4c_{\theta }\theta {\dot{\phi }}+2(c_{\theta } +1)\theta ^{2}-6c_{\theta }\ddot{\phi }\nonumber \\&\quad +2c_{\theta }{\dot{\phi }}^{2}+6\zeta \dot{\phi }^{2}+6{\dot{\theta }}+3\rho _{m}e^{-\frac{\phi }{2}}-6V(\phi )=0, \end{aligned}$$
(51)
$$\begin{aligned}&3(c_{\theta }-1)\rho _{m}-2e^{\frac{\phi }{2}}\left( 3\zeta {\dot{\phi }}\left( c_{\theta }\theta -c_{\theta }{\dot{\phi }}+\theta \right) \right. \nonumber \\&\quad \left. +(3(c_{\theta } +1)\zeta -2c_{\theta })\ddot{\phi }+3(c_{\theta }+1)V^{\prime }(\phi )\right) =0, \end{aligned}$$
(52)

plus the conservation equation (31).

Fig. 5
figure 5

Region plots in the space of parameters \(\left\{ c_{\theta },\zeta \right\} \), where the stationary point \(\mathbf {x^{1}}\) describes and accelerated universe when \(\lambda =2\), that is \(w_{eff}\left( \mathbf {x^{1}}\left( c_{\theta } ,\zeta ,2\right) \right) <-\frac{1}{3}\)

Defining \(x\equiv \frac{2c_{\phi }}{1+c_{\phi }}\frac{{\dot{\phi }}}{\theta }\) and \(\Omega _{m}\equiv \frac{3e^{-\frac{\phi }{2}}\rho _{m}}{\left( 1+c_{\theta }\right) \theta ^{2}},\) \(y^{2}\equiv \frac{3V(\phi )}{(c_{\theta }+1)\theta ^{2}} \), we obtain the dynamical system:

$$\begin{aligned} x^{\prime }&=-\frac{(c_{\theta }+1)^{2}(2c_{\theta }-3\zeta )(c_{\theta } +3\zeta )x^{3}}{16c_{\theta }^{2}(3(c_{\theta }+1)\zeta -2c_{\theta })}\nonumber \\&\quad +x\left( \frac{-3c_{\theta }^{2}-3c_{\theta }\zeta +c_{\theta }-3\zeta }{6(c_{\theta }+1)\zeta -4c_{\theta }}\right. \nonumber \\&\quad \left. +\frac{(c_{\theta }+1)y^{2}(c_{\theta }(2\lambda -1)+3\zeta )}{4c_{\theta }-6(c_{\theta }+1)\zeta }\right) \nonumber \\&\quad +\frac{(c_{\theta }-1)c_{\theta }}{3(c_{\theta }+1)\zeta -2c_{\theta }}\nonumber \\&\quad +\frac{(c_{\theta }+1)(c_{\theta }(6c_{\theta }+9\zeta -2)+3\zeta )x^{2} }{8c_{\theta }(3(c_{\theta }+1)\zeta -2c_{\theta })}\nonumber \\&\quad +\frac{c_{\theta }(2(c_{\theta }+1)\lambda -c_{\theta }+1)y^{2}}{3(c_{\theta }+1)\zeta -2c_{\theta }} , \end{aligned}$$
(53)
$$\begin{aligned} 2y^{\prime }&=-\frac{(c_{\theta }+1)^{2}(2c_{\theta }-3\zeta )(c_{\theta } +3\zeta )x^{2}y}{8c_{\theta }^{2}(3(c_{\theta }+1)\zeta -2c_{\theta })}\nonumber \\&\quad +xy\left( \frac{c_{\theta }(c_{\theta }+1)}{3(c_{\theta }+1)\zeta -2c_{\theta }} -\frac{(c_{\theta }+1)\lambda }{2c_{\theta }}\right) \nonumber \\&\quad +\frac{(c_{\theta }+1)y^{3}(-2c_{\theta }\lambda +c_{\theta }-3\zeta )}{3(c_{\theta }+1)\zeta -2c_{\theta }}\nonumber \\&\quad -\frac{(c_{\theta }+1)(c_{\theta }-3\zeta )y}{3(c_{\theta }+1)\zeta -2c_{\theta }}, \end{aligned}$$
(54)

where \(\lambda =-\frac{V_{\phi }}{V}\) and for the scalar field potential we assume the exponential potential \(V\left( \phi \right) =V_{0}e^{-\lambda \phi }\) with constant \(\lambda \). The fractional energy density of matter is

$$\begin{aligned} \Omega _{m}=\frac{(c_{\theta }+1)(2c_{\theta }-3\zeta )x^{2}}{8c_{\theta }^{2} }-x-y^{2}+1. \end{aligned}$$
(55)
Fig. 6
figure 6

Region plots in the space of parameters \(\left\{ c_{\theta },\zeta \right\} \), where the stationary point \({\mathbf {x}}^{2}\) is physically acceptable, for \(\lambda =2\) (left fig.) and \(\lambda =0\) (right fig.)

Fig. 7
figure 7

Region plots in the space of parameters \(\left\{ c_{\theta },\zeta \right\} \), where the stationary points \(\mathbf {x^{0}}\)\(\mathbf {x^{\pm }}\)\(\mathbf {x^{1}}\) and \(\mathbf {x^{2}}\) are attractors, for \(\lambda =2\)

Fig. 8
figure 8

Region plots in the space of parameters \(\left\{ c_{\theta },\zeta \right\} \), where the stationary points \(\mathbf {x^{\pm }}\)\(\mathbf {x^{1}}\) and \(\mathbf {x^{2} }\) are attractors, for \(\lambda =0.\) Point \({\mathbf {x}}^{{\mathbf {0}}}\) always describes an unstable solution

Therefore, the phase plane is defined by

$$\begin{aligned} \left\{ (x,y):0\le x-\frac{(c_{\theta }+1)(2c_{\theta }-3\zeta )x^{2} }{8c_{\theta }^{2}}+y^{2}\le 1,y\ge 0\right\} . \end{aligned}$$
(56)

The stationary points of the dynamical system (53), (54) are the points \(\mathbf {x^{0}}\)\(\mathbf {x^{\pm }}\)\(\mathbf {x^{1}}\) and a new point \({\mathbf {x}}^{2}\), with coordinates

$$\begin{aligned} {\mathbf {x}}^{2}=\left( \frac{4c_{\theta }}{\left( 1+2\lambda \right) \left( c_{\theta }+1\right) },\frac{\sqrt{-2c_{\theta }\lambda +c_{\theta } +6\zeta +2\lambda +1}}{\sqrt{c_{\theta }+1}(2\lambda +1)}\right) . \end{aligned}$$

For the exact solution at the point \({\mathbf {x}}^{2}\) we find \(\Omega _{m}\left( {\mathbf {x}}^{2}\right) =\frac{2\left( \lambda \left( 1-c_{\theta }+2\left( 1+c_{\theta }\right) \lambda \right) -6\zeta \right) }{\left( 1+c_{\theta }\right) \left( 1+2\lambda \right) ^{2}}\) and

$$\begin{aligned}&3({c_{\theta }}+1)(2\lambda +1)^{2}(3({c_{\theta }}+1)\zeta -2{c_{\theta }})w_{eff}\left( \mathbf {x^{2}}\left( c_{\theta } ,\zeta ,\lambda \right) \right) \nonumber \\&\quad =36({c_{\theta }}+1)\zeta ^{2}+4({c_{\theta }}+1)\lambda ^{2}({c_{\theta }}(3{c_{\theta }}-1)-3({c_{\theta }}+1)\zeta )\nonumber \\&\qquad -2\lambda (6({c_{\theta }}+1)({c_{\theta }}+2)\zeta +{c_{\theta }}({c_{\theta }}(3{c_{\theta }} -10)-5))+\nonumber \\&\qquad -3({c_{\theta }}(15{c_{\theta }}+2)+3)\zeta +6{c_{\theta } }({c_{\theta }}1), \end{aligned}$$
(57)

which means that the stationary point \({\mathbf {x}}^{2}\) describes a scaling solution. The point is physically acceptable when \(0\le \Omega _{m}\left( {\mathbf {x}}^{2}\right) \le 1\). For the specific cases, \(\lambda =0\) and \(\lambda =2,\) in Fig. 6 we present the range of the parameters \(\left\{ c_{\theta },\zeta \right\} \) for which the point is physically acceptable. We continue with the stability analysis of the stationary points.

For point \(\mathbf {x^{0}}\), the eigenvalues of the linearized system are derived

$$\begin{aligned}&e_{1}\left( \mathbf {x^{0}}\left( c_{\theta },\zeta ,\lambda \right) \right) \\&\quad =\frac{-18(c_{\theta }+1)\zeta ^{2}+3((c_{\theta }-10)c_{\theta } +1)\zeta -2c_{\theta }(c_{\theta }+1)}{4(c_{\theta }+3\zeta )(3(c_{\theta } +1)\zeta -2c_{\theta })},\\&e_{2}\left( \mathbf {x^{0}}\left( c_{\theta },\zeta ,\lambda \right) \right) =\frac{-2c_{\theta }\lambda +c_{\theta }+6\zeta +2\lambda +1}{4c_{\theta }+12\zeta }. \end{aligned}$$

Similarly, for the rest of the points we find

$$\begin{aligned}&e_{1}\left( \mathbf {x^{+}}\left( c_{\theta },\zeta ,\lambda \right) \right) \\&\quad =\sqrt{6(c_{\theta }+1)\zeta -4c_{\theta }}\left( \frac{3}{3\zeta -2c_{\theta }} -\frac{4}{3(c_{\theta }+1)\zeta -2c_{\theta }}\right) \\&\qquad +\frac{3c_{\theta } }{3\zeta -2c_{\theta }}+1,\\&e_{2}\left( \mathbf {x^{+}}\left( c_{\theta },\zeta ,\lambda \right) \right) \\&\quad =\frac{\left( -8c_{\theta }^{2}(c_{\theta }-3\zeta )-\frac{16c_{\theta } ^{3}\left( 2c_{\theta }^{2}+c_{\theta }(2-3\zeta )\lambda -3\zeta \lambda \right) }{\sqrt{6c_{\theta }^{2}(c_{\theta }+1)\zeta -4c_{\theta }^{3}}-2c_{\theta }^{2} }-\frac{16(c_{\theta }+1)c_{\theta }^{4}(2c_{\theta }-3\zeta )(c_{\theta }+3\zeta )}{\left( \sqrt{6c_{\theta }^{2}(c_{\theta }+1)\zeta -4c_{\theta }^{3} }-2c_{\theta }^{2}\right) ^{2}}\right) }{16c_{\theta }^{2}(3(c_{\theta }+1)\zeta -2c_{\theta })(c_{\theta }+1)^{-1}},\\&e_{1}\left( \mathbf {x^{-}}\left( c_{\theta },\zeta ,\lambda \right) \right) \\&\quad =\sqrt{6(c_{\theta }+1)\zeta -4c_{\theta }}\left( \frac{3}{3\zeta -2c_{\theta }}+\frac{4}{3(c_{\theta }+1)\zeta -2c_{\theta }}\right) \\&\qquad +\frac{3c_{\theta }}{3\zeta -2c_{\theta }}+1,\\&e_{2}\left( \mathbf {x^{-}}\left( c_{\theta },\zeta ,\lambda \right) \right) \\&\quad =\frac{\left( -8c_{\theta }^{2}(c_{\theta }-3\zeta )+\frac{16c_{\theta } ^{3}\left( 2c_{\theta }^{2}+c_{\theta }(2-3\zeta )\lambda -3\zeta \lambda \right) }{2c_{\theta }^{2}+\sqrt{6c_{\theta }^{2}(c_{\theta }+1)\zeta -4c_{\theta }^{3}} }-\frac{16(c_{\theta }+1)c_{\theta }^{4}(2c_{\theta }-3\zeta )(c_{\theta }+3\zeta )}{\left( 2c_{\theta }^{2}+\sqrt{6c_{\theta }^{2}(c_{\theta }+1)\zeta -4c_{\theta }^{3}}\right) ^{2}}\right) }{16c_{\theta }^{2}(3(c_{\theta }+1)\zeta -2c_{\theta })(c_{\theta }+1)^{-1}},\\&e_{1}\left( \mathbf {x^{1}}\left( c_{\theta },\zeta ,\lambda \right) \right) \\&\quad =\frac{(c_{\theta }+1)\left( \lambda ^{2}(3(c_{\theta }+1)\zeta -2c_{\theta })-18\zeta ^{2}\right) }{2(3(c_{\theta }+1)\zeta -2c_{\theta })(c_{\theta } \lambda +3\zeta )},\\&e_{2}\left( \mathbf {x^{1}}\left( c_{\theta },\zeta ,\lambda \right) \right) =\frac{\lambda (2(c_{\theta }+1)\lambda -c_{\theta }+1)-6\zeta }{2c_{\theta }\lambda +6\zeta }, \end{aligned}$$

and

$$\begin{aligned} e_{1}\left( \mathbf {x^{2}}\left( c_{\theta },\zeta ,\lambda \right) \right)&=-\frac{3\zeta (c_{\theta }(\lambda -1)+\lambda +1)+2c_{\theta }\lambda }{2(2\lambda +1)(3(c_{\theta }+1)\zeta -2c_{\theta })}\\&\quad +\frac{\Delta }{2(2\lambda +1)(3(c_{\theta }+1)\zeta -2c_{\theta })},\\ e_{2}\left( \mathbf {x^{2}}\left( c_{\theta },\zeta ,\lambda \right) \right)&=-\frac{\Delta }{2(2\lambda +1)(3(c_{\theta }+1)\zeta -2c_{\theta })}\\&\quad -\frac{3\zeta (c_{\theta }(\lambda -1)+\lambda +1)+2c_{\theta }\lambda }{2(2\lambda +1)(3(c_{\theta }+1)\zeta -2c_{\theta })}, \end{aligned}$$

where \(\Delta =\surd \left( 3\zeta \left( -8c_{\theta }(1+c_{\theta })+15(-1+c_{\theta })^{2}\zeta +72(1+c_{\theta })\zeta ^{2}\right) +\right. 2(-1+c_{\theta })\left( -2c_{\theta }(1+c_{\theta })+3(1+c_{\theta }(4+c_{\theta }))\zeta -27(1+c_{\theta })\zeta ^{2}\right) \lambda +\left( 4c_{\theta }\left( 4+c_{\theta }+4c_{\theta } ^{2}\right) -12(-2+c_{\theta })(1+c_{\theta })(-1+2c_{\theta })\zeta -63(1+c_{\theta })^{2}\zeta ^{2}\right) \lambda ^{2}+\left. 8\left( -1+c_{\theta }^{2}\right) (-2c_{\theta }+3(1+c_{\theta })\zeta )\lambda ^{3}\right) \).

In Figs. 7 and 8 we present the regions in the space of parameters \(\left\{ c_{\theta },\zeta \right\} \) where the stationary points of the dynamical system (53), (54) are attractors for \(\lambda =2\) and \(\lambda =0\).

7 Conclusions

In this work, we considered the extension of Einstein-aether theory in Weyl geometry. In particular, we considered the Einstein-aether action integral in the Weyl integrable geometry, where a geometric scalar field is introduced. The scalar field plays a significant role in the geometry since it defines the deviation of the Weyl affine connection from that of the Levi–Civita connection.

The action integral is considered in an isotropic and homogeneous (FLRW) background spacetime. We observe that the scalar which defines the Weyl affine connection is introduced into the gravitational field equation and it is dynamically coupled to the aether field. The scalar field is introduced two-fold from the Weyl Ricci scalar, and from the symmetric component of the covariant derivatives for the aether field; from the latter terms, the coupling then follows. This approach is an alternative way to create an Einstein-aether scalar field cosmological model. In the case of a spatially-flat vacuum FLRW spacetime, the field equations admit an exact solution where the scale factor is power-law and the parameter for the equation of state can describe acceleration.

We studied the cosmological models where a scalar field potential is included in the field equations, or dust fluid source contributes to the gravitational field equations, and when these two terms exist together. For these three additional systems, we study the dynamics and we find all the possible asymptotic behaviour of the field equations. We demonstrate our results by presenting some specific applications. We remark that the results of this work differ from the previous studies on the Einstein-aether scalar field cosmological models; however, it is clear that this model can describe an alternative inflationary behaviour. In a future work, we plan to study the stability of the small inhomogeneities for this scalar field model.