1 Introduction

One of the most interesting findings of modern cosmology is that our Universe is in a phase of accelerated expansion [1,2,3,4]. To explain this fundamental observation, one has to adopt a theory of gravity that admits an accelerating solution. Historically, such a solution was first obtained by the addition of the cosmological constant to the Einstein gravitational field equations. But we now know that gravitational models assuming the existence of the cosmological constant suffer from some theoretical problems, such as, for example, the fine tuning problem [5]. An alternative way to explain the recent acceleration of the Universe is to modify the dynamics of general relativity, such that the self-accelerated expansion is obtained from the modified Einstein gravitational field equations. There are two ways to achieve this goal. The first one is to generalize the Einstein–Hilbert action, such as in f(R) theories [6, 7], or, more generally, in f(R, Q) theories where Q can be any scalar combination of the curvature tensor/energy-momentum tensor [8,9,10,11,12,13,14]. However, these theories are generically unstable [15, 16]. The second possibility is enriching the graviton itself, such that it admits more than two degrees of freedom [17,18,19,20], and constructing the massive gravity theory. In the latter case, the helicity-0 mode of the massive graviton acts as a repulsive gravity and produce an accelerated expanding solution.

An alternative way to produce an accelerated expanding Universe is to add some light degrees of freedom to the theory of general relativity. This can be considered as an addition of some matter degrees of freedom to the theory. The simplest possibility is to add a scalar field such that its equation of state mimics the equation of state of the cosmological constant, e.g. \(p+\rho \approx 0\). This can be achieved by adding some higher powers of first derivative kinetic terms [21, 22], or by adding a potential term of the scalar field to the theory [25, 26]. Much work has been done in the context of scalar field cosmology, including inflationary [27,28,29] and dark energy models [30,31,32]. Recently, an interesting scalar field theory was proposed, which has higher than second order terms in the time derivative interactions in the action, although their equations of motion remains at most second order [33,34,35]. The cosmological implications of the Horndeski type theories have been extensively investigated in the literature [36,37,38,39,40,41,42,43].

Another interesting possibility is to add a vector degree of freedom to the theory of Einstein general relativity. The simplest example is to add a Maxwell kinetic term to the Einstein–Hilbert action and obtain an Einstein–Maxwell system [44]. One can also consider a massive vector field as an Einstein–Proca system. Much work has been done in the context of Einstein–Maxwell and Einstein–Proca theories [45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64]. Another possibility is to consider Weyl gravity, in which the metric compatibility condition does not hold anymore. As a result the covariant derivative of the metric tensor can be characterized by a vector field as \(\nabla _\mu g_{\nu \rho }=w_\mu g_{\nu \rho }\) [65,66,67,68,69,70,71,72,73,74,75,76]. One can also obtain the Einstein–Proca system from the higher dimensional Gauss–Bonnet theory in the context of Weyl geometry [78]. Cosmological implications of such a theory has been considered in [79]. One can also think of a generalization of the Proca action (as in the Horndeski theory), in which the action has higher order derivatives, with at most second order time derivative terms in the equations of motion similar to the case of scalar galileon. But assuming the U(1) symmetry, one can prove that such a vector galileon theory does not exist [80], meaning that the U(1)-vector galileon theory cannot exist. However, relaxing this constraint, one can generalize the Proca action, in a way that the zero helicity of the vector field mimics the scalar galileon interactions [81, 82]. The cosmology of this vector galileon theory, as well as a possible generalization, has been investigated in the literature [83,84,85,86,87,88].

Other higher derivative interaction terms which are not included in the vector galileon theory may also be considered. These terms in general produce Ostrogradski instabilities that cause the propagation of ghost degrees of freedom at large scales. However, one can tune the coupling constants in a way that the ghost becomes non-dynamical at scales comparable to the Hubble radius. This will make the theory effectively reliable in these scales [89].

Other possibility of generalizing Einstein’s general relativity within a purely geometric approach is to consider the effects of the torsion tensor [90,91,92,93,94,95]. The torsion tensor can in general be decomposed into a vector field \(Q_{\mu }\), an axial vector field \(S^{\mu }\) and a tensor field \(t_{\mu \nu \rho }\) satisfying the conditions \(t_{\mu \nu \rho }+t_{\nu \rho \mu }+t_{\rho \mu \nu }=0\), and \(g_{\mu \nu }t^{\mu \nu \rho }=0=g_{\mu \rho }t^{\mu \nu \rho }\), respectively. Assuming that the tensor \(t_{\mu \nu \rho }\) is zero, one can obtain a vector-tensor theory of gravity, which is well known in the context of supergravity theory. The cosmological implications of the Gauss–Bonnet theory coupled to a Weyl vector in Cartan space-time were considered in [96].

An interesting vector-tensor theory of gravity is the Einstein-aether theory [97, 98]. In this theory, one imposes a dynamical condition through a Lagrange multiplier in order to make the vector field always time-like. In this case the Universe has a predefined preferred time direction, which breaks the Lorentz invariance. Such a theory can explain the self-accelerated expansion of the Universe. The relation between the Einstein-aether theory and Horava–Lifshitz gravity [99, 100], as well as its relation to the scalar-tensor theory [101,102,103], has been carefully investigated in the physical literature.

Scalar field dark energy models have provided a very successful description of the observational properties of the Universe, including the explanation of its late acceleration. However, a priori one cannot simply reject the interesting idea that dark energy may have a more complex field structure. One such possibility is to describe dark energy in terms of a vector or Yang–Mills type field, which is also allowed to directly couple to gravity. The simplest possible action for a Yang–Mills type dark energy model minimally coupled to gravity is [104]

$$\begin{aligned} S_{V}= & {} -\int \mathrm{d}^{4}x\sqrt{-g}\Bigg \{\frac{R}{2}+\sum _{a=1}^{3}\Bigg [\frac{1 }{16\pi }F_{\mu \nu }^{a}F^{a}{}^{\mu \nu } +V(A^{2})\Bigg ] \nonumber \\&+L_{m}\Bigg \}, \end{aligned}$$
(1)

where \(F_{\mu \nu }^{a}=\nabla _{\mu }A_{\nu }^{a}-\nabla _{\nu }A_{\mu }^{a} \), \(a=1,2,3\), \(\nabla _{\mu }\) represents the covariant derivative with respect to the metric, while \(V(A^2)\) with \(A^{2}=g^{\mu \nu }A_{\mu }^{a}A_{\nu }^{a}\) represents a self-interaction potential explicitly violating gauge invariance. In the action given by Eq. (1), the dark energy component is represented by three vector fields, and thus the action (1) generalizes the Einstein–Maxwell single vector field dark energy model. The astrophysical and cosmological implications of the vector type dark energy models and of their generalizations have been extensively investigated in [105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,120,121,122].

Extended vector field dark energy models either, in which the vector field non-minimally couples to the gravitational field, have also been studied. Such models have been proposed, and their cosmological properties have been investigated in, for example, [123]. The action for the non-minimally massive vector field coupled to gravity can be represented as

$$\begin{aligned} S= & {} -\int \mathrm{d}^{4}x\sqrt{-g}\Bigg [\frac{R}{2}+\frac{1}{16\pi }F_{\mu \nu }F^{\mu \nu }-\frac{1}{2}\mu _{\Lambda }^{2}A_{\mu }A^{\mu } \nonumber \\&+\,\omega A_{\mu }A^{\mu }R+\eta A^{\mu }A^{\nu }R_{\mu \nu }+L_{m}\Bigg ], \end{aligned}$$
(2)

where \(\mu _{\Lambda }\) is the mass of the massive cosmological vector field, and \(A^{\mu }\left( x^{\nu }\right) \), \(\mu ,\nu =0,1,2,3\) is its four-potential, which is allowed to couple non-minimally to gravity. Here \(\omega \) and \(\eta \) are dimensionless coupling parameters. By following a close analogy with electrodynamics, the dark energy vector type field tensor is defined again by \(F_{\mu \nu }=\nabla _{\mu }A_{\nu }-\nabla _{\nu }A_{\mu }\).

A so called superconducting type dark energy model was recently introduced in [124, 125]. Inspired by some condensed matter concepts, the starting point of this approach is represented by the deep connection of the gravitational actions for scalar field models, given by \(S_{\phi }=-\int {\mathrm{d}^{4}x\sqrt{-g}\left[ \frac{R}{2}-\frac{1}{2}\nabla ^{\alpha }\phi \nabla _{\alpha }\phi +V(\phi )\right] }\), where \(V(\phi )\) the self-interaction potential of the scalar field, and (1), respectively. Despite having very different mathematical forms, the two actions for scalar and vector fields, respectively, can be interpreted and described as the two limiting cases of a single unified fundamental physical model that describes the spontaneous breaking of the U(1) symmetry of the “electromagnetic” type dark energy, with the corresponding action given by

$$\begin{aligned} S= & {} -\int \Bigg [ \frac{R}{2}+\frac{1}{16\pi }F_{\mu \nu }F^{\mu \nu }-\frac{ \lambda }{2} g^{\mu \nu } \left( A_{\mu }-\nabla _{\mu }\phi \right) \nonumber \\&\times \left( A_{\nu }-\nabla _{\nu }\phi \right) +V\left( A^2,\phi \right) - \frac{\alpha }{2} g^{\mu \nu }j_{\mu }\left( A_{\nu }-\nabla _{\nu }\phi \right) \nonumber \\&+L_{m}\Bigg ] \sqrt{-g}\mathrm{d}^4x , \end{aligned}$$
(3)

where \(\lambda \) and \(\alpha \) are arbitrary constants, \(L_m=L_{m}\left( g_{\mu \nu },\psi \right) \) is the Lagrangian of the total (ordinary baryonic plus dark) matter, and \(j^{\mu }=\rho u^{\mu }\) is the total mass current. In Eq. (3) \(\rho \) denotes the total matter density (including the dark matter one), while \(u^{\mu }\) is the matter four-velocity. Hence, as one can easily see, the gravitational action defined by Eq. (3) provides a unified theoretical framework for the scalar–vector interactions in a gravitational background.

In the theory of electromagnetism, one can add to the electromagnetic Lagrangian a very interesting interaction term, known as the “Chern–Simons” term, and defined as [126, 127]

$$\begin{aligned} k_\mu \epsilon ^{\mu \nu \rho \sigma }B_\nu B_{\rho \sigma }, \end{aligned}$$
(4)

