Vector dark energy models with quadratic terms in the Maxwell tensor derivatives

We consider a vector-tensor gravitational model with terms quadratic in the Maxwell tensor derivatives, called the Bopp-Podolsky term. The gravitational field equations of the model and the equations describing the evolution of the vector field are obtained and their Newtonian limit is investigated. The cosmological implications of a Bopp-Podolsky type dark energy term are investigated for a Bianchi type I homogeneous and anisotropic geometry for two models, corresponding to the absence and presence of the self-interacting potential of the field, respectively. The time evolutions of the Hubble function, of the matter energy density, of the shear scalar, of the mean anisotropy parameter, and of the deceleration parameter, respectively, as well as the field potentials are obtained for both cases by numerically integrating the cosmological evolution equations. In the presence of the vector type dark energy with quadratic terms in the Maxwell tensor derivatives, depending on the numerical values of the model parameters, the Bianchi type I Universe experiences a complex dynamical evolution, with the dust Universes ending in an isotropic phase. The presence of the self-interacting potential of the vector field significantly shortens the time interval necessary for the full isotropization of the Universe.


I. INTRODUCTION
Recent cosmological observations, based initially on the study of the distant Type Ia Supernovae, have shown that the cosmological paradigm according to which the Universe must decelerate due to its own gravitational attraction is not correct and that the Universe has experienced a transition to a late time, de Sitter type accelerated phase [1][2][3][4]. These observations have triggered a deep revision of our understanding of the cosmological dynamics and of its theoretical basis, general relativity. To explain the current observations in cosmology many theoretical ideas and suggestions have been put forward to address the intriguing facts revealed by the complex observational study of the Universe. From both theoretical and observational points of view there is a general consensus, which may be referred to as "the standard explanation of the late time acceleration," according to which observations can be easily explained once we assume the existence of a mysterious and dominant component in the Universe, called dark energy, and is fully responsible for the observed dynamics in the late phases of the evolution of the Universe.
Another important cosmological result, based on the combination of data from the observations of high redshift supernovae, the WMAP satellite, and the recently released Planck data, convincingly show that the location of the first acoustic peak in the power spectrum of the CMBR (Cosmic Microwave Background Radiation) is entirely consistent with the important prediction of the inflationary model for the total density parameter Ω of the Universe, according to which at the end of inflation Ω = 1. The important cosmological parameter w = p/ρ, where p is the total pressure and ρ is the total density of the Universe is also strongly constrained by cosmological observations which provide detailed evidence for the behavior of the equation of state of the cosmological fluid, constraining the parameter w as lying in the range −1 ≤ w < −1/3 [5].
These large number of cosmological observations have led to the formulation of the ΛCDM paradigm according to which, in order to explain the cosmological evolution, one assumes that the Universe is filled with two main components representing around 95% of its content; cold (pressureless) dark matter (CDM) and dark energy (DE), having a negative pressure. The contribution of the CDM component to the total density parameter of the Universe is of the order of Ω m ∼ 0.3 [6]. From a theoretical point of view the necessity to consider dark matter is mainly required by the necessity of explaining the unusual behavior of the galactic rotation curves as well as the formation of the large scale structure. On the other hand, dark energy is considered as representing the major component of the "chemical" composition of the Universe, giving a contribution to the total density parameter of the order of Ω DE ∼ 0.7. Dark energy is the major cause determining the recent, de Sitter type, acceleration of the Universe, as confirmed by the study of high redshift type Ia supernovae [5]. The search for an explanation of the physical (or geometrical) nature and properties of dark energy has opened a very active field of research in cosmology and theoretical physics which in turn has led to a myriad of different DE models, for reviews of DE models see, for example, [7][8][9][10][11][12][13][14].
One interesting possibility for explaining DE which has been intensively investigated is based on a number of cosmological models in which the "chemical" composition of the Universe consists of a mixture of two major components; cold dark matter and a slowly-varying, spatially inhomogeneous component, called the quintessence [15]. In this scenario the baryonic matter plays a negligible role with a minimal influence on late time cosmological dynamics. From a formal theoretical point of view and based on some particle physics results, quintessence type cosmological models can be implemented by assuming that dark energy is the energy associated with a scalar field Q with a self-interacting potential V (Q) [16]. During the cosmological evolution, when the potential energy density V (Q) of the quintessence field becomes greater than the kinetic energy density, the thermodynamic pressure p =Q 2 /2 − V (Q), associated with the quintessence Q-field, becomes negative. The cosmological and astrophysical properties of the quintessential cosmological models have been intensively investigated in the literature, for a recent review see [17]. Quintessence models differ from cosmological models of the standard general relativity including the cosmological constant since they imply that the equation of state of the quintessence field varies dynamically with cosmic time [18]. A number of alternative cosmological models, called k−essence, where the late-time acceleration of the Universe is driven by the kinetic energy of the scalar field have also been proposed [19].
Another possibility to explain the recent acceleration of the Universe and the nature of dark energy is provided by scalar fields φ that are minimally coupled to gravity via a negative kinetic energy. An interesting property of these fields is that they allow for values of the equation of state parameter, w, of dark energy to vary in such a way as to have w < −1. These types of scalar fields are known as phantom fields, proposed as an explanation for the late time acceleration of the Universe in [20]. The energy density and the pressure of phantom scalar fields are given by , respectively. Phantom cosmological models for dark energy have been investigated in detail in [21][22][23][24]. Some recent cosmological observations seem to support the interesting result that at some instant during the evolution of the Universe the value of w representing dark energy equation of state may have crossed the standard value w = −1, hence entering a de Sitter type expansion with a cosmological constant Λ. This intriguing cosmological phenomenon is called the phantom divide line crossing [23]. The crossing of the phantom divide line was investigated in the case of scalar field models with cusped potentials in [22]. The phantom divide line crossing can also be explained in cosmological models where dark energy is represented by a scalar field, which is non-minimally coupled to gravity [22].
A different line of research on dark energy is based on the assumption that instead of interpreting dark energy as a specific physical field, the cosmological dynamics of the Universe can be understood as a modification of the gravitational force itself. By following this line of thought one can assume that at very large cosmological scales general relativity cannot describe the dynamical evolution of the Universe, and therefore the acceleration of the Universe is related to an intrinsic change of the gravitational interaction. A plethora of modified gravity models, based on different extensions of general relativity, like, for example, f (R) gravity (in which the gravitational action is an arbitrary function of the Ricci scalar R) [25] and mimetic-f (R) gravity models [26], the f (R, L m ) model (where L m is the matter Lagrangian) [27], f (R, T ) modified gravity models (where T denotes the trace of the energy-momentum tensor) [28], the Weyl-Cartan-Weitzenböck (WCW) model [29], hybrid metric-Palatini f (R, R) gravity models (where R is the Ricci scalar formed from a connection independent of the metric) [30], f (R, T, R µν T µν ) type models, where R µν is the Ricci tensor and T µν is the matter energy-momentum tensor, respectively [31], the Eddington-inspired Born-Infeld theory [32], f (T , T ) gravity [33], implying coupling between torsion scalarT and trace of the matter energy-momentum tensor, or vector Gauss-Bonnet theory [34], have been recently proposed in the literature. The cosmological and astrophysical properties of these models have been extensively investigated. For a recent review of the generalized gravitational models with non-minimal curvature-matter coupling f (R, L m ) and f (R, T ) type see [35]. For a review of hybrid metric-Palatini gravity see [36]. Modified gravity models can provide convincing theoretical explanations for the late time acceleration of the Universe without advocating the existence of dark energy and can also offer some alternative explanations for the nature of dark matter.
From a field theoretical point of view however, despite the great success of the scalar field dark energy models, the possibility that dark energy has a more complex structure than allowed by the simple scalar field model cannot be ignored a priori. One promising direction in the analysis of dark energy is represented by models in which dark energy is described by a vector or Yang-Mills type field which may also couple, minimally or non-minimally, to gravity. The simplest action for a Yang-Mills type dark energy model is [37,38] where A a µ , a = 1, 2, ..., n are the potentials of the Yang-Mills field, F a µν = ∇ µ A a ν − ∇ ν A a µ , ∇ µ is the covariant derivative with respect to the metric, A 2 is defined as A 2 = g µν A a µ A a ν and V (A 2 ) represents a self-interacting potential, explicitly violating gauge invariance. In the action given by Eq. (1) there are three vector fields describing dark energy. Hence Eq. (1) generalizes the Einstein-Maxwell type single vector field dark energy model. The astrophysical and cosmological applications of the single or Yang-Mills type vector dark energy models have been comprehensively investigated in [39].
Extended vector field dark energy models where the vector field is non-minimally coupled to the gravitational field can also be constructed [40]. The action for such a non-minimally coupled vector dark energy model is given by where A µ (x ν ), µ, ν = 0, 1, 2, 3 is the four-potential of the vector type dark energy, which couples non-minimally to gravity, and µ Λ is the mass of the massive cosmological vector field, respectively. The constants ω and η are dimensionless coupling parameters, while the vector dark energy field tensor is defined as Inspired by the possible analogy between dark energy and some condensed matter concepts, a so called superconducting type dark energy model was proposed in [41]. This model describes the spontaneous breaking of U(1) symmetry of the "electromagnetic" type dark energy, and is described by the action where λ and α are constants, L m (g µν , ψ) is the Lagrangian of the total (ordinary baryonic plus dark) matter, and j µ = ρu µ is the total mass current, where ρ is the total matter density (including dark matter), and u µ is the matter four-velocity. This model can also be interpreted and understood as unifying, in a single formalism, the scalar and vector dark energy models. The predictions of the superconducting dark energy model have been compared with observations in [42]. It is the goal of the present paper to consider a vector-tensor type model of dark energy, based on the analogy with Bopp-Podolsky electrodynamics. The Bopp-Podolsky theory was first suggested by Bopp [43], and was independently reobtained by Podolsky [44]. The Bopp-Podolsky theory retains linearity of the field equations but introduces higherderivative terms proportional to the parameter m 2 where m, having the physical dimensions of mass, is a new hypothetical fundamental constant of Nature. For m → ∞ the Maxwell-Lorentz theory, and the Maxwell equations, are retained. The Bopp-Podolsky theory is formulated in terms of an action functional from which the field equations, which are of fourth order in the electromagnetic potential, are derived. However, as noted by both Bopp and Podolsky, in a certain gauge these fourth-order equations are equivalent to a pair of second-order equations [43,44]. Different aspects of the Bopp-Podolsky type extension of classical electrodynamics were investigated in [45].
To this and other ends, we start from the analogy with the Bopp-Podolsky electrodynamics and introduce a vectortensor gravitational model where the action for the minimally coupled vector field also contains additional terms, quadratic in the Maxwell tensor derivatives. These terms correspond to the covariant form of the action of the Bopp-Podolsky electrodynamics. Moreover, a term describing the non-minimal coupling between the cosmological mass current and the four-potential of the vector field is also added to the action. The possible existence of a selfinteraction potential of the vector field is also considered. From a cosmological point of view we propose to interpret the vector field as describing the dark energy component of the Universe, which is responsible for the late, de Sitter type acceleration of the Universe. We obtain the gravitational field equations of this vector dark energy model as well as the equations describing the evolution of the vector field. We investigate the Newtonian limit of the model and show that the Poisson equation as well as the Bopp-Podolsky electrodynamics can be recovered for weak fields.
The cosmological implications of this vector type dark energy model are investigated for a Bianchi type I homogeneous and anisotropic geometry. Two cases are investigated in detail, the evolution of the Universe with and without the self-interacting potential of the field, respectively. In both cases the evolution of the Hubble function, of the matter energy density, of the shear scalar, of the anisotropy parameter, of the deceleration parameter, and of the field potential are analyzed in detail. To escribe the matter content of the Universe we adopt the radiation fluid and the dust matter equations of state. We find that in the presence of the vector type dark energy with quadratic terms in Maxwell tensor derivatives the anisotropic Universe experiences a complex dynamical behavior, with the dust Universes ending in an isotropic stage, a result which is independent on the presence or absence of the self-interaction potential of the field.
The present paper is organized as follows. The field equations of the Bopp-Podolsky type vector-tensor gravitational model are derived in Section II and their Newtonian limit is also investigated. The cosmological implications of the model are investigated in Section III, where the cosmological dynamics of a Bianchi type I geometry is analyzed for both models with and without self-interaction potential of the vector field, and for two different equations of state of the cosmic matter. We discuss and conclude our results in Section IV.

