1 Introduction

After introducing the general theory of relativity [1], Einstein in 1919 tried to formulate the elementary particles interactions through general relativity [2]. He then considered an electromagnetic like matter field with energy–momentum tensor \(S_{\mu \nu }\) representing elementary particles [2]. Since this tensor is trace-less, he modified the Einstein equation to describe the new interaction as

$$\begin{aligned} R_{\mu \nu }-\frac{1}{4}Rg_{\mu \nu }=\kappa ^2 S_{\mu \nu }. \end{aligned}$$
(1)

Moreover, Einstein assumed that elementary particle field satisfies Maxwell’s equation. So taking the covariant derivative of the above equation yields

$$\begin{aligned} \frac{1}{4} \nabla _\mu R=-\frac{\kappa ^2}{c} F_{\mu \nu }j^\nu , \end{aligned}$$

where \(j^\nu \) is the electric current. It is evident that in the volume outside the elementary particles and in the absence of charges we obtain from the above equation that

$$\begin{aligned} R=R_0=constant. \end{aligned}$$

Einstein, also assumed that the standard equation of general relativity still holds, so we also have

$$\begin{aligned} R_{\mu \nu }-\frac{1}{2}Rg_{\mu \nu }+\Lambda g_{\mu \nu }=\kappa ^2T_{\mu \nu }, \end{aligned}$$
(2)

where \(T_{\mu \nu }\) is the energy–momentum tensor of the baryonic matter. In vacuum, from the above equation we obtain

$$\begin{aligned} R=4\Lambda , \end{aligned}$$

which by comparison with the result from the elementary particles equation we obtain

$$\begin{aligned} \Lambda =\frac{R_0}{4}. \end{aligned}$$
(3)

This shows that the cosmological constant can be considered as an integration constant from the elementary particles field equation [2]. By using Eqs. (1)–(3) one can obtain the so-called matter-geometry symmetric Einstein equation

$$\begin{aligned} R_{\mu \nu }-\frac{1}{4}Rg_{\mu \nu }=\kappa ^2\left( T_{\mu \nu }-\frac{1}{4}Tg_{\mu \nu }\right) , \end{aligned}$$
(4)

and

$$\begin{aligned} S_{\mu \nu }=T_{\mu \nu }-\frac{1}{4}Tg_{\mu \nu }, \end{aligned}$$

where \(S_{\mu \nu }\) as was mentioned before is the energy–momentum tensor of the elementary particles in the form of electromagnetic fields. For a dust universe we obtain a very interesting result from the above equations. In this case \(T^0_0=-\rho \) and then we obtain \(S^0_0=-\,3/4\rho \). This means that the Einstein model predicts that the ingredient of matter in the universe is \(75\%\) electromagnetic and \(25\%\) gravitational. This theory has also been introduced as the unimodular gravity, which can be formulated in different ways. In this theory the determinant of metric tensor is constrained to be a number or scalar density. As a consequence the cosmological constant appears as a constant of integration. The action of unimodular gravity with a fixed metric determinant i.e. \(\sqrt{-g}=\epsilon _0\) [3] is

$$\begin{aligned} S_{UG}=\int d^4x \left[ \sqrt{-g}\kappa ^2 R-\lambda \left( \sqrt{-g}-\epsilon _0\right) \right] +S_m, \end{aligned}$$
(5)

where \(\lambda \) is the Lagrange multiplier and \(S_m\) is the action of matter fields. It should be noted that the action \(S_{UG}\) is invariant under a restricted group of diffeomorphisms in which the determinant of metric tensor is unchanged. The field equation of this model is the same as Einstein–Hilbert equations with cosmological constant, together with the constraint \(\sqrt{-g}=\epsilon _0\). By assuming the conservation of the energy momentum tensor one can obtain \(\nabla _\mu \lambda =0\). Henneaux and Teitelboin introduced a unimodular action which is fully diffeomorphism invariant [4]. The action of this model is as follow

$$\begin{aligned} S_{HT}=\int d^4x \left[ \sqrt{-g}\kappa ^2 R-\lambda \left( \sqrt{-g}-\partial _\mu \tau ^\mu \right) \right] +S_m, \end{aligned}$$
(6)

where \(\tau ^\mu \) is a vector density. The equation of motion for the metric tensor is the same as Einstein–Hilbert field equations with cosmological constant but the determinant of the metric is constrained to be \(\sqrt{-g}=\partial _\mu \tau ^\mu \). The actions \(S_{UG}\) and \(S_{HT}\) are classically equivalent and by a change of coordinates they are related to each other. There is an alternative action which is fully diffeomorphism invariant [5] as

$$\begin{aligned} S_{DUG}=\int d^4x \sqrt{-g}\left[ \kappa ^2 R-\lambda +V^\mu \nabla _\mu \lambda \right] +S_m, \end{aligned}$$
(7)

where the vector field \(V^\mu \) is the Lagrange multiplier to keep constraint \(\nabla _\mu \lambda =0\). In this action the constraint on the value of the metric determinant is replaced by \(\nabla _\mu V^\mu =1\). By integrating by part of the action \(S_{HT}\), one can easily show that the two actions \(S_{HT}\) and \(S_{DUG}\) are equivalent at classical level. In [3] the authors use the path integral to show the differences between the above mentioned models of unimodular gravity at quantum level. However, the introduction of unimodular gravity in a more general way has been reviewed in [6,7,8].

Rastall [9] is the first person who considered the modified equation (4) with the baryonic energy–momentum tensor instead of elementary \(S_{\mu \nu }\) tensor. He also assumed a general coupling for the Ricci scalar

$$\begin{aligned} R_{\mu \nu }+\lambda Rg_{\mu \nu }=\kappa ^2T_{\mu \nu }. \end{aligned}$$

This equation predicts that the energy–momentum tensor is no longer conserved and there is a chance to transform directly to geometry \(\nabla _\nu T^{\mu \nu }\propto \nabla _\nu R\).

Rastall theory have been investigated vastly in the literature [10,11,12,13,14,15,16,17,18,19,20]. Also, many generalizations and modifications of the idea has been proposed. One of the most interesting of them is to consider a theory containing a non-minimal coupling between matter Lagrangian and geometry [21,22,23,24]. The action can then be written in the form