where \(B_{\mu \nu }=\partial _\mu B_\nu -\partial _\nu B_\mu \), and \(k_\mu \) is a constant vector acting as a coupling constant. This term was initially defined as a topological mass term for gauge fields in (2+1) dimensions [128, 129]. One can see that because of the appearance of the Levi-Civita tensor, the CPT symmetry will be broken. On the other hand, because of a constant vector \(k_\mu \), the Poincaré invariance will be broken by this term. For example, for the time-like vector \(k_\mu \), the time translational symmetry together with a boost will be broken, while the spatial translation and rotations are preserved. For the space-like vector, the spatial rotations in a plane orthogonal to the vector is preserved. Many works has been done in the context of Chern–Simons electrodynamics, including applications to quantum electrodynamics [130,131,132]. In [133] the authors obtain the Chern–Simons term from dimensional reduction of the Carroll–Field–Jackiw Higgs model. Also in [134] the Chern–Simons theory with boundaries has been considered in more details.

The Chern–Simons modified gravity represents an interesting modification of general relativity, in which the Einstein–Hilbert action is extended by adding a parity-violating Chern–Simons term [135,136,137]. This term couples to gravity via a scalar field. The Chern–Simons correction enhances parity violation through a pure curvature term, as opposed to the matter term, as is considered in standard general relativity. It is important to note that Chern–Simons modified gravity can be obtained explicitly from superstring theory, where upon four-dimensional compactification the Chern–Simons term, appearing the Lagrangian density, plays an essential role due to the Green–Schwarz anomaly-canceling mechanism [138, 139]. Two distinct formulations of Chern–Simons modified gravity have been proposed, namely, the non-dynamical formulation and the dynamical formulation (see [137] for a review of the early results). In the first formulation, the Chern–Simons field is an arbitrary function, prescribed a priori, with its effective evolution equation being equivalent to a differential constraint on the space of its allowed solutions. In the second approach, the Chern–Simons field is treated as a dynamical field, with its own effective stress-energy tensor, and obeying a dynamical evolution equation. The possibility of observationally testing the dynamical Chern–Simons modified gravity by using the accretion disk properties around slowly rotating black holes was considered in [140]. Specific signatures do appear in the electromagnetic spectrum, thus leading to the possibility of directly testing Chern–Simons modified gravity by using astrophysical observations of the emission spectra from accretion disks.

In this paper, we are going to investigate the cosmological implications of the Chern–Simons term as a representative of dark energy effect. To this end, we will restore the Lorentz invariance of the theory by promoting the constant vector \(k_\mu \) to a dynamical field and impose a constant norm constraint by adding this property by a Lagrange multiplier. The addition of such a constraint by Lagrange multiplier to the action will spontaneously breaks the Lorentz symmetry of the theory. The Chern–Simons term will be generalized to

$$\begin{aligned} \alpha \epsilon ^{\mu \nu \rho \sigma }B_\mu A_\nu B_{\rho \sigma }+\lambda (A^\mu A_\mu +\eta ), \end{aligned}$$

where \(\lambda \) is a Lagrange multiplier, \(\eta \) and \(\alpha \) are constant and \(\epsilon ^{\mu \nu \rho \sigma }\) is the Levi-Civita tensor. The new vector field \(A_\mu \) plays the role of the constant vector \(k_\mu \) in the theory and it is constrained to have a constant norm. In the Minkowski background, if one wants to quantize the resulting electrodynamics theory, one should assume that the constant \(k_\mu \) be purely spatial and has a form \(k_\mu =(0,\vec {k})\) [141, 142]. In this paper, we will not be interested in the quantization of the theory and as a result we will not impose any constraint on the sign of the norm of the vector field \(A_\mu \). Also, as we will see in the following, in the case of FRW cosmology, one should assume that the vector field be time-like and has only a (t)-component, because otherwise the spatial isotropy will be broken. However, for the sake of completeness, we will consider the case of anisotropic cosmological implications of the theory and assume that the vector field \(A_\mu \) be space-like in this case.

The present paper is organized as follows. In the next section we will introduce the basic theoretical model, and obtain the necessary results for considering the cosmological implications of the model. In Sect. 3 we will consider the isotropic cosmology of the Chern–Simons model by obtaining the de Sitter solution, and performing the cosmological perturbations around it. In Sect. 4 we will assume that the vector \(A_\mu \) becomes space-like, and then consider the case of anisotropic cosmology with the geometry described by the Bianchi type I metric. The field equations and the basic physical parameters of the model are also introduced. In Sect. 5 we present a detailed numerical analysis of the evolution equations for different equations of state of the matter content of the Universe. In the last section we discuss our results and conclude.

2 The Maxwell–Chern–Simons gravity model

We propose an action functional of the form

$$\begin{aligned} S&=\int \mathrm{d}^4x \sqrt{-g}\bigg [\kappa ^2 R-\frac{1}{4 } A_{\mu \nu }A^{\mu \nu }- \frac{1}{4}B_{\mu \nu }B^{\mu \nu } \nonumber \\&\quad +\alpha \epsilon ^{\mu \nu \rho \sigma }B_\mu A_\nu B_{\rho \sigma }+\lambda (A^\mu A_\mu +\eta ) +V(A^2)\bigg ]+S_m, \end{aligned}$$
(5)

where \(A_\mu \) and \(B_\mu \) are two vectors fields, \(A^2=A_\mu A^\mu \), \(S_m\) is the matter action, \(\eta \) is the norm of vector field \(A_\mu \) with dimension of mass squared and \(\lambda \) is the Lagrange multiplier which dynamically enforces that the vector field \(A_\mu \) has a constant norm. This action consists of an Einstein–Maxwell system coupled to a constant norm vector field \(A_\mu \). The two vector fields then interact via the topological Chern–Simons term.

Varying the action with respect to the Lagrange multiplier results in

$$\begin{aligned} A_\mu A^\mu +\eta =0, \end{aligned}$$
(6)

which states that the vector field \(A_\mu \) should have a constant norm. One should note that because the norm of the vector field \(A_\mu \) is constant due to the constraint (6), the potential term \(V(A^2)\) in (5) is effectively the cosmological constant. However, certain choices of this potential may affect the dynamics of the model. Therefore, we will keep this term in the theory.

The equation of motion for \(B_\mu \) can be written as

$$\begin{aligned} \nabla _\beta B^{\beta \alpha }-\alpha \epsilon ^{\alpha \beta \gamma \delta }B_\beta A_{\gamma \delta }+2\alpha \epsilon ^{\alpha \beta \gamma \delta }A_\beta B_{\gamma \delta }=0. \end{aligned}$$
(7)

The equation of motion for \(A_\mu \) is

$$\begin{aligned} \nabla _\beta A^{\beta \alpha }+2\lambda A^\alpha +2V^\prime A^\alpha =\alpha \epsilon ^{\alpha \beta \gamma \delta }B_\beta B_{\gamma \delta }, \end{aligned}$$
(8)

where we have defined \(V^\prime =\mathrm{d}V/\mathrm{d}A^2\). One can see that the \(A_\mu \) field is generated by the \(B_\mu \) vector field through the Chern–Simons term. Also the \(B_\mu \) vector field is generated by the \(A_\mu \) field through the Chern–Simons term. In the case of vanishing \(\alpha \), the two vector fields evolve independently.

The equation of motion for the metric is

$$\begin{aligned}&\kappa ^2 G_{\alpha \beta }+\frac{1}{8 } g_{\alpha \beta }A_{\mu \nu }A^{\mu \nu }- \frac{1}{2}A_{\mu \alpha }A^\mu _{~\beta }+\frac{1}{8 } g_{\alpha \beta }B_{\mu \nu }B^{\mu \nu }\nonumber \\&\quad -\frac{1}{2}B_{\mu \alpha }B^\mu _{~\beta }+\lambda A_\alpha A_\beta \nonumber \\&\quad -\frac{1}{2}\lambda g_{\alpha \beta }(A_\mu A^\mu +\eta )-\frac{1}{2} V g_{\alpha \beta }\nonumber \\&\quad +V^\prime A_\alpha A_\beta =\frac{1}{2}T_{\alpha \beta }. \end{aligned}$$
(9)

In order to obtain the conservation of energy-momentum tensor in this model, first note that from the constraint equation (6), one has \(A^\mu \nabla _\alpha A_\mu =0\). Also, by taking the divergence of equation (8), one can obtain

$$\begin{aligned}&A^\alpha \nabla _\alpha \lambda +\lambda \nabla _\alpha A^\alpha +V^\prime \nabla _\alpha A^\alpha =\frac{1}{4}\epsilon ^{\alpha \beta \gamma \delta }B_{\alpha \beta }B_{\gamma \delta }. \end{aligned}$$
(10)

Now, taking the derivative of Eq. (9), we arrive at

$$\begin{aligned} \nabla ^\alpha T_{\alpha \beta }&=\frac{1}{4}\alpha \epsilon ^{\mu \rho \gamma \delta }\Big [2B_\rho B_{\gamma \delta }A_{\mu \beta }+B_{\mu \rho }B_{\gamma \delta }A_\beta \nonumber \\&\quad +2B_\rho A_{\gamma \delta }B_{\mu \beta }-4A_\rho B_{\gamma \delta }B_{\mu \beta }\Big ]. \end{aligned}$$
(11)

From the action functional (5), one can see that the ordinary matter sector is minimally coupled to the vector fields and also to gravity. From the diffeomorphism invariance of the action, since there is no energy transfer from matter to the vector fields, the matter sector should be conserved independently [143]. Hence, we obtain

$$\begin{aligned} \nabla ^\alpha T_{\alpha \beta }=0, \end{aligned}$$
(12)

and then

$$\begin{aligned}&\frac{1}{4}\alpha \epsilon ^{\mu \rho \gamma \delta }\Big [2B_\rho B_{\gamma \delta }A_{\mu \beta }+B_{\mu \rho }B_{\gamma \delta }A_\beta \nonumber \\&\quad +2B_\rho A_{\gamma \delta }B_{\mu \beta }-4A_\rho B_{\gamma \delta }B_{\mu \beta }\Big ]=0. \end{aligned}$$
(13)

One can see from the above expression that the vector fields can transfer energy to each other, which is due to the non-minimal Chern–Simons coupling between them. One should note that Eq. (11) is obtained by taking the covariant divergence of the metric field equation, and should be satisfied automatically for the given system of field equations (for more details see [144]).

3 Isotropic Cosmology of the Maxwell–Chern–Simons gravity

In this paper, we will consider the cosmology of Maxwell–Chern–Simons gravity theory. We will assume that the potential term in the action (5) is

$$\begin{aligned} V(A^2)=-2\kappa ^2\Lambda +\beta A_\mu A^\mu A_\nu A^\nu , \end{aligned}$$
(14)

where \(\Lambda \) is the cosmological constant.

3.1 de Sitter solution

Let us assume that our Universe can be described by a flat FRW metric of the form