II. BOPP-PODOLSKY TYPE VECTOR DARK ENERGY MODELS
We first start by briefly introducing the basic theoretical ideas of the Bopp-Podolsky type electrodynamics in its standard formulation in Minkowski geometry. Then, by adopting, as a starting point, the view that higher order derivatives of the Maxwell tensor may play a significant role in vector type models of dark energy, we introduce the gravitational action for such a theoretical model. The gravitational field equations as well as the equations of the vector field are derived from the action together with an equation representing the covariant conservation of the energy-momentum tensor.

A. The Bopp-Podolsky model of electrodynamics
The Lagrangian density from which Maxwell's equations can be obtained by the usual variational principle is [46] where we use a system of units with c = 1. In the above equation F µν is the Maxwell electromagnetic field tensor, A µ is the four-vector potential of the field while j µ denotes the electromagnetic current. This Lagrangian is a function of the field variables and of their first derivatives. There is no reason why we should restrict ourselves to only first derivatives in the action and it therefore seems natural to try a generalization of Eq. (4) of the form [43,44] The usual variational principle applied to this Lagrangian leads to the field equations where The simplest choice for L 0 , as proposed by Bopp and Podolsky [43,44], is where m is a new fundamental constant with mass dimension 1. Using Eq. (8) in Eq. (7) we obtain so that the field Eq. (6) becomes where U µ has the property ∂ µ U µ = 0. By defining we obtain If we impose the condition ∂ µ A µ = 0, Eq. (11) becomes By introducing the four-potential a µ = A µ + U µ , we obtain To summarize, an interesting result in the Bopp-Podolsky theory is that the electromagnetic field equations, the potentials and the fie1d strengths can be written as the difference, respectively, of the potentials and field strengths of two distinct fields These two fields are described by two sets of separate field equations, with the first set corresponding to the standard Maxwell equations, while the second set represents Proca type field equations for particles with mass m. This result also provides a physical interpretation of the new fundamental constant m. However, the mass term appears with a wrong sing in the equation of motion, signaling that the massive vector field is a ghost. In order to make it clear, let us rewrite the Lagrangian (8) in the Lorentz gauge with the result One can easily check that the above Lagrangian is identical to Now using the transformation A µ → A µ + B µ , one obtains It is now seen that the kinetic term of the massive vector field appears with a wrong sign, signaling that the massive vector field is a ghost. In order to make the theory healthy at the background level, one should make the massive ghost non-dynamical. In this paper, we will consider the cosmology of this model, so the Maximum energy scale of our theory is H 0 . By assuming that the ghost mass is larger than the energy scale of our theory, the ghost mass does not have any dynamics at length scales smaller than H −1 0 which is the desired range of doing cosmology. This means that the Kinetic energy of the ghost field is much less than its potential energy. Noting that the ghost mass squared is equal to m 2 , one can see that for values m ≫ H 0 , the ghost becomes non-dynamical. This is what we consider in what follows.