$$\begin{aligned} S=\int \sqrt{-g}d^4x\left( \kappa ^2R+f(R,\mathcal {L}_m)+\mathcal {L}_m\right) , \end{aligned}$$
(8)

where \(\mathcal {L}_m\) is the matter Lagrangian. This theory has a general property of the Rastall theory which is the non-conservation of the energy–momentum tensor. This causes the matter fields to be converted directly to geometry. Cosmological implications of this theory is vastly investigated [21,22,23,24]. Other generalizations of this idea includes non-minimal coupling between matter energy–momentum tensor and geometry such as f(RT) [25,26,27,28,29,30,31,32,33,34,35,36], \(f(R,T,R_{\mu \nu }T^{\mu \nu })\) [37, 38] gravity theories. Also, one can consider non-standard interactions between matter fields such as \(f(T,T_{\mu \nu }T^{\mu \nu })\) theories where T is the trace of energy–momentum tensor [39,40,41], or derivative matter couplings [42,43,44] where one considers interactions of the form \(\nabla _\mu \mathcal {L}_m\nabla ^\mu \mathcal {L}_m\).

In this paper, we are going to consider the cosmological implications of the Einstein dark energy theory introduced in [45], where the spirit of Einstein idea is putted together with the properties of Rastall gravity. In this way, we will consider an Einstein–Hilbert action coupled with a dark energy vector field denoted as \(\Lambda _\mu \) which is minimally coupled to geometry, but non-minimally coupled to the baryonic energy–momentum tensor through the interaction of the form \(\mathcal {L}_m\Lambda _\mu \Lambda ^\mu \). We also considered the Rastall theory by adding a term proportional to the trace of the energy–momentum tensor. It should be noted that the present theory is different from the unimodular gravity since in the Einstein dark energy theory the vector field does not constraint the metric field and it is an independent dynamical field. Also, since the vector field has a non-minimal coupling with matter Lagrangian, the matter fields do not conserve in this theory as opposed to the unimodular gravity.

In [45], the authors have shown that the theory can describe the late time accelerated expansion of the universe. However, the theory was not fully satisfactory with recent observational data. In this paper, by adding a non-minimal interaction term between matter fields and the dark energy vector field, we will analyze the cosmological implications in more details and show that the modified theory is in fact capable of explaining the recent observational data in both background and first order perturbation levels. We also estimate the values of the model parameters to obtain the best fit of the theory with experiments by two sets of data corresponding to the Hubble [46] and \(f\sigma _8\) [47] functions. We have also considered the dynamical system analysis of the model and show that the theory has three fixed points. The theory we are considering does not have a conservation of the energy–momentum tensor. This implies that the behavior of the energy density in dust/radiation dominated universes are not the same as in general relativity. As a result two of the fixed points of the Einstein dark energy model are the would be dust and radiation dominated fixed points which now behaves differently due the the presence of non standard matter couplings in the theory. The third fixed point corresponds to the de Sitter expansion of the universe. Similar to general relativity, this fixed point is stable. So, the theory can in principle explain the thermal history of the universe. Also, we will analyze the theory at the perturbative level and obtain the evolution equation corresponding to the growth rate of matter perturbations.

The structure of the paper is as follows. In the next section we review the Einstein dark energy model and obtain main dynamical equations of the model. In Sect. 3, we consider the background cosmological implications of the model and in Sect. 4 we investigate its dynamical analysis. In Sect. 5, the matter scalar perturbations on top of flat FRW universe is considered and at the end we will conclude the paper.

2 The model

In this section we present the field equations for the Einstein dark energy model, and derive some of its basic theoretical consequences [45]. Let us assume that the universe is filled with a cosmological dark energy vector field \(\Lambda _{\mu }\left( x^{\nu }\right) \). We define the dark energy strength tensor as

$$\begin{aligned} C_{\mu \nu }=\nabla _{\mu }\Lambda _{\nu }-\nabla _{\nu }\Lambda _{\mu }. \end{aligned}$$
(9)

The dark energy strength tensor identically satisfies the Maxwell type equations

$$\begin{aligned} \nabla _{\lambda }C_{\mu \nu }+\nabla _{\mu }C_{\nu \lambda }+\nabla _{\nu }C_{\lambda \mu }=0. \end{aligned}$$
(10)

We define the energy–momentum tensor \(T_{\mu \nu }\) of the baryonic matter fields as

$$\begin{aligned} T_{\mu \nu }=-\frac{2}{\sqrt{-g}} \frac{\partial \left( \sqrt{-g} {\mathcal {L}}_{m} \right) }{ \partial g^{\mu \nu }}, \end{aligned}$$
(11)

where \({\mathcal {L}}_{m}\) is the Lagrangian of the total (ordinary baryonic plus dark) matter. In the following, by T we denote the trace of the matter energy–momentum tensor.

The Einstein dark energy model is described by the following action

$$\begin{aligned} S&=\int d\,^4x\sqrt{-g}\,\Bigg [\kappa ^2\left( 1-\beta _1 \right) R+ \frac{\beta _2 }{2}\,T-\frac{1}{4}C_{\mu \nu }C^{\mu \nu } \nonumber \nonumber \\&\quad +\alpha \,{\mathcal {L}}_m\,\Lambda _{\mu }\Lambda ^\mu +V(\Lambda ^2)+{\mathcal {L}}_{m} \Bigg ] , \end{aligned}$$
(12)

where \(\beta _1\), \(\beta _2\) are two arbitrary dimensionless constants and \(\alpha \) is a coupling constant with mass dimension \(M^{-2}\) representing the interaction between matter and the dark energy vector \(\Lambda _\mu \). Also, the potential term V is an arbitrary function of \(\Lambda ^2=\Lambda _\mu \Lambda ^\mu \). In this paper, we will consider the constant potential corresponding to the cosmological constant, and also a power-law case.

The energy–momentum tensor \(S_{\mu \nu }\) of the dark energy field can be obtained by varying its kinetic term with respect to the metric, which gives

$$\begin{aligned} S_{\mu \nu }=C_{\mu \alpha }C_{\nu }^{~\alpha }-\frac{1}{4}g_{\mu \nu }C_{\alpha \beta }C^{\alpha \beta }, \end{aligned}$$
(13)

with the property \(S_{\mu }^{\mu }=0\).