$$\begin{aligned} \mathrm{d}s^2=-\mathrm{d}t^2+a^2(t)\big (\mathrm{d}x^2+\mathrm{d}y^2+\mathrm{d}z^2\big ). \end{aligned}$$
(15)

Moreover, because we are interested in the self-accelerated solutions to the theory, we assume that the energy-momentum tensor vanishes. In this case, the most general ansatz for the vector field which respects the spatial isotropy is

$$\begin{aligned} A_\mu&=(A_0,0,0,0), \end{aligned}$$
(16)
$$\begin{aligned} B_\mu&=(B_0,0,0,0). \end{aligned}$$
(17)

From the above equation, one can see that the vector field \(A_\mu \) should be time-like in this case, and we will assume that \(\eta =m^2>0\). For the space-like vector field \(\eta <0\), there is no FRW solution unless \(A_\mu =0\), which contradicts Eq. (6).

By substituting the above ansatz to the field equations (6)–(9), one can see that the constraint equation (6) gives \(A_0=m\). Also, the \( B_\mu \) equation of motion (7) is automatically satisfied in this case. The equations of motion for the vector field \(A_\mu \) and for the metric can then be written as

$$\begin{aligned}&2 m(\lambda -2\beta m^2)=0, \end{aligned}$$
(18)
$$\begin{aligned}&-3\kappa ^2 H^2+\kappa ^2\Lambda +\frac{3}{2}\beta m^4-m^2\lambda =0, \end{aligned}$$
(19)
$$\begin{aligned}&-2\kappa ^2\dot{H}-3\kappa ^2 H^2 +\kappa ^2 \Lambda -\frac{1}{2}\beta m^4=0, \end{aligned}$$
(20)

which can be solved as

$$\begin{aligned} \lambda =2\beta m^2,\quad H^2=\frac{\Lambda }{3}-\frac{\beta m^4}{6\kappa ^2}. \end{aligned}$$
(21)

We will also assume that \(B_0=\mathrm{constant}\) for the sake of simplicity. In order to have a consistent de Sitter solution, one should has \(\beta <2\kappa ^2 \Lambda /m^4\). Note that, in the absence of cosmological constant \(\Lambda \), the coupling constant \(\beta \) should be negative in order to have a de Sitter solution.

3.2 Cosmological perturbations

In this section we will perform the cosmological perturbation analysis on top of the de Sitter solution obtained in the previous section. For the metric perturbation, we have

$$\begin{aligned} \mathrm{d}s^2&=-(1+2\phi )\mathrm{d}t^2+2a(S_i+\partial _i B)\mathrm{d}t\mathrm{d}x^i\nonumber \\&\quad +a^2\big ((1+2\psi )\delta _{ij}+\partial _i\partial _j E+\partial _{(i}F_{j)}+h_{ij} \big )\mathrm{d}x^i \mathrm{d}x^j, \end{aligned}$$
(22)

where \(\phi \), \(\psi \), E and B are the scalar perturbations, \(F_i\) and \(S_i\) are the vector perturbations with the property \(\partial _i F_i=0=\partial _i S_i\), and \(h_{ij}\) is associated with the tensor perturbation, which is transverse and traceless \(h_{ii}=0=\partial _i h_{ij}\). Also, the spatial indices are raised and lowered by \(\delta _{ij}\).

For the vector fields, we decompose the perturbed vector field as

$$\begin{aligned} A_\mu&=(A_0+\delta A_0,\xi _i+\partial _i \delta A), \end{aligned}$$
(23)
$$\begin{aligned} B_\mu&=(B_0+\delta B_0,\epsilon _i+\partial _i \delta B), \end{aligned}$$
(24)

where \(\delta A_0\), \(\delta B_0\), \(\delta A\) and \(\delta B\) are the scalar perturbations, and \(\xi _i\) and \(\epsilon _i\) are the vector perturbations, respectively, with the property \(\partial _i\xi _i=0=\partial _i\epsilon _i\). The perturbation of the Lagrange multiplier is \(\lambda =\lambda _0+\delta \lambda \).

Under the general linear coordinate transformation \(x^\mu \rightarrow x^\mu +\delta x^\mu \), the metric perturbations transforms as

$$\begin{aligned}&\phi \rightarrow \phi -\partial _t\delta x^0,\end{aligned}$$
(25a)
$$\begin{aligned}&B\rightarrow B+\frac{1}{a}\delta x^0-a\partial _t\delta x,\end{aligned}$$
(25b)
$$\begin{aligned}&\psi \rightarrow \psi -H\delta x^0,\end{aligned}$$
(25c)
$$\begin{aligned}&E\rightarrow E-2\delta x,\end{aligned}$$
(25d)
$$\begin{aligned}&S_i\rightarrow S_i-a\partial _t\eta _i,\end{aligned}$$
(25e)
$$\begin{aligned}&F_i\rightarrow F_i-2\eta _i,\end{aligned}$$
(25f)
$$\begin{aligned}&h_{ij}\rightarrow h_{ij}, \end{aligned}$$
(25g)

where we have decomposed \(\delta x^\mu \) as \((\delta x^0,\eta _i+\partial _i\delta x)\), with \(\partial _i\eta _i=0\). Also the various components of a vector field \(A_\mu \) will be transformed according to

$$\begin{aligned}&\delta A_0\rightarrow \delta A_0-A_0\partial _t\delta x^0,\end{aligned}$$
(26a)
$$\begin{aligned}&\delta A\rightarrow \delta A-A_0\delta x^0,\end{aligned}$$
(26b)
$$\begin{aligned}&\xi _i\rightarrow \xi _i. \end{aligned}$$
(26c)

Likewise, we have a similar transformation rule for the vector field \(B^\mu \). Also, one can find that the Lagrange multiplier will not change under this transformation, and hence

$$\begin{aligned} \delta \lambda \rightarrow \delta \lambda . \end{aligned}$$
(27)

We should note that in obtaining the above relations we have assumed that the background values \(A_0\), \(B_0\) and \(\lambda _0\) are constant, as follows from the results of the previous section.

With the help of the above transformations, one can construct seven independent gauge invariant scalar perturbations, three independent gauge invariant vector perturbations, and one gauge invariant tensor perturbation, as follows:

$$\begin{aligned} \Phi&=\phi +\partial _t\left( aB-\frac{a^2}{2}\partial _t E\right) ,\nonumber \\ \Psi&=\psi +H\left( aB-\frac{a^2}{2}\partial _t E\right) ,\end{aligned}$$
(28)
$$\begin{aligned} \delta \mathcal {A}_0&=\delta A_0+A_0\partial _t\left( aB-\frac{a^2}{2}\partial _t E\right) ,\nonumber \\ \delta \mathcal {A}&=\delta A+A_0\left( aB-\frac{a^2}{2}\partial _t E\right) ,\end{aligned}$$
(29)
$$\begin{aligned} \delta \mathcal {B}_0&=\delta B_0+B_0\partial _t\left( aB-\frac{a^2}{2}\partial _t E\right) ,\nonumber \\ \delta \mathcal {B}&=\delta B+B_0\left( aB-\frac{a^2}{2}\partial _t E\right) . \end{aligned}$$
(30)

The perturbation of the Lagrange multiplier is already gauge invariant. The vector perturbations \(\xi _i\) and \(\epsilon _i\) are also gauge invariant, and the remaining gauge invariant vector perturbation can be constructed as

$$\begin{aligned} \beta _i=S_i-\frac{1}{2}a\partial _t F_i. \end{aligned}$$
(31)

Moreover, the tensor perturbation \(h_{ij}\) is already gauge invariant, and we end up with 11 gauge invariant perturbation variables.

After substituting the above expressions in (5), and expanding the resulting action up to second order in perturbations, one can see that the tensor, vector and scalar parts are completely decoupled from each other. In the following we will consider them separately.

3.2.1 Tensor perturbation

The tensor perturbation \(h_{ij}\) is transverse and traceless, and can be described by two polarization modes \(h_{+}\) and \(h_\times \). After Fourier decomposition, one can obtain the second order action of tensor perturbation as

$$\begin{aligned} S^{(2)}_{\mathrm{tensor}}=\frac{1}{2 } \kappa ^2\int \mathrm{d}t\mathrm{d}^3\vec {k} a^3\sum _{\lambda =+, \times }\left( |\dot{h}_\lambda |^2-\frac{\vec {k}^2}{a^2}|h_\lambda |^2\right) , \end{aligned}$$
(32)

where \(\vec {k}\) is the comoving wave vector, and the dot represents the time derivative. Also, it is understood that due to the definition of the Dirac delta function, in any term in the second order action, the argument of one of the perturbation variables is k, while for the other variable is \(-k\). One can see that the Chern–Simons coupling term does not contribute to the tensor perturbations. The action (32) is equivalent to that of Einstein–Hilbert theory. This conclusion can in fact be obtained without performing the explicit calculations, since the gravity sector of the model is the same as the one for the Einstein–Hilbert theory, and the vector fields are minimally coupled to gravity. Hence the theory has two tensor propagating modes associated to the massless graviton.

3.2.2 Vector perturbation

For the vector sector of the theory, we have three transverse independent gauge invariant vector perturbations \(\beta _i\), \(\xi _i\) and \(\epsilon _i\). After expanding the action up to the second order in the perturbed quantities and performing the Fourier transform, one can obtain

$$\begin{aligned} S^{(2)}_{\mathrm{vector}}&=\frac{1}{2 } \int \mathrm{d}t\mathrm{d}^3\vec {k} a \sum _{i=1}^2\left[ \kappa ^2 \vec {k}^2 \beta _i^2 +\dot{\xi }_i^2-\frac{\vec {k}^2}{a^2} \xi _i^2 +\dot{\epsilon }_i^2\right. \nonumber \\&\quad -\frac{\vec {k}^2}{a^2}\epsilon _i^2+4\frac{i\alpha B_0}{a^4}\vec {k}.\big (\vec {\xi }(-k)\times \vec {\epsilon }(k)\big ) \nonumber \\&\quad \left. -4\frac{i\alpha A_0}{a^4}\vec {k}.\big (\vec {\epsilon }(-k)\times \vec {\epsilon }(k)\big )\right] , \end{aligned}$$
(33)

where we have added the arguments of the last two terms for clarity. The vector perturbation of the metric \(\beta _i\) does not couple to the other vector fields, and is non-dynamical. One should note that the above expression for vector perturbation is real. For real perturbed quantities \(\epsilon _i(x)\) and \(\xi _i(x)\), one has \(\xi _i(k)=\xi ^\star _i(-k)\) and \(\epsilon _i(k)=\epsilon ^\star _i(-k)\). Now, calculating the complex conjugate of the last two terms in the action (33), one can easily show that these terms are real.