B. Bopp-Podolsky type vector dark energy models
In the following we assume that the vector dark energy can be described by a Bopp-Podolsky type model, with the action given by where A µ is the dark energy potential related to dark energy field strength by where A 2 = A µ A µ , j µ = ρu µ is the cosmological matter 4-vector current, m is a constant with dimension of mass and S m is the action for ordinary matter. We have also added to the gravitational action the self-interacting potential V A 2 of the vector field, and we have allowed for the possibility of a direct coupling between the matter current j µ and dark energy vector potential A µ , with the strength of the coupling described by the constant β.
Varying the action (22) with respect to dark energy potential A µ and the metric g µν we have where a prime indicates derivative with respect to the argument and At this point, a note about the variation of the term A α j α ≡ ρu α A α is in order. The variation of the energy momentum tensor can be written as (see Appendix A) while the variation of the four-velocity of the particle is Putting all these results in the variation of A µ j µ , one can see that ρ dependence vanishes from the metric equation of motion. In the following we will assume that the energy momentum tensor of ordinary matter is that of a perfect fluid In order to write the equation of the vector field in a form similar to the one in Bopp-Podolsky electrodynamics, we introduce a new auxiliary vector field U µ , defined as Then the vector field equation (23) reduces to two coupled differential equations for A µ and U µ as and respectively. The conservation of the energy momentum tensor is now obtained by taking the covariant divergence of the metric field equation. After some algebra, one finds where θ = ∇ µ u µ is the expansion parameter, a µ = u ν ∇ ν u µ is the acceleration and h µν = g µν + u µ u ν . By taking the covariant derivative of equation (23) we obtain C. The Newtonian limit In the following we consider the weak field limit of the Bopp-Podolsky type vector dark energy model for a static source, i.e., the Newtonian limit. In this case the only non-zero component of the energy-momentum tensor is T 00 = ρ and one may easily find that R = −2∇ 2 Φ, where Φ is the Newtonian potential which is related to the metric component through g 00 = −1 + 2Φ. Note that in this paper we are considering the vector field A µ as the dark energy sector of the universe which should be very small in the Newtonian limit of the theory. So, we consider A µ as a first order perturbed field, the same order as Φ.
The trace of the equation (24) can be reduced to where T is the trace of the matter energy momentum tensor. We have to keep only the first order terms in Φ and A µ . This implies that the terms in the metric equation (24) which are quadratic in A µ do not contribute to this limit. With these assumptions one obtains the generalized Poisson equation where O(1) means that we only keep terms which are linear in A µ . This implies that only the V (A 2 ) = const. affect the Poisson equation, which is exactly the cosmological constant. Note that the minus sign behind the energy density is because of our convention in defining the Newtonian potential in g 00 .
Let us now consider the vector field equation (23) in the Newtonian limit. In this limit, the covariant derivatives should be replaced by partial derivatives, since we have assumed that the vector field A µ is a small quantity. One can then show that the vector field equation reduces to which is exactly the original Bopp-Podolsky equation.