By varying the gravitational action with respect to the metric tensor, it follows that the cosmological evolution of the universe in the presence of a vector type dark energy is described by the generalized Einstein gravitational field equations,

$$\begin{aligned}&\kappa ^2(1-\beta _1)G_{\mu \nu }-\frac{1}{2}S_{\mu \nu }\nonumber \\&\qquad +\alpha \left( {\mathcal {L}}_m\Lambda _\mu \Lambda _\nu -\frac{1}{2}\Lambda ^2T_{\mu \nu }\right) -\frac{1}{2}g_{\mu \nu }V+\Lambda _\mu \Lambda _\nu V^\prime \nonumber \\&\quad =\frac{1}{2}(1+\beta _2)T_{\mu \nu }-\frac{1}{2}\beta _2\,\left( {\mathcal {L}}_m-\frac{1}{2}T\right) g_{\mu \nu }, \end{aligned}$$
(14)

where prime denotes derivative with respect to the argument.

By varying the action (12) with respect to the vector potential, we obtain the equation

$$\begin{aligned} \nabla _{\nu }C^{\mu \nu }=2\Lambda ^\mu (V^\prime +\alpha {\mathcal {L}}_m). \end{aligned}$$
(15)

It should be noted that due to the non-minimal coupling between matter and geometry, the matter field is the source for the dark energy vector field.

By taking the divergence of the metric field equation (14) and using Eq. (15) one obtains the conservation equation of the energy–momentum tensor as

$$\begin{aligned} \nabla ^\mu T_{\mu \nu }&=\frac{1}{1+\beta _2+\alpha \Lambda ^2}\Bigg [\beta _2\nabla _\nu \left( {\mathcal {L}}_m-\frac{1}{2}T\right) \nonumber \\&\quad +\alpha ({\mathcal {L}}_mg_{\mu \nu }-T_{\mu \nu })\nabla ^\mu (\Lambda ^2)\Bigg ]. \end{aligned}$$
(16)

It can be seen from the above equation that there are two sources for the non-conservation of the energy–momentum tensor. The first one is due to the presence of the trace of the energy–momentum tensor T and the second one comes from the non-minimal coupling between matter and the dark energy vector potential in the action. In the case \(\alpha =0=\beta _2\) the energy–momentum tensor becomes conserved. We will defined a vector field

$$\begin{aligned} f_\nu&\equiv \frac{1}{1+\beta _2+\alpha \Lambda ^2}\Bigg [\beta _2\nabla _\nu \left( {\mathcal {L}}_m-\frac{1}{2}T\right) \nonumber \\ {}&\quad +\alpha ({\mathcal {L}}_mg_{\mu \nu }-T_{\mu \nu })\nabla ^\mu (\Lambda ^2)\Bigg ], \end{aligned}$$
(17)

which is the right-hand side of Eq. 16 and represents the amount of non-conservation of the energy–momentum tensor. In the case of \(f_\nu =0\), the energy–momentum tensor becomes conserved.

3 Cosmological implications

Let us consider the flat FRW universe with conformal time t

$$\begin{aligned} ds^2=a^2\eta _{\mu \nu }dx^\mu dx^\nu , \end{aligned}$$
(18)

where \(a=a(t)\) is the scale factor. The Hubble parameter can be defined as \(H=\dot{a}/a\), where dot represents derivative with respect to the conformal time.

For the dark energy vector field, we assume that only the temporal component is non-vanishing

$$\begin{aligned} \Lambda _\mu = a\big [\Lambda _0(t),{\vec {0}}\big ], \end{aligned}$$
(19)

which is dictated by the isotropy and homogeneity of the FRW space time in the Cartesian coordinates. We also assume that the universe is filled with a perfect fluid with Lagrangian density \({\mathcal {L}}_m=-\rho \) and energy–momentum tensor

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

where \(\rho \) is the energy density and p is the thermodynamics pressure.

The Friedmann and Raychaudhuri equations can be obtained from (14) as

$$\begin{aligned}&\frac{6\kappa ^2}{a^2}(1-\beta _1)H^2=(1+\alpha \Lambda _0^2)\rho \nonumber \\&\qquad +\frac{1}{2}\beta _2(\rho -3p)-(V+2\Lambda _0^2V^\prime ), \end{aligned}$$
(21)
$$\begin{aligned}&\frac{2\kappa ^2\left( 1-\beta _1\right) }{a^2}\left( H^2+2 \dot{H}\right) \nonumber \\&\quad =(\alpha \Lambda _0^2-1)p-\frac{1}{2}\beta _2(\rho +5p)-V.\end{aligned}$$
(22)

The field equation of the vector field is

$$\begin{aligned} V^\prime =\alpha \rho . \end{aligned}$$
(23)

From the above equation, one can see that in the case of a constant potential \(V=const.\), the non-minimal coupling between the matter and the dark energy vector field should vanish.

The non-conservation equation of the matter field (16) is reduced to

$$\begin{aligned}&\bigg (1+\frac{1}{2}\beta _2-\alpha \Lambda _0^2\bigg )\dot{\rho }-\frac{3}{2}\beta _2\dot{p}\nonumber \\&\quad +3(1+\beta _2-\alpha \Lambda _0^2)H(\rho +p)=0. \end{aligned}$$
(24)

Now, let us assume a specific form of the potential as

$$\begin{aligned} V(\Lambda ^2)=-\beta _3\left( -\frac{\Lambda ^2}{\kappa ^2}\right) ^\eta , \end{aligned}$$
(25)

where \(\eta \) is a dimensionless constant and \(\beta _3\) is a constant with mass dimension \(M^4\).

In the case \(\eta =0\), one has \(V=-\beta _3\) which mimics the cosmological constant. So, we define the modified cosmological constant in this model as \(\lambda =\beta _3/2\kappa ^2\). In the case of \(\eta =0\), \(\lambda \) coincides with the standard cosmological constant. However, the value of \(\lambda \) will differ from the cosmological constant for \(\eta \ne 0\).

Let us assume that the universe is filled with pressure-less dust and radiation. The energy density and pressure becomes

$$\begin{aligned} \rho =\rho _m+\rho _r,\quad p=\frac{1}{3}\rho _r, \end{aligned}$$
(26)

where m/r denotes dust/radiation respectively.

Defining the following dimensionless parameters