After obtaining the equation of motion of \(\beta _i\), and substituting it back to the action, one can see that \(\beta _i\) will vanish from the action. Therefore we are left with two vector perturbations associated to the two vector fields \(A_\mu \) and \(B_\mu \). Hence it follows that these vector perturbations interact non-trivially with each other through the Chern–Simons coupling.

Let us now consider the high curvature regime \(k\rightarrow \infty \) of the action (33). Since the vector perturbations \(\xi _i\) and \(\epsilon _i\) have the same order of magnitude, it turns out that the Chern–Simons terms can be neglected when compared with the terms containing \(\vec {k}^2\). In this case, one can write the action as

$$\begin{aligned} S^{(2)}_{\mathrm{vector},k\rightarrow \infty }=\frac{1}{2 } \int \mathrm{d}t\mathrm{d}^3\vec {k} a \sum _{i=1}^2\left[ \dot{\xi }_i^2-\frac{\vec {k}^2}{a^2} \xi _i^2 +\dot{\epsilon }_i^2-\frac{\vec {k}^2}{a^2}\epsilon _i^2\right] , \end{aligned}$$
(34)

which shows that the theory has two non-interacting vector modes. Hence, we remain with four healthy vector degrees of freedom for the theory.

3.2.3 Scalar perturbation

Finally, we will obtain the scalar perturbation part of the theory. There are five scalar modes, corresponding to the metric, the vector fields and the Lagrange multiplier. Expanding the action up to second order and performing the Fourier decomposition, one can write the action in terms of gauge invariant quantities as

$$\begin{aligned} S^{(2)}_{\mathrm{scalar}}&=\int \mathrm{d}t\mathrm{d}^3\vec {k}a^3\Big [-6\kappa ^2\dot{\Psi }^2+12\kappa ^2 H\Phi \dot{\Psi }\nonumber \\&\quad +2\kappa ^2\frac{\vec {k}^2}{a^2}\left( 2\Phi \Psi +\Psi ^2\right) +2 \left( 2\beta m^4-3\kappa ^2 H^2\right) \Phi ^2 \nonumber \\&\quad -8\beta m^3\Phi \delta \mathcal {A}_0+\frac{1}{2}\frac{\vec {k}^2}{a^2}\left( \dot{\delta \mathcal {A}}^2-2\dot{\delta \mathcal {A}}\delta \mathcal {A}_0+\delta \mathcal {A}_0^2\right) \nonumber \\&\quad +4\beta m^2 \delta \mathcal {A}_0^2+\frac{1}{2}\frac{\vec {k}^2}{a^2}\left( \dot{\delta \mathcal {B}}^2-2\dot{\delta \mathcal {B}}\delta \mathcal {B}_0+\delta \mathcal {B}_0^2\right) \nonumber \\&\quad -2m\delta \mathcal {A}_0\delta \lambda +2m^2\Phi \delta \lambda \Big ]. \end{aligned}$$
(35)

From the above equation one can see that \(\delta \lambda \), \(\delta \mathcal {B}_0\), \( \delta \mathcal {A}_0\), and \(\Phi \) are non-dynamical. Varying the action with respect to these variables leads to

$$\begin{aligned}&2m(m\Phi -\delta \mathcal {A}_0)=0, \end{aligned}$$
(36)
$$\begin{aligned}&\frac{\vec {k}^2}{a^2}(\delta \mathcal {B}_0-\dot{\delta }\mathcal {B})=0, \end{aligned}$$
(37)
$$\begin{aligned}&2m\delta \lambda -\left( \frac{\vec {k}^2}{a^2}+8\beta m^2\right) \delta \mathcal {A}_0+8\beta m^3\Phi +\frac{\vec {k}^2}{a^2}\dot{\delta \mathcal {A}}=0, \end{aligned}$$
(38)
$$\begin{aligned}&m^2\delta \lambda +2\kappa ^2\frac{\vec {k}^2}{a^2}\Psi +4\beta m^3(m\Phi -\delta \mathcal {A}_0)\nonumber \\&\quad +6\kappa ^2 H(\dot{\Psi }-H\Phi )=0. \end{aligned}$$
(39)

Solving the above equations for the non-dynamical variables, and substituting them back into the action gives

$$\begin{aligned} S^{(2)}_{\mathrm{scalar}}&=\int \mathrm{d}t\mathrm{d}^3\vec {k}\frac{2\vec {k}^2\kappa ^2 a}{H\left( 12\kappa ^2 H^2-m^2\frac{\vec {k}^2}{a^2}\right) } \nonumber \\&\quad \times \left[ H\left( 12\kappa ^2H^2-m^2\frac{ \vec {k}^2}{a^2}+4 \kappa ^2\frac{\vec {k}^2}{a^2}\right) \Psi ^2\right. \nonumber \\&\quad -2\left( m^2\frac{\vec {k}^2}{a^2}-12\kappa ^2 H^2\right) \Psi \dot{\Psi } \nonumber \\&\quad \left. -2m\frac{\vec {k}^2}{a^2}\Psi \dot{\Xi }+3H\dot{\Xi }^2\right] , \end{aligned}$$
(40)

where we have defined a new variable \(\Xi =H{\delta \mathcal {A}}-m{\Psi }\). Varying the action with respect to \(\Psi \) results in

$$\begin{aligned} \Psi =\frac{m}{4H\kappa ^2}\dot{\Xi }. \end{aligned}$$
(41)

Upon substituting the above relation into the action, one obtains

$$\begin{aligned} S^{(2)}_{\mathrm{scalar}}=\int \mathrm{d}t\mathrm{d}^3\vec {k}\frac{\vec {k}^2}{2H^2}a\dot{\Xi }^2, \end{aligned}$$
(42)

which can be solved for \(\Xi \) with the result

$$\begin{aligned} \Xi =c_1-\frac{c_2}{H} e^{-Ht}. \end{aligned}$$

Therefore, the scalar part of the perturbed action has one degree of freedom, which decays exponentially in time. The scalar degree of freedom does exist because the U(1) symmetry is broken by the addition of the Lagrange multiplier constraint. This shows that the vector field \(A_\mu \) is massive, while the vector field \(B_\mu \) is massless.

Hence, to summarize the results of the perturbative analysis of the present theory, we have found that it contains two gravitational wave modes, four vector modes, associated to the two vector fields, and one scalar degree of freedom, respectively.

4 Anisotropic cosmology: Bianchi type I Universe

In this section, we want to consider the case that the vector field \(A_\mu \) is space-like. So, from now on we will set \(\eta =-\mu ^2\). In this case, at least one of the spatial components of \(A_\mu \) should be non-zero. For simplicity, we choose the coordinate system such that the direction of the x axis coincides with the direction of \(A_\mu \). In this case the vector field \(A_\mu \) has the from

$$\begin{aligned} A_\mu =(A_0(t),A_1(t),0,0), \end{aligned}$$
(43)

In this case the vector field \(B_\mu \) may have any direction in space and so we set

$$\begin{aligned} B_\mu =(B_0(t),B_1(t),B_2(t),B_3(t)). \end{aligned}$$
(44)

The above ansatz for the vector field will break the spatial isotropy of the space-time. In order to consider the cosmology of this case, let us assume that the Universe is described by the Bianchi type I metric of the form

$$\begin{aligned} \mathrm{d}s^2=-\mathrm{d}t^2+\sum _{i=1}^{3}a_i^2(t)(\mathrm{d}x^{i})^{2}, \end{aligned}$$
(45)

where the \(a_i\) are the directional scale factors. We also assume that the matter content of the Universe has the form

$$\begin{aligned} T^{\mu }_\nu =\mathrm {diag}\big (-\rho (t),p_1(t),p_2(t),p_3(t)\big ), \end{aligned}$$
(46)

where \(\rho \) is the energy density and \(p_i\), \(i=1,2,3\), is the pressure in the direction i. Let us define the Hubble parameter, the deceleration parameter and the anisotropy parameter as [145,146,147]

$$\begin{aligned}&H_i=\frac{\dot{a_i}}{a_i},i=1,2,3,\qquad V=a_1a_2a_3, \quad H=\frac{1}{3} \sum _{i=1}^{3} H_i, \end{aligned}$$
(47)
$$\begin{aligned}&A=\frac{1}{3}\sum _{i=1}^{3}\left( \frac{\Delta H_i}{H}\right) ^2, \quad \Delta H_i=H-H_i, i=1,2,3,\nonumber \\ \end{aligned}$$
(48)
$$\begin{aligned}&\Sigma ^2=\frac{1}{2}\left( \sum _{i=1}^{3} H_i^2-3H^2\right) , \quad q=\frac{\mathrm{d}}{\mathrm{d}t}\left( \frac{1}{H}\right) -1. \end{aligned}$$
(49)

From Eqs. (47) we obtain immediately the important relation

$$\begin{aligned} \frac{\dot{V}}{V}=3H. \end{aligned}$$
(50)

The anisotropy parameter can be represented in an equivalent form as

$$\begin{aligned} A=\frac{1}{3H^2}\left( \sum _{i=1}^3{H_i^2}-3H^2\right) , \end{aligned}$$
(51)

giving

$$\begin{aligned} \sum _{i=1}^3{H_i^2}=3\left( 1+A\right) H^2. \end{aligned}$$
(52)

4.1 Gravitational field equations

Before writing the equations of motion of the theory, we should note that the (t)-component of the vector field \(B_{\mu }\) does not contribute to the field equations. From now on, we will assume that \(B_{0}=0\). On the other hand, the (t)-component of the \(A_{\mu }\) equation of motion (8), implies that \(A_{0}=0\). So, the constraint equation (6) gives the remaining component of the vector field \(A_{\mu }\) as \(A_{1}=\mu a_{1}(t).\) Also, the off-diagonal elements (ij), \(i\ne j\), of the metric field equation gives \(\dot{B}_{i}\dot{B}_{j}=0.\) This implies that at least two out of three components of the vector field \(B_{\mu }\) should be constant. We will assume that \(B_{1}\) and \(B_{2}\) is constant. Now, the (y) -component of the \(A_{\mu }\) field equation gives \(B_{1}\dot{B}_{3}=0\), which implies either \(B_{1}=0\) or \(B_{3}\) is constant. We will choose the first possibility to make the vector field \(B_{\mu }\) evolves in time. In this case, one can see that the vector fields \(A_\mu \) and \(B_\mu \) are orthogonal. With these in hand, the only remaining component of the \(A_{\mu }\) field equation becomes

$$\begin{aligned} 2\left( \lambda +2\beta \mu ^{2}\right) +H_{1}^{2}+2\frac{\alpha }{\mu } b_{2}(b_{3}H_{3}+\dot{b}_{3})-\frac{1}{V}\frac{\mathrm{d}}{\mathrm{d}t}\left( VH_{1}\right) =0. \end{aligned}$$
(53)