A. The Isotropic Cosmology
In this section we want to consider the cosmological implications of the theory. First, let us assume that the geometry of the Universe is described by the Fiedmann-Robertson-Walker metric. With this choice the possible form of the vector field A µ should have the form to preserve homogeneity and isotropy. However, with this choice the vector field strength tensor F µν and therefore the Bopp-Podolsky term vanishes in our theory. One can easily find that the Friedmann and the vector field equations in the absence of matter fields in this case can be written as The simplest possibility to satisfy the last equation is that the potential becomes constant. This is the standard de Sitter type theory, with constant Hubble parameter H = H 0 = 8πV /3. One can however drive a self accelerated expanding universe by choosing other forms for the potential V (A 2 ). In these cases the (0)-component of the vector field should be constant in order to satisfy the equation V ′ (A 2 0 ). For example, in the case that V (A 2 ) = αA 2 + βA 4 , one should have A 0 = α/2β.

B. Anisotropic Cosmology -Bianchi I model
In order to make the theory non-trivial, we should assume that the vector field has a spatial component. So, we will assume that the Universe can be described by the Bianchi-I type metric of the form and the vector field can then be written as One should note that we have assumed that the (0) component of the vector field is zero. This is because this component does not contribute to the strength tensor F µν . Also, we assume that the matter content of the Universe consists of a perfect cosmological fluid, with energy momentum tensor where ρ is the total matter density (dark plus baryonic), and p is the matter thermodynamic pressure.
For later convenience, we will define the directional Hubble factors H i , the mean Hubble factor H, the anisotropy parameter A, the shear scalarΣ 2 and the deceleration parameter q as [47] With the above definitions, one can easily see that the quantity ∇ µ (V ′ A µ ) vanishes. In this case, the time component of the conservation equation (28) leads to the usual conservation equation of the forṁ while the (x)-component of the conservation equation gives βρḂ = 0. We can then assume that B = constant, a condition which further implies that the Bopp-Podolsky term vanishes, or one should conclude that β = 0, i.e. no matter/vector field coupling. We will choose the second choice and in the following we will assume that β = 0 and then the conservation equation for the ordinary matter field hold. With these assumptions, one can see that equation (29) is satisfied identically.