$$\begin{aligned}&\tau =H_0 t,\quad H=H_0 h, \quad \bar{\rho }_i=\frac{\rho _i}{6\kappa ^2H_0^2},\nonumber \\&\quad \beta _3=6\kappa ^2H_0^2\Omega _\lambda ,\quad \alpha _1=\alpha \kappa ^2,\quad \Lambda _0=\kappa \Lambda _1, \end{aligned}$$
(27)

where \(H_0\) is the current Hubble parameter, one can write the metric field equations as

$$\begin{aligned}&(1-\beta _1)h^2=a^2\Bigg [\Omega _\lambda +\left( 1+\frac{1}{2}\beta _2\right) \bar{\rho }_m\nonumber \\&\quad +\alpha _1\Lambda _1^2(\bar{\rho }_m+\bar{\rho }_r)+(2\eta -1)\Omega _\lambda \Lambda _1^{2\eta }\Bigg ], \end{aligned}$$
(28)
$$\begin{aligned}&(1-\beta _1)(h^2+2h^\prime )=-a^2\Bigg [\frac{3}{2}\beta _2\bar{\rho }_m-3\Omega _\lambda \Lambda _1^{2\eta }\nonumber \\&\quad +\big (1+4\beta _2-\alpha _1\Lambda _1^2\big )\bar{\rho }_r\Bigg ], \end{aligned}$$
(29)

where prime here denotes derivative with respect to the dimensionless time \(\tau \). The vector field equation can be written as

$$\begin{aligned} \eta \Omega _\lambda \Lambda _1^{2(\eta -1)}=\alpha _1(\bar{\rho }_m+\bar{\rho }_m). \end{aligned}$$
(30)

As we have discussed before, in the case of \(\eta =0\) which corresponds to the case of cosmological constant, one should impose \(\alpha _1=0\). So \(\Lambda _1\) does not contribute to the background cosmological evolution of the universe. In the case \(\eta \ne 0\), one can obtain the value of \(\Lambda _1\) from the vector field equation (30) as

$$\begin{aligned} \Lambda _1=\left[ \frac{\alpha _1(\bar{\rho }_m+\bar{\rho }_r)}{\eta \Omega _\lambda }\right] ^{\frac{1}{2(\eta -1)}}. \end{aligned}$$
(31)

The baryonic matter conservation equation can be written in dimensionless coordinates as

$$\begin{aligned}&\Big [(2+\beta _2-2\alpha _1\Lambda _1^2)\bar{\rho }_m^\prime +6(1+\beta _2-\alpha _1\Lambda _1^2)h\bar{\rho }_m\Big ]\nonumber \\&\quad +\Big [2(1-\alpha _1\Lambda _1^2)\bar{\rho }_r^\prime +8(1+\beta _2-\alpha _1\Lambda _1^2)h\bar{\rho }_r\Big ]=0. \end{aligned}$$
(32)

From the structure of the above equation, we will assume that each bracket in (32) vanishes independently and as a result we have two conservation equations for dust and radiation as

$$\begin{aligned} \bar{\rho }_m^\prime +6h\left[ \frac{1+\beta _2-\alpha _1\Lambda _1^2}{2+\beta _2-2\alpha _1\Lambda _1^2}\right] \bar{\rho }_m=0, \end{aligned}$$
(33)
$$\begin{aligned} \bar{\rho }_r^\prime +4h\left[ \frac{1+\beta _2-\alpha _1\Lambda _1^2}{1-\alpha _1\Lambda _1^2}\right] \bar{\rho }_r=0. \end{aligned}$$
(34)

It should be noted that in the case \(\beta _2=0=\alpha _1\), the above equations becomes standard conservation equations for dust and radiation.

Fig. 1
figure 1

The Hubble and deceleration parameters as a function of the redshift z for three different values for the constant \(\eta =-\,0.02\) (dot-dashed), \(\eta =0\) (dashed) and \(\eta =0.1\) (dotted). We use the best fit values for the model parameters presented in Table 2. For the Hubble parameter, we have also plotted te experimental data together with their errors [46]. The \(\Lambda \)CDM curve is depicted as a red solid curve

Table 1 The fixed points of the dynamical system (43) and (44) together with eigenvalues and also their behaviors

Let us define redshift parameter as

$$\begin{aligned} 1+z=\frac{1}{a}. \end{aligned}$$
(35)

Derivatives with respect to \(\tau \) can be converted to the derivatives with respect to the redshift as

$$\begin{aligned} \frac{d}{d\tau }=-(1+z)h(z)\frac{d}{dz}. \end{aligned}$$
(36)

Using Eqs. (28) and (30) and transforming to the redshift coordinates, one can obtain the dimensionless Hubble parameter as

$$\begin{aligned} h(z)= & {} \frac{1}{\sqrt{1-\beta }}(1+z)^{-1}\Bigg [\left( 1+\frac{1}{2}\beta _2\right) \bar{\rho }_m(z)+\bar{\rho }_r(z)\nonumber \\&\quad +(1-\eta )\left[ \frac{1}{\Omega _\lambda }\left( \frac{\alpha _1(\bar{\rho }_m(z)+\bar{\rho }_r(z))}{\eta }\right) ^\eta \right] ^{\frac{1}{\eta -1}}\Bigg ]^{1/2}.\nonumber \\ \end{aligned}$$
(37)

The deceleration parameter in terms of redshift can be written as

$$\begin{aligned} q=\left( 1+z\right) \frac{d \ln h}{dz}. \end{aligned}$$
(38)

Noting that \(h(0)=1\), \(\bar{\rho }_m(0)=\Omega _{m0}=0.305\) and \(\bar{\rho }_r(0)=\Omega _{r0}=5.3\times 10^{-5}\) in redshift coordinates [48], one can see from Eq. (37) that the modified cosmological constant density parameter can be expressed by other parameters as

$$\begin{aligned} \Omega _\lambda =\left( \frac{\eta }{\alpha _1}\right) ^\eta \left[ \frac{2(\eta -1)(\Omega _{m0}+\Omega _{r0})^{\frac{\eta }{\eta -1}}}{2(\beta _1-1)+(2+\beta _2\Omega _{m0}+2\Omega _{r0})}\right] ^{\eta -1}. \end{aligned}$$
(39)

Also, we should note that on top of FRW space-time, the temporal component of the vector field \(f_\nu \) is non-vanishing and in dimensionless coordinates is given by