The remaining components of the \(B_{\mu }\) equation of motion are

$$\begin{aligned} 2\left( b_{3}H_{3}+\dot{b}_{3}\right) +b_{3}H_{1}=0, \end{aligned}$$
(54)

or, equivalently,

$$\begin{aligned} H_1=2H_3+2\frac{\dot{b_3}}{b_3}, \end{aligned}$$
(55)

and

$$\begin{aligned} -2\mu \alpha b_{2}H_{1}+b_{3}H_{3}^{2}-3H\dot{b}_{3}-\ddot{b}_{3}-\frac{b_{3} }{V}\frac{\mathrm{d}}{\mathrm{d}t}\left( VH_{3}\right) =0.\qquad \end{aligned}$$
(56)

The Friedmann equation can be written as

$$\begin{aligned}&\kappa ^{2}\left( 3\dot{H}+\sum _{i=1}^{3}H_{i}^{2}-\Lambda \right) +\frac{1}{ 4}\left( \frac{1}{4}b_{3}^{2}+\mu ^{2}\right) H_{1}^{2}\nonumber \\&\quad -\frac{1}{2}\mu ^{2}(\lambda +\beta \mu ^{2})=-\frac{1}{4}\left( \rho +\sum _{i=1}^{3}p_{i}\right) , \end{aligned}$$
(57)

and the Raychaudhuri equations are

$$\begin{aligned}&\kappa ^{2}\left( \frac{1}{V}\frac{\mathrm{d}}{\mathrm{d}t}\left( VH_{1}\right) -\Lambda \right) -\frac{1}{4}\left( \frac{1}{4}b_{3}^{2}-\mu ^{2}\right) H_{1}^{2}\nonumber \\&\quad + \frac{1}{2}\mu ^{2}(\lambda +3\beta \mu ^{2})=\frac{1}{4}(\rho +p_{1}-p_{2}-p_{3}), \end{aligned}$$
(58)
$$\begin{aligned}&\kappa ^{2}\left( \frac{1}{V}\frac{\mathrm{d}}{\mathrm{d}t}\left( VH_{2}\right) -\Lambda \right) -\frac{1}{4}\left( \frac{1}{4}b_{3}^{2}+\mu ^{2}\right) H_{1}^{2}\nonumber \\&\quad - \frac{1}{2}\mu ^{2}(\lambda +\beta \mu ^{2})=\frac{1}{4}(\rho -p_{1}+p_{2}-p_{3}), \end{aligned}$$
(59)

and

$$\begin{aligned}&\kappa ^{2}\left( \frac{1}{V}\frac{\mathrm{d}}{\mathrm{d}t}\left( VH_{3}\right) -\Lambda \right) +\frac{1}{4}\left( \frac{1}{4}b_{3}^{2}-\mu ^{2}\right) H_{1}^{2}\nonumber \\&\quad -\frac{1}{2}\mu ^{2}(\lambda +\beta \mu ^{2})=\frac{1}{4}(\rho -p_{1}-p_{2}+p_{3}). \end{aligned}$$
(60)

In the above equations we have defined

$$\begin{aligned} b_{2}(t)=\frac{B_{2}(t)}{a_{2}(t)}, \quad b_{3}(t)=\frac{B_{3}(t)}{a_{3}(t)}. \end{aligned}$$
(61)

Equation (55) can be solved for \(b_{3}\) with the result

$$\begin{aligned} b_{3}(t)=\frac{c_{3}}{a_{3}\sqrt{a_{1}}},\qquad c_{3}=\mathrm{const}. \end{aligned}$$
(62)

From Eq. (12), one has the conservation of the energy-momentum tensor

$$\begin{aligned} \dot{\rho }+3H\rho +\sum _{i=1}^{3}H_{i}p_{i}=0. \end{aligned}$$
(63)

Also, Eq. (13) satisfies automatically.

4.2 Deceleration parameter and anisotropy

In the following we will restrict our analysis to the case of a geometrically anisotropic Universe filled with a cosmological fluid with an isotropic pressure distribution, with \(p_{1}=p_{2}=p_{3}=p\).

By adding Eqs. (58)–(60) we obtain

$$\begin{aligned}&\kappa ^{2}\frac{1}{V}\frac{\mathrm{d}}{\mathrm{d}t}\left( 3HV\right) =3\kappa ^{2}\left( \dot{ H}+3H^{2}\right) =\kappa ^2 \frac{\ddot{V}}{V}=3\kappa ^{2}\Lambda \nonumber \\&\quad +\frac{1}{4 }\left( \frac{1}{4}b_{3}^{2}+\mu ^{2}\right) H_{1}^{2}+\frac{\mu ^{2}}{2} \left( \lambda -\beta \mu ^{2}\right) +\frac{3}{4}\left( \rho -p\right) .\nonumber \\ \end{aligned}$$
(64)

By substituting \(3\dot{H}\) from the above equation into Eq. (57) we find the consistency condition

$$\begin{aligned}&\kappa ^{2}\left( -9H^{2}+\sum _{i=1}^{3}H_{i}^{2}+2\Lambda \right) +\frac{1}{ 2}\left( \frac{1}{4}b_{3}^{2}+\mu ^{2}\right) \nonumber \\&\quad \times H_{1}^{2}-\beta \mu ^{4}=-\rho . \end{aligned}$$
(65)

With the use of Eqs. (54), (53) becomes

$$\begin{aligned} \frac{1}{V}\frac{\mathrm{d}}{\mathrm{d}t}\left( VH_{1}\right) =2\left( \lambda +2\beta \mu ^{2}\right) +H_{1}^{2}-\frac{\alpha }{\mu }b_{2}b_{3}H_{1}. \end{aligned}$$
(66)

By combining Eqs. (58) and (66) we find

$$\begin{aligned}&\left[ 1-\frac{1}{4\kappa ^{2}}\left( \frac{1}{4}b_{3}^{2}-\mu ^{2}\right) \right] H_{1}^{2}-\frac{\alpha }{\mu }b_{2}b_{3}H_{1} \nonumber \\&\quad -\Lambda +\lambda \left( 2+\frac{\mu ^{2}}{2\kappa ^{2}}\right) +\beta \mu ^{2}\left( 4+\frac{ 3\mu ^{2}}{2\kappa ^{2}}\right) =\frac{1}{4\kappa ^{2}}(\rho -p).\nonumber \\ \end{aligned}$$
(67)

With the help of Eqs. (55), (56) can be reformulated as

$$\begin{aligned} \frac{1}{V}\frac{\mathrm{d}}{\mathrm{d}t}\left( VH_{3}\right)= & {} \frac{1}{4}H_{1}^{2}-\frac{1}{ b_{3}}\left( \dot{b}_{3}+2\mu \alpha b_{2}\right) H_{1}\nonumber \\&-\frac{1}{V}\frac{\mathrm{d}}{ \mathrm{d}t}\left( V\frac{\dot{b}_{3}}{b_{3}}\right) , \end{aligned}$$
(68)

giving, after substitution into Eq. (60), the equation

$$\begin{aligned}&\frac{1}{4}\left[ 1+\frac{1}{\kappa ^{2}}\left( \frac{1}{4}b_{3}^{2}-\mu ^{2}\right) \right] H_{1}^{2}-\frac{1}{b_{3}}\left( \dot{b}_{3}+2\mu \alpha b_{2}\right) H_{1} \nonumber \\&\quad -\frac{1}{V}\frac{\mathrm{d}}{\mathrm{d}t}\left( V\frac{\dot{b}_{3}}{b_{3}} \right) -\Lambda -\frac{1}{2\kappa ^{2}}\mu ^{2}(\lambda +\beta \mu ^{2})= \frac{1}{4\kappa ^{2}}(\rho -p)\nonumber \\ \end{aligned}$$
(69)

Adding Eqs. (67) and (69) we obtain

$$\begin{aligned}&\frac{5}{4}H_{1}^{2}-\frac{1}{b_{3}}\left( \frac{\alpha }{\mu } b_{2}b_{3}^{2}+\dot{b}_{3}+2\mu \alpha b_{2}\right) H_{1}-\frac{1}{V}\frac{\mathrm{d} }{\mathrm{d}t}\left( V\frac{\dot{b}_{3}}{b_{3}}\right) \nonumber \\&\quad -2\Lambda +\left( 4\beta \mu ^{2}+2\lambda +\beta \frac{\mu ^{4}}{\kappa ^{2}}\right) =\frac{1 }{2\kappa ^{2}}(\rho -p), \end{aligned}$$
(70)

while equating Eqs. (67) and (69) gives

$$\begin{aligned}&\frac{1}{2}\left[ \frac{3}{2}-\frac{1}{\kappa ^2}\left( \frac{1}{4}b_{3}^{2}-\mu ^{2}\right) \right] H_{1}^{2} \nonumber \\&\quad -\frac{1}{b_{3}}\left[ \frac{\alpha }{\mu } b_{2}b_{3}^{2}-\left( \dot{b}_{3}+2\mu \alpha b_{2}\right) \right] H_{1}+\frac{1}{V}\frac{\mathrm{d}}{\mathrm{d}t}\left( V\frac{\dot{b}_{3}}{b_{3}}\right) \nonumber \\&\quad +\frac{1}{ \kappa ^{2}}\left( 2\beta \mu ^{2}+\lambda \right) \left( 2\kappa ^{2}+\mu ^{2}\right) =0. \end{aligned}$$
(71)

From Eq. (64) we obtain immediately the deceleration parameter as

$$\begin{aligned} q= & {} 2-\frac{1}{H^2}\left[ \Lambda +\frac{1}{12\kappa ^{2}}\left( \frac{1}{4} b_{3}^{2}+\mu ^{2}\right) H_{1}^{2}\right. \nonumber \\&\left. +\frac{\mu ^{2}}{6\kappa ^{2}}\left( \lambda -\beta \mu ^{2}\right) +\frac{1}{4\kappa ^{2}}\left( \rho -p\right) \right] , \end{aligned}$$
(72)

while from Eq. (65) we obtain the anisotropy parameter in the form

$$\begin{aligned} A= & {} 2+\beta \frac{\mu ^{4}}{3H^2\kappa ^{2}}-\frac{2}{3}\frac{\Lambda }{H^{2}}- \frac{1}{6\kappa ^{2}}\frac{1}{H^{2}}\left( \frac{1}{4}b_{3}^{2}+\mu ^{2}\right) H_{1}^{2} \nonumber \\&-\frac{1}{3H^2\kappa ^{2}}\rho . \end{aligned}$$
(73)