The cosmological field equations
Let us introduce a new variable F , defined as With the above assumptions, only the (x) component of the vector field equation of motion becomes non-zero, which can be written as ...
The Friedmann equations can then be simplified to and where we have defined Σ = 2/ √ 3Σ.

The case of massless vector field
Let us now investigate the cosmological implications of the Bopp-Podolsky theory with a massless vector field. In this case the potential term V (A 2 ) vanishes. In order to simplify the mathematical formalism, let us introduce a set of dimensionless variables (τ, r, P, f, h, µ, σ), defined as where H 0 is the present day value of the Hubble function. By using the above set of variables, the cosmological evolution equations for the Bopp-Podolsky type vector dark energy model can be written as and and the energy conservation equation becomesṙ where now, "dot" represents derivative with respect to τ . One should note that because we want to make the ghost degree of freedom non-dynamical, one should assume µ ≫ 1. From Eq. (52) we can obtainf as After substituting this expression off into Eq. (53), we can solve Eqs. (53) and (54) to obtain the expressions oḟ h andσ, respectively. Therefore the system of equations describing the evolution of the anisotropic Bianchi type I Universe in the presence of Bopp-Podolsky type vector dark energy can be written aṡ After adopting an equation of state for the cosmological matter, the system of differential equations Eqs. (57)-(61) must be solved by choosing some appropriate initial conditions, which we take as f (0) = f 0 , u(0) = u 0 , h(0) = h 0 , σ(0) = σ 0 , and r(0) = r 0 , respectively. In the dimensionless variables introduced above the deceleration parameter is given by In the following we will assume that the age of the Universe is of the order of t max = (2/3)t H , which gives for the dimensionless time τ the maximum value of τ max = 2/3 = 0.66.
Eqs. (65)-(67) can be solved to give the matter energy density and the pressure as (70) By assuming 2h + σ = 0, Eq. (64) can be immediately integrated to give where we have denoted and we have used the initial condition h (τ 0 ) = h 0 and σ (τ 0 ) = σ 0 , respectively. Then the requirement of the equality of the pressures in Eqs. (69) and (70) gives for h the evolution equatioṅ Eq. (73) is a Riccati type equation, which generally cannot be solved exactly. For the matter energy density and pressure we obtain The solutions are periodic, with the period T = 2 √ 2π/ (µH 0 ) = 2 √ 2πt H /µ, where t H = 1/H 0 is the present day age of the Universe. Hence one period describes roughly the entire cosmological history. In the rescaled dimensionless time τ this corresponds to a time interval 2 √ 2π/µ. An approximate simple solutions of Eq. (73) can be obtained by assuming the conditions µ 2 = 8πf 2 0 , and µ (τ − τ 0 ) / √ 2 << δ. Then Eq. (73) takes the forṁ with the general solution where c 1 is an arbitrary constant of integration, and we have assumed (2h 0 + σ 0 ) 2 > 8µ 2 . In the limit of large times the mean Hubble function tends to a constant, lim τ →∞ h = (2h 0 + σ 0 ) 2 − 8µ 2 + 2h 0 + σ 0 /6 = constant, while the volume V of the Universe increases according to For the deceleration parameter we obtain In the large time limit lim τ →∞ q = −1, and thus the anisotropic Universe ends in a de Sitter exponentially accelerating phase. However, super-accelerating phases of evolution with q < −1 are also possible. In the same limit we obtain In the large time limit the matter energy density and pressure also reach some constant values, strongly dependent on the model parameters.