$$\begin{aligned} f_0=\frac{3\beta _2(1+z)h(\bar{\rho }_m^\prime +2\bar{\rho }_r^\prime ) }{1+\beta _2-\alpha _1\Lambda _1^2}. \end{aligned}$$
(40)

In Fig. 1 we have plotted the evolution of the Hubble parameter h together with the deceleration parameter q as a function of redshift. We have assumed three different values for the constant \(\eta =-\,0.02,0,0.1\). It should be emphasized that in Sect. 5, we will obtain the best fit values of the model parameters \(\beta _1\) and \(\beta _2\) using two independent data sets of Hubble parameter and \(f\sigma _8\). The best fit values of the parameters, together with their \(1\sigma \) and \(2\sigma \) confidence intervals are shown in Table 1, for three values of \(\eta =-\,0.02,0,0.1\). In the Fig. 1, we have used the best fit values. It can be shown from the figure that the Einstein dark energy model can explain the observational data on the Hubble parameter very well. In the case \(\eta \ne 0\), one can see that the Hubble parameter becomes larger than \(\Lambda \)CDM value for redshifts greater than \(z\sim 1.5\). In the case of cosmological constant \(V(\Lambda ^2)=-\beta _3\) however, the Hubble parameter is very close to the \(\Lambda \)CDM curve. However, the EDE prediction of the Hubble parameter for redshifts \(z>1.5\) is a little smaller than \(\Lambda \)CDM value. This shows that the size of the universe is smaller for non-vanishing values of \(\eta \). The evolution of the deceleration parameter shows that the universe have more deceleration compared to the \(\Lambda \)CDM model at redshift larger than \(z\sim 1.5\). It should be mentioned that the Einstein dark energy model with a cosmological constant behaves a little different from the other cases with non-vanishing \(\eta \). This is due to the fact that in the case of \(\eta =0\), the non-minimal coupling between matter and dark energy vector field vanishes, which makes the universe to gain more acceleration.

In Fig. 2, we have plotted the evolution of the density abundance \(\Omega _m\), defined as

$$\begin{aligned} \Omega _m=\frac{\bar{\rho }_m a^2}{h^2}. \end{aligned}$$
Fig. 2
figure 2

The evolution of the matter density abundance as a function of the redshift z for three different values for the constant \(\eta =-\,0.02\) (dot-dashed), \(\eta =0\) (dashed) and \(\eta =0.1\) (dotted). We use the best fit values for the model parameters presented in Table 2. The \(\Lambda \)CDM curve is depicted as a red curve

As can be seen from the figure, the baryonic matter density becomes larger than the conservative \(\Lambda \)CDM model. The difference can be seen as an amount of non-conservation of the energy momentum tensor. This can also be seen from Eq. 40 which we have plotted in Fig. 3.

Fig. 3
figure 3

The evolution of the non-conservation of the energy–momentum tensor \(f_0\) as a function of the redshift z for three different values for the constant \(\eta =-\,0.02\) (dot-dashed), \(\eta =0\) (dashed) and \(\eta =0.1\) (dotted). We use the best fit values for the model parameters presented in Table 2

For redshifts greater than \(z>0.2\) the non-conservation of the energy–momentum tensor \(f_0\) becomes non-zero which causes the matter density abundance to behave differently from the \(\Lambda \)CDM theory.

We have also plotted the temporal component of the dark energy vector field in Fig. 4 for non-vanishing values of the parameter \(\eta \). In the case of vanishing \(\eta \), this quantity vanishes as we have discussed earlier. It should be noted that from the structure of the equations, the coupling constant \(\alpha \) does not appear alone in the field equations. The only appearance of this constant is in the expression for the temporal component of the dark energy vector field (31). So, in order to plot the temporal component of the dark energy vector, we have modulated \(\Lambda _1\) by \((\alpha _1/\eta )^{1/2}\). It can be seen from the figure that the dark energy vector field tends to zero as the redshift increases. The maximum value of the vector occurs at present time with \(z=0\). Also, it should be mentioned that the qualitative behavior of the vector field is the same for different but non-vanishing values of \(\eta \).

4 Dynamical system analysis

Let us rewrite the Friedmann equation (37) in the form

$$\begin{aligned}&\left[ \Omega _\lambda \left( \frac{\eta }{\alpha _1}\right) ^{\eta }\frac{a^2}{h^2}\right] ^{\frac{1}{1-\eta }}=\frac{1}{\eta -1}\left[ \frac{a^2}{h^2}(\bar{\rho }_m+\bar{\rho }_r)\right] ^{\frac{\eta }{\eta -1}}\nonumber \\&\quad \times \left[ \beta _1-1+\left( 1+\frac{1}{2}\beta _2\right) \left( \frac{\bar{\rho }_m a^2}{h^2}\right) +2\left( \frac{\bar{\rho }_r a^2}{h^2}\right) \right] . \end{aligned}$$
(41)

The above equation suggests that the theory could be analyzed by two dynamical variables

$$\begin{aligned} \Omega _m=\frac{\bar{\rho }_m a^2}{h^2},\quad \Omega _r=\frac{\bar{\rho }_r a^2}{h^2}, \end{aligned}$$
(42)

which are the standard dust and radiation density abundances. Using the conservation equations (33) and (34), one can write the dynamical system of the model as

$$\begin{aligned} \Omega _m^\prime&=\frac{\Omega _m}{2}\Bigg (-2+f(\Omega _m,\Omega _r)\!+\!\frac{6\beta _2(\eta -1)(\Omega _m+\Omega _r)}{(2+\beta _2)(\Omega _m+\Omega _r)-\eta (2-2\beta _1+\beta _2\Omega _r)}\Bigg ), \end{aligned}$$
(43)
$$\begin{aligned} \Omega _r^\prime&=\frac{\Omega _r}{2}\Bigg (-4+f(\Omega _m,\Omega _r)+\frac{6\beta _2(\eta -1)(\Omega _m+\Omega _r)}{2)(\Omega _m+\Omega _r)-\eta (2-2\beta _1-\beta _2\Omega _m)}\Bigg ), \end{aligned}$$
(44)

where prime represents derivative with respect to \(\ln a\) and we have defined