4.3 The evolution equations for the Bianchi type I cosmological model

Since the field \(b_{2}\) does not appear in the expressions of q and A, we can take \(b_{2}=0\) without any lack of generality. Then from Eq. (67) we obtain for \(H_{1}\) the expression

$$\begin{aligned} H_{1}^{2}=\frac{\frac{1}{4\kappa ^{2}}(\rho -p)+\Lambda _{\mathrm{eff}}}{1-\frac{1}{ 4\kappa ^{2}}\left( \frac{1}{4}b_{3}^{2}-\mu ^{2}\right) }, \end{aligned}$$
(74)

where we have denoted

$$\begin{aligned} \Lambda _{\mathrm{eff}}=\Lambda -\lambda \left( 2+\frac{\mu ^{2}}{2\kappa ^{2}} \right) -\beta \mu ^{2}\left( 4+\frac{3\mu ^{2}}{2\kappa ^{2}}\right) . \end{aligned}$$
(75)

By substituting this expression of \(H_{1}^{2}\) into Eq. (64) we obtain the evolution equation for V as

$$\begin{aligned} \frac{\ddot{V}}{V}= & {} 3\Lambda +\frac{\mu ^{2}}{2\kappa ^{2}}\left( \lambda -\beta \mu ^{2}\right) +\frac{1}{4\kappa ^{2}}\frac{\left( \frac{1}{4} b_{3}^{2}+\mu ^{2}\right) }{1-\frac{1}{4\kappa ^{2}}\left( \frac{1}{4} b_{3}^{2}-\mu ^{2}\right) } \nonumber \\&\times \left[ \frac{1}{4\kappa ^{2}}(\rho -p)+\Lambda _{\mathrm{eff}}\right] +\frac{3}{4\kappa ^{2}}\left( \rho -p\right) . \end{aligned}$$
(76)

Equation (69) gives the time evolution of the field \(b_{3}\) as

$$\begin{aligned} \frac{1}{V}\frac{\mathrm{d}}{\mathrm{d}t}\left( V\frac{\dot{b}_{3}}{b_{3}}\right)= & {} \frac{1}{4 }\left[ 1+\frac{1}{\kappa ^{2}}\left( \frac{1}{4}b_{3}^{2}-\mu ^{2}\right) \right] \nonumber \\&\times \frac{\frac{1}{4\kappa ^{2}}(\rho -p)+\Lambda _{\mathrm{eff}}}{1-\frac{1}{ 4\kappa ^{2}}\left( \frac{1}{4}b_{3}^{2}-\mu ^{2}\right) }\nonumber \\&-\frac{\dot{b}_{3} }{b_{3}}\sqrt{\frac{\frac{1}{4\kappa ^{2}}(\rho -p)+\Lambda _{\mathrm{eff}}}{1-\frac{1 }{4\kappa ^{2}}\left( \frac{1}{4}b_{3}^{2}-\mu ^{2}\right) }}\nonumber \\&-\Lambda -\frac{1}{2\kappa ^{2}}\mu ^{2}(\lambda +\beta \mu ^{2})-\frac{1}{ 4\kappa ^{2}}(\rho -p).\nonumber \\ \end{aligned}$$
(77)

4.4 Dimensionless form of the cosmological evolution equations

In order to simplify the mathematical formalism we rescale the physical and geometrical quantities by introducing the set of dimensionless variables \(\left( \tau ,r,P,B_{0},\mu _{0},\beta _{0},h\right) \), defined as

$$\begin{aligned} t&= \frac{1}{\sqrt{3\Lambda +\frac{\mu ^{2}}{2\kappa ^{2}}\left( \lambda -\beta \mu ^{2}\right) }}\tau ,\nonumber \\ H&=\sqrt{3\Lambda +\frac{\mu ^{2}}{ 2\kappa ^{2}}\left( \lambda -\beta \mu ^{2}\right) }\,h, \nonumber \\ \rho&= 4\kappa ^{2}\Lambda _{\mathrm{eff}}r,\quad p=4\kappa ^{2}\Lambda _{\mathrm{eff}}P,\quad b_{3}=2\kappa B_{0},\nonumber \\&\mu =\kappa \mu _{0}, \quad \beta =\frac{\beta _{0}}{\kappa ^{2}}. \end{aligned}$$
(78)

Then the system of equations (76) and (77) becomes

$$\begin{aligned} \frac{1}{V}\frac{\mathrm{d}^{2}V}{\mathrm{d}\tau ^{2}}= & {} 1+\lambda _{1}\frac{\left( B_{0}^{2}+\mu _{0}^{2}\right) }{1-\frac{1}{4}\left( B_{0}^{2}-\mu _{0}^{2}\right) }\left[ (r-P)+1\right] \nonumber \\&+12\lambda _{1}\left( r-P\right) , \end{aligned}$$
(79)
$$\begin{aligned}&\frac{1}{V}\frac{\mathrm{d}}{\mathrm{d}\tau }\left( V\frac{1}{B_{0}}\frac{\mathrm{d}B_{0}}{\mathrm{d}\tau } \right) \nonumber \\&\quad =\lambda _{1}\left[ 1+\left( B_{0}^{2}-\mu _{0}^{2}\right) \right] \frac{(r-P)+1}{1-\frac{1}{4}\left( B_{0}^{2}-\mu _{0}^{2}\right) }\nonumber \\&\qquad -2\sqrt{ \lambda _{1}}\frac{1}{B_{0}}\frac{\mathrm{d}B_{0}}{\mathrm{d}\tau }\sqrt{\frac{(r-P)+1}{1- \frac{1}{4}\left( B_{0}^{2}-\mu _{0}^{2}\right) }} \nonumber \\&\qquad -4\lambda _{1}(r-P)-\lambda _{2}, \end{aligned}$$
(80)

where

$$\begin{aligned} \lambda _{1}=\frac{\Lambda -\lambda \left( 2+\frac{\mu _{0}^{2}}{2}\right) -\beta _{0}\mu _{0}^{2}\left( 4+\frac{3\mu _{0}^{2}}{2}\right) }{4\left[ 3\Lambda +\frac{\mu _{0}^{2}}{2}\left( \lambda -\beta _{0}\mu _{0}^{2}\right) \right] } \end{aligned}$$
(81)

and

$$\begin{aligned} \lambda _{2}=\frac{\Lambda +\frac{1}{2}\mu _{0}^{2}(\lambda +\beta _{0}\mu _{0}^{2})}{3\Lambda +\frac{\mu _{0}^{2}}{2}\left( \lambda -\beta _{0}\mu _{0}^{2}\right) }. \end{aligned}$$
(82)

The energy conservation equation can be written as

$$\begin{aligned} \frac{\mathrm{d}r}{\mathrm{d}\tau }+3h\left( r+P\right) =0. \end{aligned}$$
(83)

By assuming an equation of state of the form

$$\begin{aligned} P=(\gamma -1)r, \quad \gamma =\mathrm {constant}, \quad 1\le \gamma \le 2, \end{aligned}$$
(84)

we obtain for the energy density of the cosmological matter the expression

$$\begin{aligned} r=\frac{r_{0}}{V^{\gamma }}, \end{aligned}$$
(85)

where \(r_{0}\) is an arbitrary integration constant. The numerical value of \(r_0\) can be taken as one without any loss of generality, so that the present density of the Universe, at the time \(t=t_{\mathrm{pres}}\), is given by \(\rho \left( t_{\mathrm{pres}}\right) =4\kappa ^2\Lambda _{\mathrm{eff}}/V^{\gamma }\left( t_{\mathrm{pres}}\right) \). Therefore the system of equations (79) and (80), describing the evolution of a Bianchi type I Universe in the Maxwell–Chern–Simons theory can be reformulated as a first order dynamical system given by

$$\begin{aligned} \frac{\mathrm{d}V}{\mathrm{d}\tau }= & {} u, \end{aligned}$$
(86)
$$\begin{aligned} \frac{\mathrm{d}u}{\mathrm{d}\tau }= & {} \left\{ 1+\lambda _{1}\frac{\left( B_{0}^{2}+\mu _{0}^{2}\right) }{1-\frac{1}{4}\left( B_{0}^{2}-\mu _{0}^{2}\right) }\left[ \frac{\left( 2-\gamma \right) r_{0}}{V^{\gamma }}+1\right] \right. \nonumber \\&\left. +12\lambda _{1} \frac{\left( 2-\gamma \right) r_{0}}{V^{\gamma }}\right\} V, \end{aligned}$$
(87)
$$\begin{aligned} \frac{\mathrm{d}B_{0}}{\mathrm{d}\tau }= & {} b_{0}, \end{aligned}$$
(88)
$$\begin{aligned} \frac{\mathrm{d}b_{0}}{\mathrm{d}\tau }= & {} \frac{1}{B_{0}}b_{0}^{2}+\lambda _{1}\left[ 1+\left( B_{0}^{2}-\mu _{0}^{2}\right) \right] \frac{\left( 2-\gamma \right) r_{0}V^{-\gamma }+1}{1-\frac{1}{4}\left( B_{0}^{2}-\mu _{0}^{2}\right) } B_{0}\nonumber \\&-\left[ 2\sqrt{\lambda _{1}}\sqrt{\frac{\left( 2-\gamma \right) r_{0}V^{-\gamma }+1}{1-\frac{1}{4}\left( B_{0}^{2}-\mu _{0}^{2}\right) }}+ \frac{u}{V}\right] b_{0} \nonumber \\&-\left[ 4\lambda _{1}\frac{\left( 2-\gamma \right) r_{0}}{V^{\gamma }} +\lambda _{2}\right] B_{0}. \end{aligned}$$
(89)

The system of differential equations (86)–(89) must be integrated with the initial conditions \(V(0)=V_{0}\), \(u(0)=u_{0}\), \( B_{0}(0)=B_{0}^{(0)}\), and \(b_{0}(0)=b_{0}^{(0)}\), respectively. The behavior of the solutions of the dynamical system is determined by two arbitrary parameters, \(\lambda _{1}\) and \(\lambda _{2}\), representing the combination of the four free model parameters \(\left( \Lambda ,\lambda ,\beta ,\mu \right) \).

The deceleration parameter can be written as

$$\begin{aligned} q= & {} 2-9\frac{V^{2}}{u^{2}}\left[ \frac{1}{3}+\frac{1}{3}\lambda _{1}\left( B_{0}^{2}+\mu _{0}^{2}\right) \frac{\left( 2-\gamma \right) r_{0}V^{-\gamma }+1}{1-\frac{1}{4}\left( B_{0}^{2}-\mu _{0}^{2}\right) }\right. \nonumber \\&\left. +\,4\lambda _{1}\frac{ \left( 2-\gamma \right) r_{0}}{V^{\gamma }}\right] . \end{aligned}$$
(90)