Cosmological evolution of the anisotropic radiation fluid Universe
The radiation epoch, in which the Universe consisted of a plasma of nuclei, electrons and photons, is one of the most important periods in the evolution of the Universe. During this period the temperature was in the range of 10 9 −10 3 K, and the temperatures remained too high for the binding of electrons to nuclei. Therefore during this phase the Universe was filled with a radiation fluid, described by the equation of state P = r/3. The radiation era lasted from around t = 10 s to t = 10 13 s [50], giving for the dimensionless time τ the range τ ∈ τ in = 2.18 × 10 −17 , τ f in = 2.18 × 10 −5 , where for the present day value of the Hubble constant we have adopted the numerical value H 0 = 67.31 km/Mpc s = 2.185 × 10 −18 s [48]. We approximate the initial value of the Hubble function at the beginning of the radiation era as being given by h (τ in ) ≈ 1/2τ in = 2.28 × 10 16 , with the initial dimensionless density of the matter in the Universe given as r (τ in ) ≈ 3h 2 (τ in ) = 1.57 × 10 33 . The initial value of the shear scalar can be obtained as σ (τ in ) = (3/2)A (τ in ) h 2 (τ in ). We define (arbitrarily) the initial value of the anisotropy parameter as A (τ in ) = 1, which gives σ (τ in ) = (3/2)h (τ in ) In order to obtain the evolution of the anisotropic Universe in the presence of a Bopp-Podolsky type dark energy, we have integrated numerically Eqs. (57)-(61) by using the following initial conditions: f (τ in ) = 0.55, u (τ in ) = −1, σ (τ in ) = 2.80 × 10 16 , h (τ in ) = 2.28 × 10 16 , and r (τ in ) = 1.57 × 10 33 , respectively. The time variations of the mean Hubble function, matter energy density, shear scalar, anisotropy parameter, deceleration parameter, and of the ratio of the time variation of the Bopp-Podolsky vector potential and scale factor, respectively, are presented in Figs. 1-3, for different values of the field mass term µ.
As one can see from Fig. 1, the mean Hubble function is a monotonically decreasing function of time, indicating an expansionary evolution of the anisotropic Bianchi type I cosmological model. The cosmological evolution rate is practically independent on the numerical values of µ. In the large time limit very small differences in the expansion rate, determined by the variation of µ, may appear. The radiation fluid energy density, shown in the right panel of Fig. 1, is also a monotonically decreasing function of the cosmological time, and its time evolution is not affected significantly by the variations of the numerical values of µ. The shear scalar σ, shown in the left panel of Fig. 2, decreases rapidly during the cosmological evolution, showing an almost linear dependence on τ . The evolution of the shear scalar is influenced by the numerical values of µ only in the large time limit. A similar dynamics can be seen for the time evolution of the anisotropy parameter A, which decreases significantly, indicating the tendency of the anisotropic Universe to evolve towards an isotropic stage. The evolution of A depends on the numerical values of µ only at the late phases of the radiation era. The mean deceleration parameter, presented in Fig. 3, has only positive values, indicating a decelerating expansion, which essentially depends on the numerical values of µ. The time variation of the Bopp-Podolsky vector potential f , depicted in the right panel of Fig. 3, shows a strongly µ-dependent dynamical evolution, with the function f monotonically decreasing in time.

