The Maxwell-Chern-Simons gravity and its cosmological implications

We consider the cosmological implications of a gravitational theory containing two vector fields coupled minimally to gravity as well as a generalized Chern-Simons term that couples the two vector fields. One of the vector fields is the usual Maxwell field, while the other is a constrained vector field with constant norm included in the action via a Lagrange multiplier. The theory admits a de Sitter type solution, with healthy cosmological perturbations. We will show that there is 6 degrees of freedom propagate on top of de Sitter space-time, two tensor polarizations and four degrees of freedom related to two massless vector fields interacting with each other via Chern-Simons interaction term. We also investigate in detail the behavior of the geometric and physical parameters of a homogeneous and anisotropic Bianchi type I Universe, by using both analytical and numerical methods, by assuming that the matter content of the Universe can be described by the stiff causal and pressureless dust fluid equations of state. The time evolution of the Bianchi type I Universe strongly depends on the initial conditions of the physical and geometrical quantities, as well as on the numerical values of the model parameters. Two important observational parameters, the mean anisotropy parameter, and the deceleration parameter, are also studied in detail, and we show that independently of the matter equation of state the cosmological evolution of the Bianchi type I Universe always ends in an isotropic and exponentially accelerating, de Sitter type, phase.


I. INTRODUCTION
One of the most interesting findings of modern cosmology is that our Universe is in a phase of accelerated expansion [1]. To explain this fundamental observation, one has to adopt a theory of gravity that admits a self-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 phenomenological problems, such as, for example, the fine tuning problem [2]. 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 [3], or, more generally, in f (R, Q) theories where Q can be any scalar combination of the curvature tensor/energy-momentum tensor [4]. The second possibility is enriching the graviton itself, such that it admits more than 2 degrees of freedom [5], and constructing the massive gravity theory. In the latter case, the accelerated expansion of the Universe can be obtained not by imposing a repulsive gravity (like the cosmological constant), but with the reduction of attractive gravity.
An alternative way to produce a self-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 + ρ ≈ 0. This can be achieved by adding some higher derivative kinetic terms [6], or by adding a potential term of the scalar field to the theory [7]. Many works have been done in the context of scalar field cosmology, including inflationary [8] and dark energy models [9]. 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 [10]. The cosmological implications of the Horndeski type theories have been extensively investigated in the literature [11].
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 [12]. One can also consider a massive vector field as an Einstein-Proca system. Many works has been done in the context of Einstein-Maxwell and Einstein-Proca theories [13]. Another possibility is to consider the 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 ∇ µ g νρ = w µ g νρ [14,15]. One can also obtain the Einstein-Proca system from the higher dimensional Gauss-Bonnet theory in the context of Weyl geometry [16]. Cosmological implications of such a theory has been considered in [17]. 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 [18], 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 helicity-0 of the vector field mimics the scalar galileon interactions [19]. The cosmology of this vector galileon theory, as well as possible generalization, have been investigated in the literature [20,21].
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 [22].
Other possibility of generalizing Einstein's general relativity within a purely geometric approach is to consider the effects of the torsion tensor [23]. The torsion tensor can in general be decomposed into a vector field Q µ , a axial vector field S µ and a tensor field t µνρ satisfying the conditions t µνρ + t νρµ + t ρµν = 0, and g µν t µνρ = 0 = g µρ t µνρ , respectively. Assuming that the tensor t µνρ 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 [24].
An interesting vector-tensor theory of gravity is the Einstein-aether theory [25]. 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 [26], as well as its relation to the scalar-tensor theory [27], 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 [28] where F a µν = ∇ µ A a ν − ∇ ν A a µ , a = 1, 2, 3, ∇ µ represents the covariant derivative with respect to the metric, while V (A 2 ) with A 2 = g µν A a µ A 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 [29].
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, [30]. The action for the non-minimally massive vector field coupled to gravity can be represented as where µ Λ is the mass of the massive cosmological vector field, and A µ (x ν ), µ, ν = 0, 1, 2, 3 is its four-potential, which is allowed to couple non-minimally to gravity. Here ω and η are dimensionless coupling parameters. By following a close analogy with electrodynamics, the dark energy vector type field tensor is defined again as F µν = ∇ µ A ν − ∇ ν A µ . A so called superconducting type dark energy model was recently introduced in [31]. 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 where V (φ) 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 where λ and α are arbitrary constants, L m = L m (g µν , ψ) is the Lagrangian of the total (ordinary baryonic plus dark) matter, and j µ = ρu µ is the total mass current. In Eq. (3) ρ denotes the total matter density (including the dark matter one), while u µ 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 [32] where B µν = ∂ µ B ν − ∂ ν B µ , and k µ 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 [33]. 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 µ , the Lorentz invariance will be broken by this term. This implies that the Chern-Simons term will reduce the Lorentz invariance symmetry to the 3 dimensional translational symmetry. Many works has been done in the context of Chern-Simons electrodynamics, including applications to quantum electrodynamics [34]. In [35] the authors obtain the Chern-Simons term from dimensional reduction of the Carroll-Field-Jackiw Higgs model. Also in [36] 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 [37][38][39]. 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 anomalycanceling mechanism [40]. Two distinct formulations of Chern-Simons modified gravity have been proposed, namely, the nondynamical formulation and the dynamical formulation (see [39] 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 [41]. 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 promote the constant vector k µ to a dynamical field and impose a constraint its norm remains constant by adding this property by a Lagrange multiplier. This will restore the Lorentz invariance of the theory. The Chern-Simons term will be generalized to where λ is a Lagrange multiplier, η and α are constant and ǫ µνρσ is the Levi-Civita tensor. The new vector field A µ plays the role of the constant vector k µ 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 µ be purely spatial and has a form k µ = (0, k) [42]. In this paper, we will not interested to 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 µ . Also, as we will be see in the following, in the case of FRW cosmology, one should assume that the vector field be time-like and has only (t)-component, because otherwise the spacial 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 µ 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 Section III 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 Section IV we will assume that the vector A µ 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 Section V 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 and conclude our results.

II. THE MODEL
We propose an action functional of the form where A µ and B µ are two vectors fields, A 2 = A µ A µ , S m is the matter action, η is the norm of vector field A µ with dimension of mass squared and λ is the Lagrange multiplier which dynamically enforces that the vector field A µ has a constant norm. This action consists of an Einstein-Maxwell system coupled to a constant norm vector field A µ . The two vector fields then interact via the topological Chern-Simons term.
Varying the action with respect to the Lagrange multiplier results in which states that the vector field A µ should have a constant norm. The equation of motion for B µ can be written as The equation of motion for A µ is where we have defined V ′ = dV /dA 2 . One can see that the A µ field is sourced by the B µ vector field through the Chern-Simons term, and acquires an effective mass 2λ. Also B µ vector field is sourced by the A µ field through the Chern-Simons term. In the case of vanishing α, the two vector field evolve independently. The equation of motion for the metric is In order to obtain the conservation of energy-momentum tensor in this model, first note that from the constraint equation (6), one has A µ ∇ α A µ = 0. Also, by taking the divergence of equation (8), one can obtain Now, taking the derivative of equation (9), we arrive at This shows that the energy-momentum tensor is not conserved due to non-linear interaction between two vector fields.

III. 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 where Λ is the cosmological constant.

A. de Sitter solution
Let us assume that our Universe can be described by a flat FRW metric of the form Moreover, because we are interested in the self-accelerated solutions to the theory, we assume that the energymomentum tensor vanishes. In this case, the most general ansatz for the vector field which respects the spatial isotropy is From the above equation, one can see that the vector field A µ should be time-like in this case, and we will assume that η = m 2 > 0. For the space-like vector field η < 0, there is no FRW solution unless A µ = 0 which contradicts equation (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 µ equation of motion (7) satisfies automatically in this case. The equations of motion for the vector field A µ and for the metric can then be written as which can be solved as We will also assume that B 0 = constant for the sake of simplicity. In order to have a consistent de Sitter solution, one should has β < 2κ 2 Λ/m 4 . Note that, in the absence of cosmological constant Λ, the coupling constant β should be negative in order to have a de Sitter solution.

B. 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 where φ, ψ, E and B are the scalar perturbations, F i and S i are the vector perturbations with the property ∂ i F i = 0 = ∂ i S i , and h ij is associated with the tensor perturbation which is transverse and traceless h ii = 0 = ∂ i h ij . Also, the spatial indices are raised and lowered by δ ij . For the vector fields, we decompose the perturbed vector field as where δA 0 , δB 0 , δA and δB are the scalar perturbations and ξ i and ǫ i are the vector perturbations with the property Under the general linear coordinate transformation x µ → x µ + δx µ , the metric perturbations transforms as where we have decomposed δx µ to (δx 0 , η i + ∂ i δx), with ∂ i η i = 0. Also the various components of a vector field A µ will be transformed as Likewise, we have a similar transformation rule for the vector field B µ . Also, one can find that the Lagrange multiplier will not change under this transformation We should note that in obtaining the above relations we have assumed that the background values A 0 , B 0 and λ 0 are constant, as is obtained from the previous section.
With the aid of above transformations, one can construct 7 independent gauge invariant scalar perturbations, 3 independent gauge invariant vector perturbations and one gauge invariant tensor perturbation as follows and the perturbation of the Lagrange multiplier is already gauge invariant. The vector perturbations ξ i and ǫ i are already gauge invariant and the remaining gauge invariant vector perturbation can be constructed as 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 to (5) and expand the resulting action up to second order in perturbations, one can see that the tensor, vector and scalar parts is decoupled completely from each other. In the following we will consider them separately.

Tensor perturbation
The tensor perturbation h ij is transverse and traceless, and can be described by two polarization modes h + and h × . After Fourier decomposition, one can obtain the second order action of tensor perturbation as where k is the comoving wave vector and dot represents time derivative. Also, it is understood that in any term in the second order action, the argument of one of the perturbation variables is k, and for the other variable is −k, due to the definition of Dirac delta function. One can see that the Chern-Simons coupling term does not contribute to the tensor perturbation. The action (30) is equivalent to that of Einstein-Hilbert theory. So, the theory has two tensor propagating modes associated to the massless graviton.

Vector perturbation
We have 3 transverse independent gauge invariant vector perturbations β i , ξ i and ǫ i . After expanding the action up to second order in perturbations and performing the Fourier transformation, one can obtain where we have added the arguments of the last two terms for clarity. The vector perturbation of the metric β i does not couple to the other vector fields and is non-dynamical. Obtaining the equation of motion of β i and substituting it back to the action, one can see that β i will be vanishes from the action. We have left with 2 vector perturbations associated with two vector fields A µ and B µ . One can see that these vector perturbations interact non-trivially with each other through the Chern-Simons coupling, together with a self interaction of the vector field B µ . At the end, we will have 4 healthy vector degrees of freedom for the theory.

Scalar perturbation
In this section we will obtain the scalar perturbation part of the theory. There are 5 scalar modes, corresponding to the metric, the vector fields and the Lagrange multiplier. Expanding the action up to second order and performing Fourier decomposition, one can write the action with gauge invariant quantities as One can see from the above equation that δλ, δB 0 , δA 0 , and Φ are non-dynamical. Varying the action with respect to these variables leads Solving the above equations for non-dynamical variables and substituting them back to the action gives where we have defined a new variable Ξ = Hδ A − mΨ. Varying the action with respect to Ξ results in Upon substituting the above relation to the action, one obtains which can be solved for Ψ with the result Ψ = 0. So, the scalar part of the perturbed action vanishes and one can conclude that there is no dynamical scalar mode in the theory.
In summary, we have 2 gravitational wave mode and 4 vector modes associated with two vector fields in the theory. One should mention the Lagrange multiplier constraint implies that the vector field A µ acquire only two degrees of freedom. This is similar to the case of Einstein-aether theory [25].

IV. ANISOTROPIC COSMOLOGY: BIANCHI TYPE-I UNIVERSE
In this section, we want to consider the case that the vector field A µ is space-like. So, from now on we will set η = −µ 2 . In this case, at least one of the spatial components of A µ 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 µ . In this case the vector field A µ has the from In this case the vector field B µ may have any direction in space and so we set The above ansatz for the vector field will break the spacial 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 where a i 's are the directional scale factors. We also assume that the matter content of the Universe has a form where ρ 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 [43] H From Eqs. (44) we obtain immediately the important relatioṅ The anisotropy parameter can be represented in an equivalent form as Before solving the cosmological equations, we should note that the Chern-Simons term can be considered as an effective cosmological constant, as long as A i 's and B i 's do not change very rapidly. This is because the Levi-Civita tensor is proportional to the volume of the Universe V , so that the epsilon terms in the action varies very slowly, mimicking the cosmological constant. However, from the vector field equations (7) and (8), one can see that the source terms (which are proportional to the Chern-Simons term) decay at late times.

A. Gravitational field equations
Before writing the equations of motion of the theory, we should note that the (t)-component of the vector field B µ 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 µ equation of motion (8), implies that A 0 = 0. So, the constraint equation (6) gives the remaining component of the vector field A µ as A 1 = µa 1 (t). Also, the off-diagonal elements (ij), i = j, of the metric field equation givesḂ iḂj = 0. This implies that at least two out of three components of the vector field B µ should be constant. We will assume that B 1 and B 2 is constant. Now, the (y)-component of the A µ field equation gives B 1Ḃ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 µ evolves in time. With these in hand, the only remaining component of the A µ field equation becomes The remaining components of the B µ equation of motion are or, equivalently, and The Friedmann equation can be written as and the Raychaudhuri equations are and In the above equations we have defined Equation (52) can be solved for b 3 with the result Also from equation (11), one can show that the matter field is conserved in this casė

B. 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. (55)-(57) we obtain By substituting 3Ḣ from the above equation into Eq. (54) we find the consistency condition With the use of Eq. (51), Eq. (50) becomes By combining Eqs. (55) and (63) we find With the help of Eq. (52), Eq. (53) can be reformulated as giving, after substitution into Eq. (57), the equation Adding Eqs. (64) and (66) we obtain while equating Eqs. (64) and (66) gives From Eq. (61) we obtain immediately the deceleration parameter as while from Eq. (62) we obtain the anisotropy parameter in the form C. 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. (64) we obtain for H 1 the expression where we have denoted By substituting this expression of H 2 1 into Eq. (61) we obtain the evolution equation for V as Eq. (66) gives the time evolution of the field b 3 as

D. 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 (τ, r, P, B 0 , µ 0 , β 0 , h), defined as Then the system of Eqs. (73) and (74) becomes where and The energy conservation equation can be written as dr dτ + 3h (r + P ) = 0.
By assuming an equation of state of the form we obtain for the energy density of the cosmological matter the expression 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 pres , is given by ρ (t pres ) = 4κ 2 Λ ef f /V γ (t pres ). Therefore the system of equations (76) and (77), 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 The system of differential equations (83)-(86) 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, λ 1 and λ 2 , representing the combination of the four free model parameters (Λ, λ, β, µ).
The deceleration parameter can be written as From Eq. (54) we obtain the following relation between the deceleration and the anisotropy parameters Hence the behavior of both q and A is determined by the set of model parameters (µ 0 , λ 1 , λ 2 , r 0 ).

V. 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 (83)-(86) we will fix the numerical values of the parameters (λ 1 , λ 2 ). Once this is done, the numerical values of the free parameters of the model, µ 0 and β 0 , can be obtained from Eqs. (78) and (79) as and respectively. It is interesting to note that the numerical value of the coefficient µ 0 is determined by λ 1 and λ 2 only, while β 0 is also determined by the arbitrary values of λ.

A. 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, ρ > 10ρ n , where ρ 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 [44]. In the non-relativistic limit both the quantum and 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 → ρ. In this high density limit the speed of sound, given by c 2 s = dp/dρ → 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 ≤ 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 The results of the numerical integration of the above system are presented in Figs. 1 and 2, respectively. In order to numerically integrate the cosmological evolution equations for the stiff fluid case we have fixed the value of the free parameter λ 2 as λ 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 τ , 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 λ 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 λ 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 λ 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 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 λ 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. The expansion of the Bianchi type I Universe begins in its very early stages with q having values around q ≈ 2. The Universe is initially in a decelerating state, with q > 0, but for τ ≈ 0.7 the deceleration parameter reaches the value q ≈ 0, and the Universe enters in an accelerating phase. In the large time limit q reaches values of the order of q ≈ −1, indicating the presence of the de Sitter expansion. The overall evolution of q is slightly dependent on the numerical values of λ 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 λ 1 . In the large time limit A → 0, indicating that when the Universe enters the de Sitter phase it is already in an isotropic state. High values of λ 1 lead to a rapid transition to the de Sitter era, as well as to the rapid isotropization of the Bianchi type I geometry.

B. 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 γ = 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 µ 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 λ 1 as λ 1 = 0.1, and we vary the numerical values of µ 0 , and of λ 2 , as given by 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 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 µ 0 . In the large time limit the matter energy density tends to zero, lim τ →∞ 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.
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 ≈ 1.5 at τ = 0. The dust Bianchi I Universe is in a decelerating phase for small values of τ , with q > 0. However, for τ = τ cr ≈ 0.9, the deceleration parameter reaches its borderline value q ≈ 0, and for τ > τ cr the Universe enters in an accelerating phase, with q < 0. In the large time limit q reaches values of the order of q ≈ −1. The time evolution of q is basically independent on the numerical values of µ 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 µ 0 . In the large time limit A becomes zero, A → 0. This situation corresponds to values of q of the order of q ≈ −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 µ 0 .

VI. 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 the Lorentz invariance of the theory. In order to restore the Lorentz 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. In this way we have obtained a bi-vector tensor theory of gravity. Another example of a bi-vector gravity theory was considered in [24] 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 µ 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 2 gravitational wave modes, together with 4 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 no scalar mode in this theory, leaving 6 degrees of freedom around the de Sitter background. It should be mentioned that the constrained vector field A µ should have 3 dynamical degrees of freedom around the de Sitter background, since we have a self-interaction potential in the action. However, because of the Lagrange multiplier term, the helicity-0 mode will become non-dynamical in this case. 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 bi-vector tensor theory with 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 [45]. 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 [46], 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σ level. On the other hand, the recent Planck Collaboration results [47] 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. The presence of such an anomaly seems to indicate the presence of a Bianchi type I anisotropic geometry of the Universe. 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.
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 λ 1 and λ 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 both the cases of the 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 physical