From Eq. (57) we obtain the following relation between the deceleration and the anisotropy parameters:

$$\begin{aligned} A= & {} q+\frac{3V^{2}}{u^{2}}\left\{ \lambda _2-\lambda _{1}\frac{\left( B_{0}^{2}+\mu _{0}^{2}\right) }{1-\frac{1}{4}\left( B_{0}^{2}-\mu _{0}^{2}\right) }\left[ \left( 2-\gamma \right) r_{0}V^{-\gamma }+1\right] \right. \nonumber \\&\left. -4\left( 3\gamma -2\right) \lambda _{1}\frac{r_{0}}{V^{\gamma }}\right\} . \end{aligned}$$
(91)

Hence the behavior of both q and A is determined by the set of model parameters \(\left( \mu _{0},\lambda _{1},\lambda _2,r_{0}\right) \).

5 Cosmological evolution of the Bianchi type I Universes in the Maxwell–Chern–Simons theory

In the present Section we will investigate the time dependence of the geometrical and thermodynamical parameters of the Bianchi type I space-times in the Maxwell–Chern–Simons theory. In order to obtain a relevant physical and cosmological picture we adopt for the description of the matter content a number of equations of state that could be relevant for the understanding of the properties of the ultra-high density matter the Universe may have contained in its very early stages.

In order to numerically integrate the evolution equations (86)–(89) we will fix the numerical values of the parameters \(\left( \lambda _1,\lambda _2\right) \). Once this is done, the numerical values of the free parameters of the model, \(\mu _0\) and \(\beta _0\), can be obtained from Eqs. (81) and (82) as

$$\begin{aligned} \mu _0=\frac{\sqrt{2-6 \lambda _2}}{\sqrt{4 \lambda _1+2 \lambda _2-1}}, \end{aligned}$$
(92)

and

$$\begin{aligned} \beta _0=\frac{\left( 4 \lambda _1+2\lambda _2-1\right) \left[ -\lambda \lambda _2+\lambda +\Lambda \left( 4 \lambda _1+2 \lambda _2-1\right) \right] }{6 \lambda _2^2+4 \lambda _2-2}, \end{aligned}$$
(93)

respectively. It is interesting to note that the numerical value of the coefficient \(\mu _0\) is determined by \(\lambda _1\) and \(\lambda _2\) only, while \(\beta _0\) is also determined by the arbitrary values of \(\lambda \).

5.1 Stiff fluid filled Bianchi type I Universe

An important equation of state, extensively used to describe the properties of high density matter, is the Zeldovich (or stiff matter) equation of state, which can be used for matter densities significantly higher than nuclear densities, \(\rho >10\rho _{n}\), where \(\rho _{n}\) is the nuclear density. The Zeldovich equation of state can be obtained theoretically from a relativistic Lagrangian that allows bare nucleons to interact attractively by exchanging a scalar meson, and to interact repulsively by exchanging a massive vector meson [148]. In the non-relativistic limit both the quantum and the classical description of strong interactions yield Yukawa-type potentials. At the highest matter densities the nuclear interactions are dominated by the vector meson exchange, and one can show, by using a mean field approximation, that in the extreme limit of infinite densities the pressure tends to the energy density, \(p\rightarrow \rho \). In this high density limit the speed of sound, given by \(c_{s}^{2}=\mathrm{d}p/\mathrm{d}\rho \rightarrow 1\). Therefore the stiff fluid equation of state satisfies the causality condition, which requires that the speed of sound is less than the speed of light, \(c_{s}\le c\). For the Zeldovich fluid with \(r=P\), the energy conservation gives the dependence of the density as a function of the comoving V in the form \(r=r_{0}/V^{2}\).

In the case of the stiff fluid the cosmological evolution equations take the simple form

$$\begin{aligned} \frac{\mathrm{d}V}{\mathrm{d}\tau }=u, \quad \frac{\mathrm{d}u}{\mathrm{d}\tau }=\left\{ 1+\lambda _{1}\frac{\left( B_{0}^{2}+\mu _{0}^{2}\right) }{1-\frac{1}{4}\left( B_{0}^{2}-\mu _{0}^{2}\right) }\right\} V, \end{aligned}$$
(94)
$$\begin{aligned} \frac{\mathrm{d}B_{0}}{\mathrm{d}\tau }= & {} b_{0},\quad \frac{\mathrm{d}b_{0}}{\mathrm{d}\tau }=\frac{1}{B_{0}}b_{0}^{2}+ \left[ \lambda _{1}\frac{1+\left( B_{0}^{2}-\mu _{0}^{2}\right) }{1-\frac{1}{ 4}\left( B_{0}^{2}-\mu _{0}^{2}\right) }-\lambda _{2}\right] B_{0}\nonumber \\&-\left[ 2 \sqrt{\frac{\lambda _{1}}{1-\frac{1}{4}\left( B_{0}^{2}-\mu _{0}^{2}\right) } }+\frac{u}{V}\right] b_{0}. \end{aligned}$$
(95)

The results of the numerical integration of the above system are presented in Figs. 1 and 2, respectively.

Fig. 1
figure 1

Variation of the comoving volume element V (upper left figure), of the mean dimensionless Hubble function h (upper right figure), of the dimensionless vector field \(B_0\) (lower left figure), and of the dimensionless matter energy density r for a Bianchi type I Universe filled with a stiff fluid in the Maxwell–Chern–Simons theory, for different values of the parameters \(\lambda _1\) and \(\mu _0\): \(\lambda _1=0.11\), \(\mu _0= 0.447\) (solid curve), \(\lambda _1 =0.12\), \(\mu _0=0.377\) (dotted curve), \(\lambda _1=0.13\), \(\mu _0=0.333\) (short dashed curve), \(\lambda _1 =0.14\), \(\mu _0=0.301\) (dashed curve), and \(\lambda _1 =0.15\), \(\mu _0=0.277\) (long dashed curve), respectively. The initial conditions used to numerically integrate the cosmological evolution equations are \(V (0)=0.8\), \(u(0)=4\), \(B_0 (0)=0.7\), and \(b_0(0)=-0.001\), respectively. In all cases the numerical value of the parameter \(\lambda _2\) has been fixed as \(\lambda _2=0.33\)

Full size image

In order to numerically integrate the cosmological evolution equations for the stiff fluid case we have fixed the value of the free parameter \(\lambda _2\) as \(\lambda _2=0.33\). As initial conditions we have adopted the values \(V(0)=0.8\), \(u(0)=4\), \(B_0 (0)=0.7\), and \(b_0(0)=-0.001\). The comoving volume element of the Bianchi type I Universe, presented in the left panel of Fig. 1, is a monotonically increasing function of the cosmological time \(\tau \), indicating an expansionary evolution of the Bianchi type I Universe. For small times the increase of V is almost linear, and the variation of the numerical values of the parameter \(\lambda _1\) influences the behavior of V only in the large time limit. The mean Hubble function, shown in the right panel of Fig. 1, is a monotonically decreasing function of time, and its behavior is slightly influenced in the large time limit by the variation of \(\lambda _1\). In the large time limit h tends to a constant value, indicating that in the presence of the Maxwell–Chern–Simons terms the stiff fluid filled Universe ends in an exponentially expanding de Sitter type era. The time evolution of the vector field \(B_0\), depicted in the left lower panel of Fig. 1 shows a strong dependence on the numerical values of \(\lambda _1\). \(B_0\) is a monotonically decreasing function of time, and in the large time limit it takes very small numerical values. The energy density of the matter, represented in the right lower panel of Fig. 1, is a monotonically decreasing function of time, whose evolution is practically independent on the variation of the numerical values of \(\lambda _1\). In the large time limit the matter energy density tends to zero, indicating that the de Sitter time expansion leads to a vacuum Universe, whose dynamics is dominated by the contributions of the Maxwell–Chern–Simons vector fields.

The time variations of the deceleration parameter of the Bianchi type I Universe, and of its anisotropy parameter, are represented in Fig. 2.

Fig. 2
figure 2

Variation of the deceleration parameter q (left figure), and of the anisotropy parameter A (right figure) for a Bianchi type I Universe filled with a stiff fluid in the Maxwell–Chern–Simons theory, for different values of the parameters \(\lambda _1\) and \(\mu _0\): \(\lambda _1=0.11\), \(\mu _0= 0.447\) (solid curve), \(\lambda _1 =0.12\), \(\mu _0=0.377\) (dotted curve), \(\lambda _1=0.13\), \(\mu _0=0.333\) (short dashed curve), \(\lambda _1 =0.14\), \(\mu _0=0.301\) (dashed curve), and \(\lambda _1 =0.15\), \(\mu _0=0.277\) (long dashed curve), respectively. The initial conditions used to numerically integrate the cosmological evolution equations are \(V (0)=0.8\), \(u(0)=4\), \(B_0 (0)=0.7\), and \(b_0(0)=-0.001\), respectively. In all cases the numerical value of the parameter \(\lambda _2\) has been fixed as \(\lambda _2=0.33\)

Full size image

The expansion of the Bianchi type I Universe begins in its very early stages with q having values around \(q\approx 2\). The Universe is initially in a decelerating state, with \(q>0\), but for \(\tau \approx 0.7\) the deceleration parameter reaches the value \(q\approx 0\), and the Universe enters in an accelerating phase. In the large time limit q reaches values of the order of \(q\approx -1\), indicating the presence of the de Sitter expansion. The overall evolution of q is slightly dependent on the numerical values of \(\lambda _1\). The time variation of the anisotropy parameter A, shown in the right panel of Fig. 2, is also strongly dependent on the model parameter \(\lambda _1\). In the large time limit \(A\rightarrow 0\), indicating that when the Universe enters the de Sitter phase it is already in an isotropic state. High values of \(\lambda _1\) lead to a rapid transition to the de Sitter era, as well as to the rapid isotropization of the Bianchi type I geometry.

5.2 The dust Bianchi type I Universe

As a next cosmological application of our model we will consider the case of the dust Universe, with the matter having negligibly small thermodynamic pressure, corresponding to the choice \(\gamma =1\). The overall dynamics is very similar to the stiff fluid case, with the Bianchi I type Universe isotropizing in the large time limit, and ending in a de Sitter type expansionary phase. In the following we will investigate the effect of the variation of the parameter \(\mu _0\) on the cosmological dynamics.

In order to numerically integrate the cosmological evolution equations in the case of the dust Universe we fix the numerical value of the free parameter \(\lambda _1\) as \(\lambda _1=0.1\), and we vary the numerical values of \(\mu _0\), and of \(\lambda _2\), as given by