$$\begin{aligned} f(\Omega _m,\Omega _r)=\frac{2(\Omega _m+\Omega _r)\Big (-2+2\beta _1+3(1+\beta _2)\Omega _m+4(1+\beta _2)\Omega _r\Big )-(\Omega _m+2\Omega _r)\Big (2-2\beta _1+2\beta _2\Omega _m+4\beta _2\Omega _r\Big )}{(\beta _1-1)(\eta -1)(\Omega _m+\Omega _r)}. \end{aligned}$$
(45)
Fig. 4
figure 4

The evolution of the modulated temporal component of the dark energy vector \(\Lambda _0\) as a function of the redshift zfor two non-vanishing values for the constant \(\eta =-\,0.02\) (dot-dashed) and \(\eta =0.1\) (dotted). We use the best fit values for the model parameters presented in Table 2

Fig. 5
figure 5

The phase space portrait of the dynamical system (43) and (44). We have plotted the \((\Omega _m,\Omega _r)\) planes for three different values of \(\eta =-\,0.02\) (left), \(\eta =0\) (center) and \(\eta =0.1\) (right). The corresponding fixed points are also shown

The model is a two dimensional autonomous dynamical system. In order to determine the behavior of the universe at the fixed points, we define an effective equation of state parameter

$$\begin{aligned} \omega _{eff}\equiv -\frac{1}{3}-\frac{2}{3}\frac{h^\prime }{h^2}=\frac{3(\beta _1-1)+(3\Omega _m+4\Omega _r)(1+\beta _2)}{3(\eta -1)(\beta _1-1)}+\frac{2\eta (\beta _1-1)\Omega _r-\eta \beta _2(\Omega _m+2\Omega _r)(3\Omega _m+4\Omega _r)}{6(\eta -1)(\beta _1-1)(\Omega _m+\Omega _r)}. \end{aligned}$$
(46)

The dynamical system (43) and (44) has three fixed points which we have summarized in Table 1. The fixed point M behaved like a matter dominated universe if \(\beta _2=0\). So, this point behaves a little different from the standard matter dominated epoch in general relativity. This also can be seen from the conservation equation (33). This fixed point is a saddle point. In Fig. 5, we have plotted the stream plot of the dynamical system (43) and (44) for three different values of \(\eta =-\,0.02,0,0.1\). We have also shown the fixed points in the figures.

The fixed point R behaves like a radiation dominated fixed point if \(\beta _2=0\). This fixed point is an unstable fixed point which plays a role of an approximate radiation dominated phase in our model. Remembering that the dynamical variables \(\Omega _m\) and \(\Omega _r\) are positive, one can see that these two fixed points behaves as unstable nodes in the positive triangle of the phase portrait.

The last fixed point \(\Lambda \) has an effective equation of state parameter \(\omega _{eff}=-\,1\) and is a de Sitter fixed point. This node is in fact a fixed line as can be seen in Fig. 5 as a solid black curve. In the case of \(\eta =0\), only the point \((\Omega _m,\Omega _r)=(0,0)\) lies in the positive quarter. This point is in fact the stable de Sitter fixed point of the standard \(\Lambda \)CDM theory. In the case of \(\eta =0.1\) we have a fixed line in the positive quarter and one can see from the figure that all the curves end up at this line. We should note that all the points in the line is in fact a fixed point. As a result in this case the phase space is smaller than that of general relativity since the curves can not escape the fixed line to end up at (0, 0) point. The case of \(\eta =-\,0.02\) is a little different since no points of the fixed line lie in the positive quarter of the phase space. In this case all the streams in this quarter will end up at the origin which is not a fixed point but as one can see from the figure that is very close to it.

In summary, in all three cases, the evolution of the universe can be started from the radiation dominated fixed point which continue to the matter dominated fixed point and then ends up at a stable de Sitter epoch. As a result the thermal history of the universe can be explained in this model.

5 Matter perturbation of the model

Table 2 The best values of the model parameters \(\beta _1\) and \(\beta _2\), together with the best values of \(\sigma _8^0\) and also the initial condition \(\xi \), for three different values of the parameter \(\eta =-\,0.02,0,0.1\)

In this section we will consider the growth of matter perturbation of the Einstein dark energy model. We will consider the scalar perturbations of the field equations in the Newtonian gauge. In this gauge, the scalar perturbation variables E, B vanishes and the perturbed conformal FRW universe can be written as

$$\begin{aligned} ds^2=a^2(t)\Big [-(1+2\varphi )dt^2+(1-2\psi )d\vec {x}^2\Big ], \end{aligned}$$
(47)

where \(\varphi \) and \(\psi \) are the Bardeen potentials. The perturbed energy momentum tensor is defined as

$$\begin{aligned} \delta T^0_0&=-\delta \rho \equiv -\rho \delta ,\nonumber \\ \delta T^0_i&=(1+w)\rho \partial _i v,\quad \nonumber \\ \delta T^i_j&=\delta ^i_jc_s^2\rho \delta , \end{aligned}$$
(48)

where, \(\delta \) is the matter density contrast defined as \(\delta =\delta \rho /\rho \), \(\rho \) is the background value of the energy density and v is the scalar mode of the velocity perturbation. Also, we have defined the sound speed as \(\delta p/\delta \rho =c_s^2\). The equation of state parameter of the background matter field can be given as \(p/\rho =w\). In the following, we will assume that the perturbed and unperturbed matter content of the universe have the equations of motion of the form \(c_s^2=0=w\).

The scalar mode of perturbed dark energy vector field is taken as

$$\begin{aligned} \Lambda _\mu = a[\Lambda _0+{\mathcal {M}},\,-\partial _i {\mathcal {N}} ] \end{aligned}$$
(49)

where \(\Lambda _0=\Lambda _0(t)\) is the non vanishing background component of the vector field \(\Lambda _\mu \).

The conservation equation (16) up to first order in perturbations can be written as

$$\begin{aligned}&2(2+\beta _2-2\alpha \Lambda _0^2)(1+\beta _2-2\alpha \Lambda _0^2)(\theta -3\dot{\psi })\nonumber \\&\quad -12\beta _2\alpha \Lambda _0 H(\Lambda _0\varphi -\mathcal {M})+(2+\beta _2-2\alpha \Lambda _0^2)^2\dot{\delta }=0, \end{aligned}$$
(50)

and