Cosmological evolution of the anisotropic dust Universe
As a second application of the anisotropic Bopp-Podolsky type cosmological model we consider the evolution of the dust, matter dominated, Bianchi type I Universe, with P = 0. We assume that the matter dominated era began when the Universe was about 400,000 years old (after the recombination era), corresponding to an initial value of the dimensionless time coordinate of τ in ≈ 3 × 10 −5 . For the numerical value of the Hubble function at the beginning of the matter dominated era we adopt the value h (τ in ) = 2.5 × 10 4 , while for the initial value of the dimensionless energy density of the matter we assume the value r (τ in ) = 1.75 × 10 9 . We assume for the initial value of the anisotropy parameter at the beginning of the matter dominated era the value A (τ in ) = 0.60, giving for the initial value of the shear scalar σ (τ in ) = 2.31 × 10 4 .
The cosmological evolution is obtained by numerically integrating Eqs. As one can see from Fig. 4, in the case of the dust Bianchi type I Universe in the presence of a Bopp-Podolsky field, both the mean Hubble function and the matter energy-density are monotonically decreasing functions of time, indicating an expansionary cosmological dynamics. Both the rate of the expansion, as described by the Hubble function, as well as the matter energy density show almost no µ in the large time limit. The shear scalar σ, shown in the left panel of Fig. 5, is a monotonically decreasing function of τ during the entire considered period of the cosmological expansion. The time evolution of σ is strongly dependent on the numerical values of µ in the initial and middle stages of evolution. In the long time limit the shear scalar reaches the (approximately) zero value around the time interval τ ≈ 0.15, corresponding to an age of the Universe of the order of t ≈ 6.864 × 10 16 s, and to a redshift of z ≈ 2. The mean anisotropy parameter A, represented in the right panel of Fig. 5, shows a similar behavior, being a monotonically decreasing function for all times. In the long time limit the anisotropy parameter tends to a constant, very small value, which can be approximated as zero for time intervals longer than τ = 0.15. For the initial range of considered time intervals the behavior of the anisotropy parameter depends on the numerical values of µ, but in the The deceleration parameter q, depicted in the left panel of Fig. 6, shows that for the considered initial and parameter values the anisotropic dust Universe starts its evolution from a decelerating state with q ≈ 0.4 at τ = 3 × 10 −5 . Due to the presence of the vector field the deceleration parameter is a monotonically decreasing function of time, reaching the value q ≈ 0.33 at τ ≈ 0.15. The overall dynamics of q depends on the numerical values of the Bopp-Podolsky parameter. The evolution of f , depicted in the right panel of Fig. 6, shows a monotonically decrease in time, with the large time dynamics also essentially dependent on the numerical values of µ. In the large time limit f becomes a constant.
Hence an initially anisotropic Bianchi type I Universe fully isotropizes in the presence of a Bopp-Podolsky type vector field, and at z ≈ 2 the expansion becomes isotropic. However, at the end of the anisotropic era the vector field f still survives as a constant component in the composition of the Universe, and it may play the role of an effective cosmological constant that may trigger the accelerated expansion of the Universe when z < 0.5.
C. The effect of the vector field mass on the cosmological evolution In the following we will consider the effect of the mass term M on the cosmological expansion. In this case, the cosmological equations governing the evolution of the anisotropic Bianchi type I Universe takes the form and respectively, where a "dot" denotes the derivative with respect to the dimensionless time parameter τ , and we have defined a dimensionless parameter n as M = nH 0 . Also, we will assume that µ ≫ 1. Hence the system of Eqs. (81)-(84) can be written asḟ σ = 1 6µ 4 a 2 a 2 16πf 2 µ 2 ((4h − σ)(2h + σ) − 3(P + r)) + µ 4 + 8πu 2 + 64π 2 f 4 (2h + σ) 2 + 32πµ 2 f u(2h + σ) + respectively, while the energy conservation equation becomeṡ The system of Eqs. In the following we will restrict our investigation to the cosmological evolution of the pressureless dust Universe models for a Bopp-Podolsky type dark energy in the presence of a mass term. In the following we consider the cosmological evolution of a Bianchi type I Universe filled with pressureless dust, with P = 0, in the presence of a massive Bopp-Podolsky type vector field. In order to facilitate the comparison with the massless case we use the same initial conditions to numerically integrate the cosmological evolution equations (85)-(90), that is, we adopt as initial conditions f (τ in ) = 2.4, u (τ in ) = −0.001, σ (τ in ) = 2.31 × 10 4 , h((τ in ) = 2.5 × 10 4 , and r (τ in ) = 1.75 × 10 9 , respectively, and we consider the same values of µ. For B and a we assume the initial conditions B (τ in ) = 0.3 and a (τ in ) = 0.1. Moreover, we fix the value of n as n = 2.3. The time evolutions of the Hubble function, matter energy density, shear scalar, mean anisotropy parameter, deceleration parameter, and Bopp-Podolsky vector potential are presented in Figs. 7-9, respectively. the cosmological evolution on the Bopp=Podolsky parameter µ, at least for the considered range of µ. The Hubble function, shown in the left panel of Fig. 7, is a decreasing function of time, with the rate of expansion of the Universe independent on the numerical values of the model parameter µ in the large time limit. The energy density of the radiation fluid decreases monotonically during the expansionary phase, and its overall dynamics is not influenced significantly by the modifications of µ. The shear scalar σ, presented in the left panel of Fig. 8, is a monotonically decreasing, linear function of time, and in the large time limit its behavior is not influenced by the modifications of the numerical values of µ. For the considered range of parameters σ becomes (approximately) zero at τ ≈ 0.10. The anisotropy parameter A(τ ), shown in the right panel of Fig. 8, decreases rapidly to zero, indicating that the Universe isotropizes in a shorter time interval, as compared to the massless case. The evolution of A is overall independent on the numerical values of µ. The anisotropic dust Universe accelerates more rapidly as compared to the massless case, and enters the isotropic phase with q ≈ 0.10, with the deceleration parameter, depicted in the left panel of Fig. 9, showing a strong dependence on the numerical values of µ. The Bopp-Podolsky vector field potential B, represented in the right panel of Fig. 9, is a slowly decreasing function of time, also showing a strong dependence on the numerical values of µ.

IV. DISCUSSIONS AND FINAL REMARKS
In this paper we have investigated a specific vector-tensor type gravitational theory which is inspired by the Bopp-Podolsky electrodynamics. The Bopp-Podolsky theory represents an interesting generalization of classical electrodynamics which has been extensively investigated at the level of elementary particle theory. Its basic idea, the addition of quadratic terms in the Maxwell tensor derivatives to the action of the vector-tensor type gravitational models, may prove to be a useful extension of the standard vector-tensor models with minimal coupling between gravity and the vector field. In our model we have also considered the possibility of the existence of some self-interaction processes in dark energy which can be described by means of a potential term V A 2 . The self interaction potential is a function of the square of the four-potential of the vector type dark energy. In the present analysis we have restricted our investigations to the case of the linear potential only with V ∝ A 2 . A possible non-minimal coupling between the matter current and the four-potential of the vector field was also considered. As compared to the vector models based on standard Maxwell electrodynamics, the addition of the new terms in the action enriches significantly the theoretical framework, thus opening the possibility of a more general approach to vector type dark energy models.
We have considered in detail the cosmological implications of the Bopp-Podolsky type vector dark energy. In our study we have concentrated on the cosmological evolution of the Bianchi type I cosmological models. We have investigated numerically two distinct scenarios, corresponding to the absence and presence of a self-interacting potential for the vector field. In view of the possible applications for the description of the early Universe we have considered first the radiation fluid model. We have investigated it for the case of the massless Bopp-Podolsky vector fields, by varying the numerical values of the parameter µ. The presence of a Bopp-Podolsky type dark energy component could induce a complex dynamical behavior of the early Universe, with the cosmological evolution decelerating (q > 0). The anisotropy parameter A slowly decreases during the radiation phase. Interestingly, the evolutions of the Hubble function, of the matter energy density and of the shear scalar are not significantly affected by the modifications of the numerical values of the model parameter µ. The evolution is also strongly dependent on the initial conditions adopted for the Hubble function, energy density, shear scalar, and of the function f , and its derivatives.
In the case of the dust anisotropic Universe, the cosmological evolution can reach an isotropic phase in both massless and massive cases. The Universe fully isotropizes in both models, so that A = 0 after a finite time interval, corresponding to z ≈ 2. The presence of the massive vector field significantly speeds up the cosmological evolution towards an isotropic phase. The nature of the cosmological evolution strongly depends on the adopted range of the numerical values of the parameter µ in the massless case, but this dependence is diminished in the presence of the massive vector field. On the other hand, the Bopp-Podolsky vector field becomes a constant at the end of the anisotropic phase, and its presence could trigger an accelerated, de Sitter type evolution for z < 0.5.
However, it is important to point out that in the present model the late-time de-Sitter evolution is not an attractor of the system, and that the longtime evolution of the Universe, extending well below the present time, is oscillatory. As shown in Section III B 3, the period of the cosmological oscillations is of the order of T ≈ (8.88/µ)t H , which for µ = 10 is of the same order as the present age of the Universe. Once the Universe reaches this age, the direction of the expansion is reversed, and the Hubble function becomes an increasing function of the cosmological time, and the Universe experiences an overall cyclical behavior [51][52][53][54][55], consisting of a succession of expanding and collapsing phases. From a theoretical point of view in the present model the classically oscillating solution is obtained by adding a vector type dark energy with quadratic terms in the Maxwell tensor derivatives to the ordinary baryonic fluid.
As a result the equation of state parameter w satisfies in the long time limit the condition w ≤ −1, indicating the possibility of future super-accelerating stages of evolution of the Universe, before its recollapse towards a decelerating phase. One interesting question is if this type of model can avoid the singularity theorems. On the other hand we expect that quantum gravity eects, arising from a full theory of quantum gravity, would become important when the size of the Universe approaches the Planck scale. A detailed analysis of the singular behavior of the model would also require the consideration of other quantum effects that may appear due to the presence of quantum fluctuations. There is also the interesting possibility that after a number of oscillations, the Universe filled with vector type dark energy with quadratic terms in the Maxwell tensor could evolve to the bounce point through quantum tunneling, and then expand again. Hence we may tentatively assume that a bouncing Universe that avoids the classical singularity is also possible in the present model.
Recently, precise observations of the cosmic microwave background radiation have provided the possibility of testing the fundamental cosmological predictions of ination on the primordial uctuations, such as, for example, scale independence and Gaussianity [48]. One of the basic ideas of inflation, also supported by the cosmic no-hair conjecture, conjectures that ination makes classical anisotropy negligibly small. However, recently a number of astrophysical observations of the large scale structure of the Universe have raised some questions about the validity of the principles of homogeneity and isotropy [49]. Together with the recent Planck results, these observations point towards the possibility of the existence of an intrinsic large scale anisotropy in the Universe. By adopting this line of thought, in [56] it was proposed that the global geometry of the Universe may be described by the homogeneous locally rotationally symmetric (LRS) class of metrics. These interesting geometries induce a preferred direction in the cosmic sky, and a CMB that is isotropic at the level of the background. Such models can be tested with CMB and supernovae data by using the distortion of the luminosity distances generated by the anisotropic geometry. Hence, by taking into account the latest Planck data, possible existence of small large scale cosmological anisotropies cannot be ruled out completely. However, their physical origins are still unknown, with the most popular explanation being that anisotropy is due to deviations from isotropy of the primordial fluctuations [48]. But presently no clearly established physical mechanism that could generate such deviations is known [57]. We suggest that the existence of the Bopp-Podolsky vector dark energy may provide the physical mechanisms that could explain the generation of the anisotropies in the early Universe, as well as their survival at the present time.
Vector field models face the important problem of their stability. Several vector-tensor gravitational theories contain instabilities in the form of ghosts, or unstable growth of the linearized perturbations [58]. The presence of these instabilities is due to the longitudinal vector polarization modes that appear in the vector-tensor models. The existence of ghosts or tachyons during the early time evolution or in the small wavelength regime of the vector-tensor type cosmological models may indicate that the vacuum of these models is unstable. The problem of stability with respect to small linear perturbations and the hyperbolicity is also of fundamental importance for the Bopp-Podolsky type dark energy models with derivatives of the Maxwell tensor in the action. A full solution of this problem can be obtained by considering perturbations of the field equations of our model, and investigating their stability. We will consider this topic in a future publication.
In this paper we have introduced a vector-tensor type model representing an extension of vector type dark energy models. Further investigations of the corresponding cosmological models may provide us with some measure for discriminating between different evolutionary scenarios suggested by theoretical structure of the theory. Moreover, this model may contribute to a better understanding of some other fundamental processes like, for example, inflation and structure formation, which have played a fundamental role in the evolution of our Universe. so that n = n µ n ν g µν g . (A3) For the perfect fluid the density is only a function of n, i.e. ρ = ρ(n). Using thermodynamics relations one obtains δρ = ∂ρ ∂n s δn = 1 n (ρ + p)δn.