$$\begin{aligned} \lambda _2=\frac{-4 \lambda _1 \mu _0^2+\mu _0^2+2}{2 \left( \mu _0^2+3\right) }. \end{aligned}$$
(96)

Similarly to the stiff fluid case, as initial conditions we adopt the values \(V(0)=0.8\), \(u(0)=4\), \(B_0 (0)=0.7\), and \(b_0(0)=-0.001\), respectively. The time variations of the comoving volume element, Hubble function, vector field, matter energy density, deceleration parameter and anisotropy parameter are plotted in Figs. 3 and 4, respectively.

Fig. 3
figure 3

Variation of the comoving volume element V (upper left figure), of the mean dimensionless Hubble function h (upper right figure), of the dimensionless vector field \(B_0\) (lower left figure), and of the dimensionless matter energy density r for a dust fluid filled Bianchi type I Universe in the Maxwell–Chern–Simons theory, for different values of the parameter \(\mu _0\): \(\mu _0 =0.10\) (solid curve), \(\mu _0 =0.20\), (dotted curve), \(\mu _0=0.30\), (short dashed curve), \(\mu _0 =0.40\) (dashed curve), and \(\mu _0 =0.50\) (long dashed curve), respectively. The initial conditions used to numerically integrate the cosmological evolution equations are \(V (0)=0.8\), \(u(0)=4\), \(B_0 (0)=0.7\), and \(b_0(0)=-0.001\), respectively. In all cases the value of the parameter \(\lambda _1\) has been fixed as \(\lambda _1=0.1\), while \(\lambda _2=\left( -4 \lambda _1 \mu _0^2+\mu _0^2+2\right) /2 \left( \mu _0^2+3\right) \approx 0.33\) for all considered values of \(\mu _0\)

Full size image

The comoving volume element of the dust fluid filled Bianchi type I Universe, plotted in the left panel of Fig. 3, is a monotonically increasing function of the cosmological time \(\tau \), indicating that the Bianchi type I Universe is globally expanding. For small time intervals the variation of V is almost linear. The mean Hubble function, depicted in the right panel of Fig. 3, is a monotonically decreasing function of time. The time variations of both V and h are practically independent on the modifications of the numerical values of the parameter \(\mu _0\), with a very slight influence manifesting itself only in the large time limit. On the other hand for large values of \(\tau \), h tends to a constant value, indicating that in the presence of the Maxwell–Chern–Simons terms the dust fluid filled Bianchi type I Universe experiences a transition to an accelerating phase, ending in a de Sitter type regime. The time evolution of the vector field \(B_0\), represented in the left lower panel of Fig. 3, indicates a strong dependence on the numerical values of \(\mu _0\). \(B_0\) is a monotonically decreasing function of time, and it continues to decrease even during the accelerated expansion of the Bianchi type I geometry, and in the isotropic phase. The energy density r of the dust matter, depicted in the right lower panel of Fig. 3, is a monotonically decreasing function of time, whose dynamics is essentially independent on the modifications of the numerical values of \(\mu _0\). In the large time limit the matter energy density tends to zero, \(\lim _{\tau \rightarrow \infty }r=0\), thus showing that the accelerated expansion of the dust Bianchi type I Universe in the presence of a Maxwell–Chern–Simons type field leads in its final stages to a vacuum Universe, in which the energy density of the ordinary matter gives a negligibly contribution to the total energy of the Universe.

The time variations of the deceleration parameter and of the anisotropy parameter of the dust Bianchi type I Universe in the presence of the Maxwell–Chern–Simons field are plotted in Fig. 4.

Fig. 4
figure 4

Variation of the deceleration parameter q (left figure), and of the anisotropy parameter A (right figure) for a dust fluid filled Bianchi type I Universe in the Maxwell–Chern–Simons theory, for different values of the parameter \(\mu _0\): \(\mu _0 =0.10\) (solid curve), \(\mu _0 =0.20\), (dotted curve), \(\mu _0=0.30\), (short dashed curve), \(\mu _0 =0.40\) (dashed curve), and \(\mu _0 =0.50\) (long dashed curve), respectively. The initial conditions used to numerically integrate the cosmological evolution equations are \(V (0)=0.8\), \(u(0)=4\), \(B_0 (0)=0.7\), and \(b_0(0)=-0.001\), respectively. In all cases the value of the parameter \(\lambda _1\) has been fixed as \(\lambda _1=0.1\), while \(\lambda _2=\left( -4 \lambda _1 \mu _0^2+\mu _0^2+2\right) /2 \left( \mu _0^2+3\right) \approx 0.33\) for all considered values of \(\mu _0\)

Full size image

The time variation of the deceleration parameter is presented in the left panel of Fig. 4. For the chosen numerical values of the model parameters the expansion of the Bianchi type I Universe starts with values of q of the order of 1.5, \(q\approx 1.5\) at \(\tau =0\). The dust Bianchi I Universe is in a decelerating phase for small values of \(\tau \), with \(q>0\). However, for \(\tau =\tau _{cr}\approx 0.9\), the deceleration parameter reaches its borderline value \(q\approx 0\), and for \(\tau >\tau _{cr}\) the Universe enters in an accelerating phase, with \(q<0\). In the large time limit q reaches values of the order of \(q\approx -1\). The time evolution of q is basically independent on the numerical values of \(\mu _0\). The time variation of the anisotropy parameter A, depicted in the right panel of Fig. 4, does show a mild dependence on the numerical values of the parameter \(\mu _0\). In the large time limit A becomes zero, \(A\rightarrow 0\). This situation corresponds to values of q of the order of \(q\approx -1\), which shows that the Universe is already isotropic before entering the de Sitter phase. The large time behavior of A, before the full isotropization period, shows a small dependence on the numerical values of \(\mu _0\).

6 Discussions and final remarks

In this paper we have considered the effects of the topological Chern–Simons term on the cosmological evolution of the Universe. The original Chern–Simons term consists of a constant vector field, which breaks a subset of the Poincaré group containing the time translation and boost. In order to restore the whole Poincaré symmetry, we have promoted this constant vector to a dynamical vector field, with constant norm, through a condition which is imposed in the action by a Lagrange multiplier term. Another example of a gravity theory with two dynamical vector fields was considered in [96] where the Weyl vector and the trace of torsion tensor are two vector fields coupled minimally to gravity. The norm of the vector field \(A_\mu \) in our model could in principle be positive or negative or even zero. However, in order to respect the spatial isotropy of the Universe in the FRW space-time this constrained vector field should be time-like. The theory has a self-accelerating de Sitter solution. We have performed the cosmological perturbation analysis of the theory around this de Sitter background, and we have found that the theory has two gravitational wave modes, together with four vector modes, corresponding to the two vector fields of the theory. The vector modes of the theory interact with each other through the Chern–Simons interaction term. This interaction makes the vector degrees of freedom of the theory non-trivial. Also there is one scalar mode in this theory, leaving seven degrees of freedom around the de Sitter background. It should be mentioned that the constrained vector field \(A_\mu \) has three dynamical degrees of freedom around the de Sitter background, since we have a Lagrange multiplier term which breaks the U(1) gauge invariance for \(A_\mu \). This is similar to the case Einstein-aether theory where a constrained time-like vector field is added to the Einstein–Hilbert action. The main difference is that here we have a gravity theory with two dynamical vector field which couple through a non-trivial Chern–Simons interaction term. The perturbation analysis of the theory then shows that the theory is stable and healthy around de Sitter background. Of course we have to perform Hamiltonian analysis in order to find the exact number of degrees of freedom of the theory.

In order to investigate the cosmological implications of the space-like constrained vector field one should consider an anisotropic space-time. In this paper, we have considered the dynamical behavior of a Bianchi type I space-time in the framework of the Maxwell–Chern–Simons gravity model. The Bianchi type I anisotropic space-time represents the simplest, and most natural, extension of the standard FRW metric, to which it reduces in the particular limit of equal directional scale factors. Bianchi type models can be considered as viable alternatives to the standard flat FRW cosmologies [149,150,151,152]. The small observed deviations from the exact isotropy and the anomalies in the Cosmic Microwave Background could be explained by the presence of an anisotropic expansion of the Universe. Bianchi type VIIh anisotropic cosmological models were considered in [153], in a tentative to explain the large scale asymmetry observed in the Cosmic Microwave Background distribution. An in-depth comparison with the WMAP first-year data on large angular scales did show that a chance alignment can be eliminated at a 3\(\sigma \) level. On the other hand, the recent Planck Collaboration results [154] did show convincingly that the Bianchi type VIIh anisotropic cosmological model cannot fit the recent observational data obtained by the Planck satellite. However, one of the large angle anomalies of the Cosmic Microwave Background, the low quadrupole moment, indicates a great amount of power suppression at large scales. In this context it is important to note that this anomaly is within the cosmic variance. On the other hand the presence of such an anomaly, if confirmed by further observations, seems to indicate the existence of a Bianchi type I anisotropic geometry of the Universe, which could indeed alleviate the low quadrupole moment problem [155,156,157,158]. The extreme smallness of the quadrupole component of the Cosmic Microwave Background temperature distribution seems to suggest that for a homogeneous but anisotropic Universe the deviation from the isotropic flat FRW geometry must be small. Therefore such a deviation can be naturally explained by the existence of a background Bianchi type I geometry. However, it is important to point out that in order to alleviate the low quadrupole moment, the anisotropy must have a specific and well-determined evolution. Moreover, the presence in the Universe of an anisotropic geometry could also raise the quadrupole moment, and this is more likely to happen during the cosmological dynamics [159]. It is also important to point out that cosmological observations do not favor the Bianchi type I geometry [154].

In the present paper we have found that in the framework of Maxwell–Chern–Simons gravity theory the Bianchi type I homogeneous but anisotropic Universe presents a complex dynamics. In our analysis we have assumed that the matter content of the Universe consists of a perfect barotropic cosmological fluid. In particular, for this type of matter source, the Bianchi type I models we have considered do always isotropize. The nature of the cosmological evolution strongly depends on the dimensionless model parameters \(\lambda _1\) and \(\lambda _2\), as well as on the adopted initial conditions for the matter energy density, Hubble function, and of the vector field itself. The transition to an isotropic phase of the Universe is associated with an accelerated, de Sitter type expansion, in which the considered cosmological models end in the large time limit. This type of behavior is independent on the nature of the cosmological fluid, and it does appear in the two cases of stiff and dust fluids, respectively.

Hence, the inclusion in the gravitational action of two vector fields, coupled via a Chern–Simons term, leads to the possibility of obtaining more general gravitational models, thus allowing for the possibility of a better and more realistic physical description of both the very early and the late evolution of our Universe.