$$\begin{aligned}&(2+\beta _2-2\alpha \Lambda _0^2)\Big [k^2(4\alpha \Lambda _0\mathcal {M}-\beta _2\delta )+\alpha \Lambda _0\dot{\Lambda }_0\theta \Big ]\nonumber \\&\quad +2(2+\beta _2-2\alpha \Lambda _0^2)(1+\beta _2\alpha \Lambda _0^2)(\dot{\theta }-k^2\varphi )\nonumber \\&\quad -4H\Big [\beta _2^2-(1-\alpha \Lambda _0^2)^2\Big ]\theta =0, \end{aligned}$$
(51)

where we have defined \(\theta =\nabla _i\nabla ^i v\) and Fourier transformed with wave vector \(\vec {k}\). Note that in order to simplify the above equations, we have used the background conservation equation.

From now on, we will work in the sub-horizon limit where the wave number is much greater than the Hubble parameter \(k\gg H\). The above equations can be combined to eliminate the variable \(\theta \) and we obtain the evolution equation of the matter density contrast, which can be simplified in the sub-horizon limit as

$$\begin{aligned}&(2+\beta _2-2\alpha \Lambda _0^2)\ddot{\delta }+2\Big [(1-\beta _2-\alpha \Lambda _0^2) H-2\alpha \Lambda _0\dot{\Lambda }_0\Big ]\dot{\delta }\nonumber \\&\quad +2k^2\Big [(1+\beta _2+\alpha \Lambda _0^2)\varphi +\beta _2\delta -2\Lambda _0\mathcal {M}\Big ]=0. \end{aligned}$$
(52)

In order to obtain the relations of the variables \(\mathcal {M}\) and \(\varphi \), we will use the Einstein and also the dark energy vector field in the sub-horizon limit.

From non-diagonal components of the Einstein equation, one obtains \(\psi =\varphi \). Also the (0) component of the vector field equation in sub-horizon limit reads

$$\begin{aligned} \mathcal {M}=\dot{\mathcal {N}}-H\mathcal {N}. \end{aligned}$$
(53)

The (i) component of the vector field equation can be simplified in the sub-horizon limit to

$$\begin{aligned} H\dot{\mathcal {N}}+\left[ \beta _3\eta \left( \frac{\Lambda _0}{\kappa }\right) ^{2\eta }\Lambda _0^{-2}a^2+H^2-\alpha a^2\rho \right] \mathcal {N}=0, \end{aligned}$$
(54)

which gives the dynamical equation for the variable \(\mathcal {N}\). The (00) component of the Einstein field equation in the sub-horizon limit can be written as

$$\begin{aligned} 8(1-\beta _1)\kappa ^2 k^2\varphi =(2+\beta _2+2\alpha \Lambda _0^2)\rho a^2\delta . \end{aligned}$$
(55)

Now, substituting the variable \(\varphi \) from the Einstein equation (55) into equation (52), one obtains

$$\begin{aligned}&(2+\beta _2-2\alpha \Lambda _0^2)\ddot{\delta }+2\Big [(1-\beta _2-\alpha \Lambda _0^2)H-2\alpha \Lambda _0\dot{\Lambda }_0\Big ]\dot{\delta }\nonumber \\&\quad +\left[ 2\beta _2 k^2+\frac{(1+\beta _2+\alpha \Lambda _0^2)(2+\beta _2+2\alpha \Lambda _0^2)}{2(\beta _1-1)\kappa ^2}a^2\rho \right] \delta \nonumber \\&\quad +4\kappa ^2\alpha \Lambda _0(\dot{\mathcal {N}}-H\mathcal {N})=0. \end{aligned}$$
(56)

Making Eq. (54) dimensionless and substituting the background variables from Eq. (31), one obtains

$$\begin{aligned} 4h(z)^2\Big [(1+z)\mathcal {N}(z)^\prime -\mathcal {N}(z)\Big ]=0, \end{aligned}$$
(57)

where we have transformed to the redshift coordinates and prime represents derivative with respect to z. The above equation has a solution

$$\begin{aligned} \mathcal {N}(z)=c_1(1+z), \end{aligned}$$
(58)

where \(c_1\) is an integration constant. One can see that the \(\mathcal {N}\) dependency in Eq. (56) disappears and one obtains the evolution equation of the density contrast in dimensionless form as

$$\begin{aligned}&(2+\beta _2-2\alpha _1\Lambda _1^2)\delta ^{\prime \prime }\nonumber \\&\quad +2\Big [(1-\beta _2-\alpha _1\Lambda _1^2)h-2\alpha _1\Lambda _1\Lambda _1^\prime \Big ]\delta ^\prime \nonumber \\&\quad +\left[ 2\beta _2 \gamma ^2+\frac{3(1+\beta _2+\alpha _1\Lambda _1^2)(2+\beta _2+2\alpha _1\Lambda _1^2)}{2(\beta _1-1)}a^2\bar{\rho }_m\right] \delta =0. \end{aligned}$$
(59)

It should be noted that in the case of vanishing \(\beta _1\), \(\beta _2\) and \(\alpha \), the above equation reduces to the standard equation of the matter density contrast. In the above equation, we have also defined \(\gamma =k/H_0\), which is the dimensionless counterpart of the wave number.

Let us solve the equation governing the evolution of the density contrast. This should be solved together with background equations (31), (33) and (37). We will use a generalized \(\Lambda \)CDM initial conditions in deep matter dominated era in which

$$\begin{aligned} \frac{d\delta }{d\ln a}|_{z_\star }=\xi \delta |_{z_\star }, \end{aligned}$$
(60)

where \(z_\star \) is some point in the deep matter dominated era which we will assume to be \(z_\star =7.1\). Also \(\xi \) is a constant which determines deviation from \(\Lambda \)CDM model. The case \(\xi =1\) corresponds to the \(\Lambda \)CDM model.

Fig. 6
figure 6

The evolution of \(f\sigma _8(z)\) as a function of the redshift z for three different values for the constant \(\eta =-\,0.02\) (dot-dashed), \(\eta =0\) (dashed) and \(\eta =0.1\) (dotted). We use the best fit values for the model parameters presented in Table 2. The \(\Lambda \)CDM curve is depicted as a red curve

To compare the Einstein dark energy model with observational data, we will use two independent data sets on the Hubble parameter in the redshift range \(z\sim (0,2)\) [46] and also the observational data on \(f\sigma _8\) [47] which is defined as

$$\begin{aligned} f \sigma _8 \equiv \sigma _8(z)f(z). \end{aligned}$$
(61)

Here \(\sigma _8(z)=\sigma _8^0\, \delta (z)\). The constant \(\sigma _8^0\) is model dependent. The growth rate of matter perturbations is defined as

$$\begin{aligned} f=\frac{d \ln \delta }{d\ln a}=-(1+z)\frac{\delta '}{\delta }. \end{aligned}$$
(62)

We estimate the values of \(\sigma _8^0\), \(\xi \) and also the model parameters \(\beta _1\) and \(\beta _2\) by maximizing the likelihood function defined as

$$\begin{aligned} {\mathcal {L}}={\mathcal {L}}_{0}e^{-\chi ^2/2}, \end{aligned}$$
(63)

where \({\mathcal {L}}_0\) is the normalization constant and the quantity \(\chi ^2\) in our case is given by

$$\begin{aligned} \chi ^2&=\chi _H^2+\chi _{f\sigma _8}^2\nonumber \\ {}&=\sum _i\left( \frac{H_{i,obser}-H_{i,theory}}{\sigma _i}\right) ^2\nonumber \\ {}&~~~~+\sum _j\left( \frac{f\sigma _{8j,obser}-f\sigma _{8j,theory}}{\sigma _j}\right) ^2, \end{aligned}$$
(64)

where \(H_{i,theory}\) and \(f\sigma _{8j,theory}\) are the theoretical values for the observables \(H_{i,obser}\) and \(f\sigma _{8j,obser}\) and \(\sigma _i\) is the error of the ith data. The total likelihood function is the multiplication of individual likelihoods of the two sets since the data sets we are using here are independent.

In Table 2, we have summarized the best fit values of the parameters \(\beta _1\), \(\beta _2\), \(\xi \) and \(\sigma _8^0\) together with their \(1\sigma \) and \(2\sigma \) confidence intervals for three different values of the parameter \(\eta =-\,0.02,0,0.1\). We have used the best fit values of Table 2 to plot the figures. It can be seen from the table that the value of the parameter \(\beta _1\) at best fit and also up to \(2\sigma \) confidence level is negative. Also, the value of the parameter \(\beta _2\) is positive at best fit value. In Fig. 6, we have plotted the evolution of the function \(f\sigma _8\) as a function of redshift for three different values of the parameter \(\eta =-\,0.02,0,0.1\). We have also plotted the \(\Lambda \)CDM curve as a red solid line. One can see from the figure that the evolution of the function \(f\sigma _8\) for \(\eta =0\), is very similar to the \(\Lambda \)CDM theory. The difference can be traced back to the presence of the term proportional to \(\beta _2\). However, for non-vanishing cases of \(\eta \) the behavior of this function differs significantly from the \(\Lambda \)CDM theory. However, all cases satisfied observational data. The present value of \(f\sigma _8\) for positive values of the parameter \(\eta \) is greater than that of the \(\Lambda \)CDM value and for negative values of \(\eta \) its value becomes smaller. One can then see that for redshifts greater than 0.2, non-vanishing values of \(\eta \) fits very well with observations. In summary, it sees that more data would be needed to decide the best model which fits the data.

6 Conclusion and final remarks

In this paper we have considered cosmological implications of the Einstein dark energy model. This model consists of an Einstein–Hilbert Lagrangian with a coupling proportional to the trace of the energy–momentum tensor, coupled to a dark energy vector field. The vector field has a non-minimal coupling with matter fields, which make the energy–momentum tensor non-conservative. This non-conservation of the energy–momentum tensor results in creation of matter out of geometry. The rate of such a creation is calculated in [45]. We have also considered a power-law potential term for dark energy vector field.

We have obtained the best fit values of the model parameters by using two sets of independent data from Hubble parameter and also the function \(f\sigma _8\). At redshifts smaller than \(z\sim 3\), the behavior of the Hubble parameter is very similar to the \(\Lambda \)CDM model for small values of the parameter \(\eta \). For large values of this parameter, the behavior of the Hubble parameter is very different and one can not find a best fit with observational data. As a result larger values of \(\eta \) is ruled out by observations and so we have considered small values in this paper. Despite the behavior of the Hubble parameter the matter density abundance behaves differently from \(\Lambda \)CDM model. At larger redshifts, the matter density abundance is larger than the \(\Lambda \)CDM value implying that there are more matter present at those redshifts. However, the matter density decreases more rapid than \(\Lambda \)CDM model implying that the present values of the matter density abundance is the same as \(\Lambda \)CDM value. This shows that the rate of changing matter content of the universe to curvature is getting smaller at late times. This can also be seen from the evolution of the function \(f_0\) in Fig. 4, since the rate of creation of matter is proportional to this function [45].

Dynamical system analysis of the theory can also show this behavior. As we have discussed in this paper, the theory has three fixed points, one of them is an exact de Sitter node. However, We have two fixed points which behave a little different from \(\Lambda \)CDM matter and radiation nodes, since the constant \(\beta _2\) is small. The difference is directly related to the term \(\beta _2 T\) in the action which comes from Rastall’s idea. This makes the theory non-conservative and consequently the fixed points become different. However, the thermal history of the universe can be achieved in this theory since all the required fixed points are present.

We have also considered the first order perturbation analysis of the model and obtain the growth rate of matter density contrast in the sub-horizon limit. The differential equation governing the behavior of matter density contrast is affected by both Rastall’s term \(\beta _2T\) and also by non-minimal coupling term \(\mathcal {L}_m\Lambda _\mu \Lambda ^\mu \) term. Since the theory is different from \(\Lambda \)CDM theory at every times, we have modified the initial conditions on the density contrast to cover this new model. From the best fit values presented in Table 2, one can see that at deep matter dominated epoch, the rate of change of \(\delta \) is bigger than that of \(\Lambda \)CDM theory by about \(30\%\). However, the value of \(\sigma _8^0\) in this theory does not change much compared to \(\Lambda \)CDM model. We have also shown that the Einstein dark energy model can be compatible with observational data for redshifts greater than 0.2.

In summary we would like to say that the Einstein’s idea of considering elementary particles as electromagnetic like fields can be put forward to make some progress in obtaining satisfactory model of dark energy. But more data would still be needed to decide which model is more friendly with observations.