# Cosmology in non-local Bopp-Podolski electrodynamics

Zahra Haghani<EMAIL_ADDRESS>School of Physics, Damghan University,
Damghan, 41167-36716, Iran.

###### Abstract

Non-local extension of the Bopp-Podolski non-linear electrodynamics is
investigated. It is well-known that the theory in flat space time, reduces to
the Proca theory. However, it will be shown that in curved space time the
resulting theory will differ from the coupled Einstein-Proca system. This
theory admits de sitter solution. The cosmological perturbations on top of the
de Sitter space-time shows that the tensor and vector modes are healthy.
However there is a scalar mode in this model which behaves like the Pais-
Uhlenbeck oscillator. This shows that this theory contains an Ostrogradski
ghost in the scalar sector. Anisotropic cosmology of the model is also
investigated and we will show that the behavior of the universe at late time
depends strongly on the initial conditions. However, independent of the
parameters of the theory, the model predicts an isotropic universe at late
times.

##  I Introduction

The late-time accelerated expansion of our universe is demonstrated by
different observational data acce . The late-time acceleration occurs at very
low energy scale $\sim 10^{-4}eV$, so to investigate this phenomenon one can
consider a modification to GR as an effective field theory of gravity. There
are many attempts in the literature to describe the observational data. The
main idea is adding new degrees of freedom to general relativity (GR) in a
consistent way so that the resulting theory does not contain instabilities and
also satisfy current observational data. Adding new degrees of freedom can be
done in different approches. The first approach is to introduce some matter
sources with a large negative pressure. This branch is the well-known as dark
energy branch. There are many types of such sources in the literature such as
scalar fields, vector fields and etc DEmodels . The $\Lambda$CDM model is the
simplest and most successful model to describe current acceleration of the
universe and this model is in a good agreement with observational data.
However this model has some theoretical and phenomenological problems,
motivated cosmologists to consider more complicated cases containing dynamical
extra degrees of freedom lamp . The second approach is to modify the gravity
sector in Riemannian geometry at large distances such as $f(R)$, scalar-tensor
theories modg . Also, one can go further and make a theory with a non-minimal
couplings between matter and geometry such as $f(R,T)$ fRT ,
$f(R,T,R_{\mu\nu}T^{\mu\nu})$ FRTR , $f(R,L_{m})$, etc matt . Cosmological
implications of such theories are well-known and investigated extensively in
the literature. Other possibilities in this branch is to consider non-
Reimannian geometries such as Weyl-Dirac theory weyl , Cartan geometry cartan
and their generalizations genwc .

In this work we will investigate a special class of a vector-tensor theory of
gravity which contains nonlocal self-interactions of the vector field. The
application of vector field theories in cosmological behavior of the universe
is well-known and the consequences of such an extension is investigated vastly
both for early times vecinf  and also for the late time cosmology veclate .
More important possibilities for the vector field Lagrangian would include the
Maxwell theory describing a massless vector field which has U(1) symmetry and
the Proca theory for the massive vector field.

The kinetic term of the Maxwell and also Proca theories
$-1/4F_{\mu\nu}F^{\mu\nu}$ has a famous problem that it diverges at the point
charge. There are some efforts in the literature to avoid this problem. Born
and Infeld BI  have considered non-linear terms of the strength tensor of the
form

$\mathcal{L}\propto 1-\sqrt{\det\left(\eta_{\mu\nu}+\beta F_{\mu\nu}\right)},$

to avoid the point charge divergence of the electric charge. The Born-Infeld
electrodynamics reduces to the Maxwell’s theory in low energy limits. The
cosmological implications of this kind of matter field has been vastly
investigated BIcos .

The other possibility for avoiding the aforementioned divergence is to add
some higher derivative self-interaction terms to the theory. This was first
done by Bopp and Podolsky BP  where their Lagrangian contains

$\displaystyle\mathcal{L}\propto-\frac{1}{m^{2}}\big{(}\partial_{\alpha}F^{\mu\nu}\big{)}^{2},$
(1)

where $m$ is a constant with mass dimension 1. After a field redefinitions it
can be shown that the Bopp-Podolsky Lagrangian describes two independent
dynamical fields corresponding to one massless and one massive vector field.
Because of the higher derivative nature of the Bopp-Podolsky interaction term,
the theory contains Ostrogradski instability which shows itself as a massive
vector ghost. One should note that both generalizations of the Maxwell theory
introduced above, respects U(1) symmetry. Cosmological implications of the
Bopp-Podolsky theory and also its ghost instability on top of de Sitter
expanding background has been investigated in BPcos .

In an effective action the origin of nonlocal terms can be seen in the one-
loop approximation via the heat kernel method. Also one can consider the
nonlocal contributions to the quantum effective action within the covariant
perturbation expansion in the field strengths and spacetime curvatures
barvinsky . From an independent viewpoint, nonlocal terms in classical action
may come from integrating out one of the healthy degrees of freedom of the
theory. In these cases the resulting nonlocal action is free from pathologies.
Also there are special higher dimensional gravity theories, in which reduction
to four dimensions gives an action including nonlocal terms, e.g.
$\sqrt{-\Box}$ in DGP model DGP . Nonlocal terms has been used to modify IR
and UV limits of GR and nonlocal cosmology is extensively investigated in the
literature noncos .

Recently, nonlocal terms has been extensively studied in massive field
theories and specially massive gravity. In these cases, one can retain the
gauge invariance of the theory by introducing some nonlocal terms. Nonlocal
terms for massive spin two fields was suggested in the context of
degravitation degra . Also, nonlocal massive gravity was studied in Maggior .
In this paper, we will consider the effect of nonlocal terms in a vector field
theory.

It is well-known that the massive vector theory can be written in a form that
it is explicitly $U(1)$ symmetric. This can be done by the additional
dynamical degree of freedom which behaves correctly under the $U(1)$
transformation. To illustrate this point further, let us consider a dynamical
massive vector field in flat space-time described by the Proca action of the
form

$\displaystyle S=\int
d^{4}x\left(-\frac{1}{4}F_{\mu\nu}F^{\mu\nu}-\frac{1}{2}m^{2}A_{\mu}A^{\mu}\right),$
(2)

where $F_{\mu\nu}=2\partial_{[\mu}A_{\nu]}$ is the field strength of the Proca
field $A_{\mu}$ with dimension $M$ and $m^{2}$ is the mass squared. This
theory has 3 degrees of freedom on flat space and breaks explicitly the $U(1)$
symmetry of the Maxwell’s theory due to existence of the mass term. Now, let
us Stueckelberg transform the field $A_{\mu}$ as

$A_{\mu}\rightarrow A_{\mu}+\partial_{\mu}\phi,$

where $\phi$ is a scalar field. The resulting action becomes

$\displaystyle S=\int d^{4}x\bigg{(}$
$\displaystyle-\frac{1}{4}F_{\mu\nu}F^{\mu\nu}-\frac{1}{2}m^{2}A_{\mu}A^{\mu}$
$\displaystyle-m^{2}\partial_{\mu}\phi
A^{\mu}-\frac{1}{2}m^{2}\partial_{\mu}\phi\partial^{\mu}\phi\bigg{)},$ (3)

which has an explicit $U(1)$ symmetry if both $A_{\mu}$ and $\phi$ fields
transform as

$\displaystyle A_{\mu}$ $\displaystyle\rightarrow
A_{\mu}+\partial_{\mu}\zeta,$ $\displaystyle\phi$
$\displaystyle\rightarrow\phi-\zeta,$ (4)

where $\zeta$ is an arbitrary function of coordinates. One can proceed further
and write the action (I) in a form that it contains only the field $A_{\mu}$
and also has an explicit $U(1)$ symmetry. This can be done by obtaining the
scalar field $\phi$ in terms of the vector field $A_{\mu}$ from the scalar
equation of motion $\Box\phi=-\partial_{\mu}A^{\mu},$ with the result

$\displaystyle\phi=-\frac{1}{\Box}(\partial_{\mu}A^{\mu}),$ (5)

and substituting back to the action (I). The result is

$\displaystyle S=\int d^{4}x\bigg{(}$
$\displaystyle-\frac{1}{4}F_{\mu\nu}F^{\mu\nu}+\frac{m^{2}}{4}F_{\mu\nu}\frac{1}{\Box}F^{\mu\nu}\bigg{)},$
(6)

which we have performed integration by parts to simplify the action. One can
see that the resulting action has a $U(1)$ symmetry, but it becomes non-local.
As a result, the Proca theory in flat space is equivalent to non-local theory
(6).

The generalization of the Proca action has also attracted considerable
attention in the literature gp . In these theories the general kinetic term
for a massive vector field beyond Maxwell term is considered which has 3
degrees of freedom. Generalization to the case of non-abelian vector field is
also considered in nonpro .

The line of the paper is as follows. In the next section, we will introduce
the model and obtain the field equations and also discuss on the local
counterpart of the theory. In section III we will investigate the isotropic
cosmology of the model and consider its cosmological perturbation on top of de
Sitter solution. Section IV will be devoted to the anisotropic implications of
the theory and in the last section we will conclude and discuss on some
possible issues and future lines.

##  II The model

In this paper, we will consider the effects of a Proca field in cosmology, but
we want to make the theory $U(1)$ symmetric. This can be achieved by using the
non-local counterpart of the Proca action. As a result, consider the following
action

$\displaystyle S=\int d^{4}x\sqrt{-g}$
$\displaystyle\bigg{(}\kappa^{2}(R-2\Lambda)$
$\displaystyle-\frac{1}{4}F_{\mu\nu}F^{\mu\nu}+\sum_{n}\alpha_{n}F_{\mu\nu}\frac{1}{\Box^{n}}F^{\mu\nu}\bigg{)},$
(7)

where we have generalized the non-local term to contain higher order
derivatives and $\Lambda$ is the cosmological constant. Here, $\alpha_{n}$ is
an arbitrary constant with dimension $M^{2n}$.

In this paper we will consider the case $n=1$ for the sake of simplicity. As
was discussed in Introduction, the above theory with $n=1$, can be considered
as a non-local version of the so-called BP electrodynamics BP ; BPcos . It is
well-known that the Einstein-BP theory suffers from Ostrogradski instability
in both vector and scalar modes because of its higher order derivative
occurring in the action farsi . However, one may expect that the non-local
version is free from instabilities (at least in the vector mode) because it
comes from the Proca theory. However, the procedure of obtaining (I) from (2),
highly depends on the geometry of space-time. In fact partial derivatives
should commute to make these two actions equal. In curved space-times the non-
local version of the Proca field (II) is different from the Einstein-Proca
system by some non-linear curvature terms which can turn on the Ostrogradski
instability of the theory. Another way to see this instability is to rewrite
the action (II) with $n=1$ as

$\displaystyle S=\int d^{4}x\sqrt{-g}$
$\displaystyle\bigg{(}\kappa^{2}(R-2\Lambda)-\frac{1}{4}F_{\mu\nu}F^{\mu\nu}$
$\displaystyle+\alpha
F^{\mu\nu}\xi_{\mu\nu}+\lambda^{\mu\nu}(\Box\xi_{\mu\nu}-F_{\mu\nu})\bigg{)},$
(8)

where $\xi_{\mu\nu}$ and $\lambda_{\mu\nu}$ are two covariant antisymmetric
tensors. Obtaining equations of motion of these fields and substituting back
to the action (II) one obtains the action (II). Here $\xi_{\mu\nu}$ is
dimensionless and $\lambda_{\mu\nu}$ has dimension $M^{2}$. One can decouple
the Maxwell field from the $\xi_{\mu\nu}$ and $\lambda_{\mu\nu}$ fields by
applying the field redefinition of the form

$\displaystyle F_{\mu\nu}\rightarrow
F_{\mu\nu}+2\alpha\xi_{\mu\nu}-2\lambda_{\mu\nu},$ (9)

which makes the action

$\displaystyle S$ $\displaystyle=\int
d^{4}x\sqrt{-g}\bigg{(}\kappa^{2}(R-2\Lambda)-\frac{1}{4}F_{\mu\nu}F^{\mu\nu}$
$\displaystyle+\alpha^{2}\xi_{\mu\nu}\xi^{\mu\nu}+\lambda_{\mu\nu}\lambda^{\mu\nu}-2\alpha\lambda_{\mu\nu}\xi^{\mu\nu}+\lambda_{\mu\nu}\Box\xi^{\mu\nu}\bigg{)}.$
(10)

The rest of the action containing $\xi_{\mu\nu}$ and $\lambda_{\mu\nu}$ can be
further diagonalized by performing the field transformation of the form

$\displaystyle\xi_{\mu\nu}$
$\displaystyle=-\frac{c}{\alpha}A_{\mu\nu}+B_{\mu\nu},$
$\displaystyle\lambda_{\mu\nu}$ $\displaystyle=cA_{\mu\nu}+\alpha B_{\mu\nu},$
(11)

with the result

$\displaystyle S$ $\displaystyle=\int
d^{4}x\sqrt{-g}\bigg{(}\kappa^{2}(R-2\Lambda)-\frac{1}{4}F_{\mu\nu}F^{\mu\nu}$
$\displaystyle+\alpha B_{\mu\nu}\Box
B^{\mu\nu}-\frac{c^{2}}{\alpha}(A_{\mu\nu}\Box A^{\mu\nu}-4\alpha
A_{\mu\nu}A^{\mu\nu})\bigg{)}.$ (12)

One can see from the above action that we have a massless antisymmetric tensor
$B_{\mu\nu}$ and a massive one $A_{\mu\nu}$. One should assume that $\alpha>0$
in order to make the field $B_{\mu\nu}$ healthy. In this case however, one
should assume $c^{2}<0$ to make $A_{\mu\nu}$ healthy. This in fact makes the
field redefinition (II) complex. The theory then consists of a massless
Maxwell field and a massless antisymmetric field $B_{\mu\nu}$ together with a
massive complex antisymmetric tensor field $A_{\mu\nu}$. We will not follow
this line in this paper, and take all of the fields real. We will show in the
following that an instability will be appeared in the scalar sector in the
perturbation analysis of the theory on top of the de Sitter solution.

Let us now obtain the equations of motion of the theory (II). The equation of
motion for the vector field $A_{\mu}$, the tensor field $\xi_{\mu\nu}$ and the
Lagrange multiplier $\lambda_{\mu\nu}$ can be written respectively as

$\displaystyle\nabla^{\beta}(F_{\alpha\beta}-2\alpha\xi_{\alpha\beta}+2\lambda_{\alpha\beta})=0,$
(13) $\displaystyle\alpha F^{\alpha\beta}+\Box\lambda^{\alpha\beta}=0,$ (14)

and

$\displaystyle\Box\xi_{\alpha\beta}=F_{\alpha\beta}.$ (15)

As one can see from the above field equations, the Maxwell strength tensor
$F_{\mu\nu}$ is the source for both antisymmetric tensors $\xi_{\mu\nu}$ and
$\lambda_{\mu\nu}$. Also it is worth mentioning that the theory has a $U(1)$
symmetry by shifting the vector field $A_{\mu}\rightarrow
A_{\mu}+\partial_{\mu}\varphi$. This implies that the $A_{\mu}$ field equation
satisfies a conservation equation of the form (13). As, a result, we have a
Noether current associated with the aforementioned symmetry.

The gravitational field equation of motion can be written as

$\displaystyle\kappa^{2}(G_{\mu\nu}+\Lambda\
g_{\mu\nu})-\frac{1}{2}F_{\mu}{}^{\alpha}F_{\nu\alpha}+\frac{1}{8}F_{\alpha\beta}F^{\alpha\beta}g_{\mu\nu}-\tfrac{1}{2}F^{\alpha\beta}g_{\mu\nu}\lambda_{\alpha\beta}-\tfrac{1}{2}\alpha
F^{\alpha\beta}\ g_{\mu\nu}\xi_{\alpha\beta}+\alpha
F_{\nu}{}^{\alpha}\xi_{\mu\alpha}$ $\displaystyle+\alpha
F_{\mu}{}^{\alpha}\xi_{\nu\alpha}-\tfrac{1}{2}\
\xi_{\nu}{}^{\alpha}\Box\lambda_{\mu\alpha}\
-\tfrac{1}{2}\xi_{\mu}{}^{\alpha}\Box\lambda_{\nu\alpha}+\tfrac{1}{2}\nabla_{\alpha}\lambda^{\alpha\beta}\nabla_{\nu}\xi_{\mu\beta}+\tfrac{1}{2}\lambda_{\nu}{}^{\alpha}\Box\xi_{\mu\alpha}+\tfrac{1}{2}\lambda_{\mu}{}^{\alpha}\Box\xi_{\
\nu\alpha}+\tfrac{1}{2}\xi_{\nu}{}^{\alpha}\nabla_{\beta}\nabla_{\mu}\lambda_{\alpha}{}^{\beta}$
$\displaystyle-\tfrac{1}{2}\xi^{\alpha\beta}\nabla_{\beta}\nabla_{\mu}\lambda_{\nu\alpha}+\tfrac{1}{2}\
\lambda_{\nu}{}^{\alpha}\nabla_{\beta}\nabla_{\mu}\xi_{\alpha}{}^{\beta}-\tfrac{1}{2}\lambda^{\alpha\
\beta}\nabla_{\beta}\nabla_{\mu}\xi_{\nu\alpha}+\tfrac{1}{2}\
\xi_{\mu}{}^{\alpha}\nabla_{\beta}\nabla_{\nu}\lambda_{\alpha}{}^{\beta}-\tfrac{1}{2}\xi^{\alpha\beta}\nabla_{\beta}\nabla_{\nu}\lambda_{\mu\alpha}+\tfrac{1}{2}\lambda_{\
\mu}{}^{\alpha}\nabla_{\beta}\nabla_{\nu}\xi_{\alpha}{}^{\beta}$
$\displaystyle-\tfrac{1}{2}\lambda^{\alpha\beta}\nabla_{\beta}\nabla_{\nu}\xi_{\mu\
\alpha}-\tfrac{1}{2}g_{\mu\nu}\nabla_{\gamma}\xi_{\alpha\beta}\
\nabla^{\gamma}\lambda^{\alpha\beta}+\tfrac{1}{2}\nabla_{\beta}\xi_{\nu\alpha}\nabla_{\mu}\lambda^{\alpha\beta}-\tfrac{1}{2}\nabla_{\beta}\xi_{\alpha}{}^{\beta}\nabla_{\mu}\lambda_{\
\nu}{}^{\alpha}+\tfrac{1}{2}\nabla^{\beta}\lambda_{\nu}{}^{\alpha}\
\nabla_{\mu}\xi_{\alpha\beta}+\tfrac{1}{2}\nabla_{\alpha}\lambda^{\alpha\beta}\nabla_{\mu}\xi_{\nu\beta}$
$\displaystyle+\tfrac{1}{2}\nabla_{\beta}\xi_{\mu\alpha}\nabla_{\nu}\lambda^{\alpha\beta}+\tfrac{1}{2}\nabla_{\mu}\xi_{\alpha\beta}\nabla_{\nu}\lambda^{\alpha\beta}-\tfrac{1}{2}\nabla_{\beta}\xi_{\alpha}{}^{\beta}\nabla_{\nu}\lambda_{\mu}{}^{\alpha}+\tfrac{1}{2}\nabla^{\beta}\lambda_{\mu}{}^{\alpha}\nabla_{\nu}\xi_{\alpha\beta}+\
\tfrac{1}{2}\nabla_{\mu}\lambda^{\alpha\beta}\nabla_{\nu}\xi_{\alpha\beta}=0.$
(16)

In the following the cosmological solution of this model is considered.

##  III Isotropic cosmology

Let us now assume that the universe can be described by a flat FRW line
element of the form

$\displaystyle ds^{2}=-dt^{2}+a^{2}(dx^{2}+dy^{2}+dz^{2}),$ (17)

where $a$ is the scale factor. The homogeneity and isotropy of the FRW metric
implies that the vector field $A_{\mu}$ should have a form

$\displaystyle A_{\mu}=(A_{0}(t),0,0,0).$ (18)

This form is the only choice which preserve isotropy and homogeneity of the
universe. The antisymmetric tensor fields $\lambda_{\mu\nu}$ and
$\xi_{\mu\nu}$ has six degrees of freedom consists of a 3-vector corresponding
to the (0i) components and a pseudo 3-vector field corresponding to (ij)
components maggoire . However, in an isotropic and homogeneous FRW universe
these tensors should vanish, because otherwise there is a non-zero 3-vector
field which breaks the isotropy of the space-time. So, in the FRW universe, we
have $\lambda_{\mu\nu}=0=\xi_{\mu\nu}$. Also note that the form (18) for the
vector field implies that the field strength $F_{\mu\nu}$ vanishes. As a
result, equations (13)-(15) are trivially satisfied and equation (II) reduces
to the Einstein’s equation in the presence of the cosmological constant which
has a de Sitter solution with the Hubble parameter of the form

$\displaystyle H=\sqrt{\frac{\Lambda}{3}}.$ (19)

Note that in an isotropic universe, the non-local BP interaction does not
contribute to the evolution of the universe. In the following, we will see
that the non-local BP affects the perturbations around de Sitter space time.
Also, in the next section, we will consider the anisotropic universe in which
the BP interaction term will affect the evolution of the universe.

###  III.1 Cosmological perturbations

Let us consider the perturbation of the action (II) on top of the de Sitter
solution obtained in the previous section. The metric perturbation can be
written as

$\displaystyle ds^{2}$
$\displaystyle=-(1+2\varphi)\,dt^{2}+2a(S_{i}+\partial_{i}B)dx^{i}\,dt$
$\displaystyle+a^{2}\big{(}(1+2\psi)\delta_{ij}+\partial_{i}\partial_{j}E+\partial_{(i}F_{j)}+h_{ij}\big{)}dx^{i}dx^{j},$
(20)

where $h_{ij}$ is the traceless and transverse tensor mode with
$h_{ii}=0=\partial_{i}h_{ij}$, $F_{i}$ and $S_{i}$ are transverse vector modes
with $\partial_{i}F_{i}=0=\partial_{i}S_{i}$, and $\psi$, $\varphi$, $B$ and
$E$ are four scalar modes. Note that the special indices $(i,j=1,2,3)$ are
raised and lowered with $\delta_{ij}$.

The vector field $A_{\mu}$ can be decomposed as

$\displaystyle A_{\mu}=(A_{0}+\delta A_{0},A^{\perp}_{i}+\partial_{i}\delta
A),$ (21)

where $\delta A_{0}$ and $\delta A$ are scalar modes and $A^{\perp}_{i}$ is
the vector mode with $\partial_{i}A^{\perp}_{i}=0$. Note that the background
value $A_{0}$ is an arbitrary function of time since it is not contribute to
the background cosmology. However in this section we will assume that $A_{0}$
is constant for the sake of simplicity. For the antisymmetric tensor fields
$\xi_{\mu\nu}$ and $\lambda_{\mu\nu}$, the background values are zero and one
has a decomposition

$\displaystyle\xi_{i0}$ $\displaystyle=\xi^{\perp}_{i}+\partial_{i}\xi,$
$\displaystyle\xi_{ij}$
$\displaystyle=\epsilon_{ij}^{~{}~{}k}(\chi^{\perp}_{k}+\partial_{k}\chi),$
(22)

and

$\displaystyle\lambda_{i0}$
$\displaystyle=\lambda^{\perp}_{i}+\partial_{i}\lambda,$
$\displaystyle\lambda_{ij}$
$\displaystyle=\epsilon_{ij}^{~{}~{}k}(\rho^{\perp}_{k}+\partial_{k}\rho),$
(23)

where $\xi^{\perp}_{i}$, $\chi^{\perp}_{i}$, $\lambda^{\perp}_{i}$ and
$\rho^{\perp}_{i}$ are vector modes with vanishing divergence, and $\xi$,
$\chi$, $\lambda$ and $\rho$ are scalar modes. Also, $\epsilon_{ij}^{~{}~{}k}$
is the 3-dimensional Levi-Civita symbol. In summary, we have one tensor mode
associated to the metric perturbation, 7 vector modes and 10 scalar modes.
Now, let us define the gauge invariant quantities associated to above
perturbation variables. For the metric perturbation, one can define two gauge
invariant scalar perturbations of the form

$\displaystyle\Phi$
$\displaystyle=\varphi+\partial_{t}\left(aB-\frac{a^{2}}{2}\partial_{t}E\right),$
$\displaystyle\Psi$
$\displaystyle=\psi+H\left(aB-\frac{a^{2}}{2}\partial_{t}E\right),$ (24)

and one gauge invariant vector perturbation

$\displaystyle\Pi_{i}=S_{i}-\frac{1}{2}a\partial_{t}F_{i}.$ (25)

Also, the tensor perturbation $h_{ij}$ is gauge invariant. For the vector
field $A_{\mu}$, one can define two gauge invariant scalar perturbations of
the form

$\displaystyle\mathcal{Y}$ $\displaystyle=\delta
A_{0}+A_{0}\partial_{t}\left(aB-\frac{a^{2}}{2}\partial_{t}E\right),$
$\displaystyle\mathcal{Z}$ $\displaystyle=\delta
A+A_{0}\left(aB-\frac{a^{2}}{2}\partial_{t}E\right).$ (26)

The vector mode $A^{\perp}_{i}$ is gauge invariant. Finally, one can easily
check that because of the zero background values for the antisymmetric tensor
fields $\xi_{\mu\nu}$ and $\lambda_{\mu\nu}$, all vector and scalar modes
associated to these fields are gauge invariant.

After expanding the action (II) up to second order in perturbation parameters
around the de Sitter background (19), one can see that tensor, vector and
scalar modes are decoupled from each other and as a result, in the following
we will consider these modes independently.

####  III.1.1 Tensor perturbations

Expanding the action (II) up to second order in the tensor mode $h_{ij}$, one
obtains

$\displaystyle
S^{(2)}_{tensor}=\frac{1}{2}\sum_{+,\times}\int\,d^{3}k\,dt\,\kappa^{2}\,a^{3}\left[\dot{h}^{2}_{ij}-\frac{\vec{k}^{2}}{a^{2}}h^{2}_{ij}\right],$
(27)

where we have performed Fourier transformation and sum over two helicity
degrees of freedom.

The above result shows that the tensor mode in Non-local BP theory is
equivalent to the standard Einstein’s theory with a cosmological constant.
This is however not surprising since we do not have a source for a tensor mode
from the non-linear BP interaction. As a result the speed of the propagation
of tensor mode in this theory is equal to the speed of light, satisfying
recent gravitational wave observations gw .

####  III.1.2 Vector perturbation

Expanding the action (II) up to second order in gauge invariant perturbation
vector modes and performing Fourier transformation, one obtains

$\displaystyle S^{(2)}_{vector}$ $\displaystyle=\int
d^{3}kdta\Bigg{[}\dot{\vec{A}}^{\perp
2}-\frac{\vec{k}^{2}}{a^{2}}\vec{A}^{\perp
2}+\kappa^{2}\vec{k}^{2}\vec{\Pi}^{2}+4\left(\frac{\vec{k}^{2}}{a^{2}}-\frac{4\Lambda}{3}\right)(\vec{\lambda}^{\perp}.\vec{\xi}^{\perp}-\vec{\rho}^{\perp}.\vec{\chi}^{\perp})+\sqrt{\frac{16\Lambda}{3}}(\vec{A}^{\perp}.\vec{\xi}^{\perp}+\dot{\vec{\xi}}^{\perp}.\vec{\lambda}^{\perp}-\dot{\vec{\chi}}^{\perp}.\vec{\rho}^{\perp})$
$\displaystyle+4(\dot{\vec{A}}^{\perp}.\vec{\lambda}^{\perp}+\ddot{\vec{\xi}}^{\perp}.\vec{\lambda}^{\perp}-\ddot{\vec{\chi}}^{\perp}.\vec{\rho}^{\perp})+4\alpha\vec{A}^{\perp}.\dot{\vec{\xi}}^{\perp}+\frac{4}{a}(\vec{\rho}^{\perp}-\alpha\vec{\chi}^{\perp}).(\vec{k}\times\vec{A}^{\perp})+\frac{8H}{a^{2}}\Big{(}\vec{\rho}^{\perp}.(\vec{k}\times\vec{\xi}^{\perp})+\vec{\lambda}^{\perp}.(\vec{k}\times\vec{\chi}^{\perp})\Big{)}\Bigg{]}.$
(28)

One can see from the above action that $\vec{\Pi}$, $\vec{\lambda}^{\perp}$
and $\vec{\rho}^{\perp}$ are non-dynamical with the equations of motion
$\vec{\Pi}=0$ and

$\displaystyle-3\dot{\vec{A}}-3\ddot{\vec{\xi}}-\sqrt{3\Lambda}\dot{\vec{\xi}}+4\Lambda\vec{\xi}-3\frac{\vec{k}^{2}}{a^{2}}\vec{\xi}-\frac{6H}{a}\vec{k}\times\vec{\chi}=0,$
$\displaystyle
3\ddot{\vec{\chi}}+\sqrt{3\Lambda}\dot{\vec{\chi}}-4\Lambda\vec{\chi}+3\frac{\vec{k}^{2}}{a^{2}}\vec{\chi}-\frac{6H}{a}\vec{k}\times\vec{\xi}-\frac{3}{a}\vec{k}\times\vec{A}=0,$
(29)

where we have dropped the $\perp$ index for simplicity. It should be noted
that these equations are in fact constraint equations due to the Lagrange
multiplier $\lambda_{\mu\nu}$. Solving the above equations, will give
$\vec{\xi}$ and $\vec{\chi}$ in terms of $\vec{A}$. As a result, we will have
one dynamical vector perturbation $\vec{A}$ with Lagrangian

$\displaystyle S^{(2)}_{vector}$ $\displaystyle=\int
d^{3}k\,dt\,a\Bigg{[}\dot{\vec{A}}^{\perp
2}-\frac{\vec{k}^{2}}{a^{2}}\vec{A}^{\perp
2}+\sqrt{\frac{16\Lambda}{3}}\alpha\vec{A}^{\perp}.\vec{\xi}^{\perp}$
$\displaystyle+4\alpha\vec{A}^{\perp}.\dot{\vec{\xi}}^{\perp}-\frac{4\alpha}{a}\vec{\chi}^{\perp}.(\vec{k}\times\vec{A}^{\perp})\Bigg{]},$
(30)

where $\vec{\xi}^{\perp}$ and $\vec{\chi}^{\perp}$ are known functions of
$\vec{A}^{\perp}$ from equations (III.1.2).

We should note that in the standard BP electrodynamics, there are two vector
modes on top of the de Sitter solution, one of which is an Ostrogradski ghost.
By promoting the BP interaction to a non-local term, the Ostrogradski ghost
becomes a Lagrange multiplier and vanishes from the action.

####  III.1.3 Scalar perturbations

Let us now consider the scalar perturbation of the theory (II). As is obtained
from the previous section, the theory has 8 gauge invariant scalar degrees of
freedom. After Fourier transforming the second order action in scalar modes,
one obtains

$\displaystyle S^{(2)}_{scalar}$
$\displaystyle=\int\,d^{3}k\,dt\,\Bigg{[}\frac{4\mathit{a}\lambda
k^{2}\sqrt{\Lambda}\dot{\xi}}{\sqrt{3}}+4\mathit{a}\lambda
k^{2}\ddot{\xi}-\frac{4\mathit{a}k^{2}\sqrt{\Lambda}\rho\dot{\chi}}{\sqrt{3}}-4\mathit{a}k^{2}\rho\ddot{\chi}-2\mathit{a}k^{2}(2\alpha\xi-2\lambda+\mathcal{Y})\dot{\mathcal{Z}}+\mathit{a}k^{2}\dot{\mathcal{Z}{}}^{2}$
$\displaystyle+8\sqrt{3}\mathit{a}^{3}\kappa^{2}\sqrt{\Lambda}\Phi\dot{\Psi}-12\mathit{a}^{3}\kappa^{2}\dot{\Psi}{}^{2}-4\mathit{a}^{3}\kappa^{2}\Lambda\Phi^{2}+\frac{4\lambda
k^{4}\xi}{\mathit{a}}-\frac{4k^{4}\rho\chi}{\mathit{a}}+8\mathit{a}\kappa^{2}k^{2}\Phi\Psi+4\mathit{a}\kappa^{2}k^{2}\Psi^{2}$
$\displaystyle-\frac{16}{3}\mathit{a}\lambda
k^{2}\Lambda\xi+\frac{16}{3}\mathit{a}k^{2}\Lambda\rho\chi+\mathit{a}k^{2}\mathcal{Y}^{2}+4\mathit{a}\alpha
k^{2}\xi\mathcal{Y}-4\mathit{a}\lambda k^{2}\mathcal{Y}\Bigg{]}.$ (31)

One can see from the above action that the scalar perturbations $\Phi$,
$\mathcal{Y}$, $\lambda$ and $\rho$ are non-dynamical, with equations of
motion

$\displaystyle\Lambda\Phi=\frac{k^{2}}{a^{2}}\Psi+\sqrt{3\Lambda}\dot{\Psi},$
$\displaystyle\mathcal{Y}=2\lambda+2\alpha\xi+\dot{\mathcal{Z}},$
$\displaystyle
2\lambda=\ddot{\xi}+\sqrt{\frac{\Lambda}{3}}\dot{\xi}-\left(\frac{k^{2}}{a^{2}}-4\Lambda-2\alpha\right)\xi,$
(32)

and

$\displaystyle\left(3\frac{k^{2}}{a^{2}}-4\Lambda\right)\chi+\sqrt{3\Lambda}\dot{\chi}+3\ddot{\chi}=0.$
(33)

It is worth mentioning that the scalar mode $\rho$ is a Lagrange multiplier
and the corresponding equation of motion determine the evolution of the
perturbed field $\chi$. Also the scalar mode $\lambda$ is an auxiliary field.

Substituting the above solutions back to the action (III.1.3), one can see
that the scalar mode $\Psi$ becomes non-dynamical with equation of motion
$\Psi=0$. This is similar to the Einstein-Hilbert theory since there is no
non-minimal interactions between the curvature tensor and the fields
$\xi_{\mu\nu}$ and $A_{\mu}$. The resulting action becomes

$\displaystyle S^{(2)}_{scalar}$
$\displaystyle=\int\,d^{3}k\,dt\,\frac{ak^{2}}{\alpha^{2}}\,\Bigg{[}\ddot{\xi}^{2}-\left(2\frac{k^{2}}{a^{2}}-\frac{4}{3}(2\Lambda-3\alpha)\right)\dot{\xi}^{2}$
$\displaystyle+\left(\frac{k^{4}}{a^{4}}+2(2\alpha-\Lambda)\frac{k^{2}}{a^{2}}+\frac{16}{9}\Lambda(\Lambda-3\alpha)\right)\xi^{2}\Bigg{]}.$
(34)

The above theory is equivalent to the Pais-Uhlenbeck oscillator PU  with
angular momentums

$\displaystyle\omega_{1,2}^{2}=\frac{k^{2}}{a^{2}}+2\alpha-\frac{4\Lambda}{3}\pm\sqrt{4\alpha^{2}-\frac{2\Lambda}{3}\frac{k^{2}}{a^{2}}},$
(35)

which is well-known to suffer from an Ostrogradski instability. This shows
that the scalar sector of the theory has a ghost degree of freedom. This is in
fact similar to the standard Bopp-Podolski electrodynamics farsi , so the non-
local theory can not remove the ghost in the scalar sector. In order to cure
this problem one can add a constraint term to reduce the dimensionality of the
phase space and then remove the Ostrogradski instability from the theory tolle
. It may be done by adding another higher derivative Maxwell term to the
action as a Lagrange multiplier which will be the scope of our future works.
From another point of view, the acion (II) with $n=1$ in Minkowski space time,
is the non local version of the Proca action. So it does not have any
pathology. However, as mentioned above it is not the case for curved space
time.

It should be noted that the $\xi$ mode corresponds to the $A_{\mu\nu}$ field
of section II. This can be seen in the subhorizon limit of the theory where
the above action reduces to

$\displaystyle S^{(2)}_{scalar}$
$\displaystyle=\int\,d^{3}k\,dt\,2a^{3}\,\left(\frac{k^{2}}{a^{2}}\right)^{2}\Bigg{[}\dot{\tilde{\xi}}^{2}-\frac{1}{2}\tilde{\xi}^{2}\Bigg{]},$
(36)

where $\tilde{\xi}=i\xi$. So, the $\tilde{\xi}$ field corresponds to the
scalar part of the complex tensor field $A_{\mu\nu}$ as was obtained in
section II.

##  IV Anisotropic cosmology

Let us now consider the anisotropic cosmology of the non-local BP
electrodynamics. Assuming that the universe can be described by the Bianchi
type-I space-time of the form

$\displaystyle ds^{2}=-dt^{2}+a^{2}dx^{2}+b^{2}(dy^{2}+dz^{2}),$ (37)

where $a(t)$ and $b(t)$ are the scale factors corresponding to spacial
directions $x$, $y$ and $z$. Let us define the following quantities which are
very useful in the analysis of anisotropic cosmology

$\displaystyle H_{1}$ $\displaystyle=\frac{\dot{a}}{a},\quad
H_{2}=\frac{\dot{b}}{b},\quad\Delta H_{i}=H_{i}-H,$ $\displaystyle 3H$
$\displaystyle=\sum_{i}H_{i}=H_{1}+2H_{2},$ $\displaystyle\quad 3A$
$\displaystyle=\sum_{i}\left(\frac{\Delta H_{i}}{H}\right)^{2},\quad
q=-1+\frac{d}{dt}\left(\frac{1}{H}\right).$ (38)

Here, $H$ is the mean Hubble parameter, $A$ is the anisotropic parameter and
$q$ is the deceleration parameter.

In order to solve the set of equations (13)-(II), let us perform some
assumptions. From equation (13), one obtains

$\displaystyle\nabla^{\beta}(F_{\alpha\beta}-2\alpha\xi_{\alpha\beta}+2\lambda_{\alpha\beta})=0,$
(39)

which has a solution

$\displaystyle\lambda_{\alpha\beta}=\alpha\xi_{\alpha\beta}-\frac{1}{2}F_{\alpha\beta},$
(40)

where we have set the integration constant to zero. Note that $F_{\mu\nu}$ is
the field strength of the vector field $A_{\mu}$. Inspiring from solution
(40), one can assume that $\xi_{\alpha\beta}$ and $\lambda_{\alpha\beta}$ are
also the field strength of some vector fields $\xi_{\mu}$ and $\lambda_{\mu}$
respectively

$\displaystyle\xi_{\alpha\beta}=\nabla_{\alpha}\xi_{\beta}-\nabla_{\beta}\xi_{\alpha},$
$\displaystyle\lambda_{\alpha\beta}=\nabla_{\alpha}\lambda_{\beta}-\nabla_{\beta}\lambda_{\alpha}.$
(41)

Substituting back to (40), one obtains

$\displaystyle\lambda_{\mu}=\alpha\,\xi_{\mu}-\frac{1}{2}A_{\mu}.$ (42)

Now, let us assume that the vectors $A_{\mu}$ and $\xi_{\mu}$ can be
decomposed as

$\displaystyle A_{\mu}$ $\displaystyle=\kappa\left(0,\int a\,U_{1}dt,\int
b\,U_{2}dt,\int b\,U_{2}dt\right),$ $\displaystyle\xi_{\mu}$
$\displaystyle=\frac{\kappa}{H_{0}^{2}}\left(0,\int a\,X_{1}dt,\int
b\,X_{2}dt,\int b\,X_{2}dt\right),$ (43)

where $H_{0}$ is the current Hubble parameter and $U_{i}$ and $X_{i}$ with
$i=1,2$ are some arbitrary dimensionless functions of time. The time component
of the above vectors are set to zero, because they do not contribute to the
field equations. Also, the complicated form of the above assumption will help
us to write the field equations in a more compact form.

Now, let us define the following dimensionless parameters

$\displaystyle\tau=H_{0}t,\quad H=H_{0}h,\quad\alpha=\beta
H_{0}^{2},\quad\Lambda=H_{0}^{2}\lambda.$ (44)

With the above assumptions, equations of motion (13)-(II) reduces to

$\displaystyle
h^{2}\left[3(\tilde{A}^{2}-1)+2(\tilde{A}^{2}+1)S_{1}X_{1}+2(2-2\tilde{A}+5\tilde{A}^{2})S_{2}X_{2}\right]+\lambda+\sum_{j=1}^{2}\bigg{[}j\left(S_{j}^{\prime}X_{j}^{\prime}-X_{j}S_{j}^{\prime\prime}+S_{j}X_{j}^{\prime\prime}-\frac{1}{4}U_{j}^{2}\right)\bigg{]}=0,$
(45) $\displaystyle h^{\prime}$
$\displaystyle\left[-2(\tilde{A}+1)+4(1-2\tilde{A})S_{2}X_{2}\right]-h^{2}\left[(\tilde{A}+1)^{2}(3+2S_{1}X_{1})+2(-4+10\tilde{A}+5\tilde{A}^{2})S_{2}X_{2}\right]$
$\displaystyle-h\left[2(1+4S_{2}X_{2})\tilde{A}^{\prime}+3X_{1}^{2}\left(\frac{S_{1}}{X_{1}}\right)^{\prime}+4(2\tilde{A}-1)\left(S_{2}X_{2}\right)^{\prime}\right]+\lambda+2U_{2}(S_{2}+\beta
X_{2})+S_{1}X_{1}^{\prime\prime}-X_{1}S_{1}^{\prime\prime}$
$\displaystyle\qquad-\sum_{j=1}^{2}\bigg{[}j\left(S_{j}^{\prime}X_{j}^{\prime}+\frac{1}{4}U_{j}^{2}\right)\bigg{]}=0,$
(46) $\displaystyle h^{2}$
$\displaystyle\left[2Z_{1}+2S_{1}X_{1}-3+\tilde{A}(1-\tilde{A})(2Z_{1}+8S_{2}X_{2}+3)\right]+h^{\prime}\left[2Z_{1}-2+\tilde{A}(2Z_{1}+1)\right]+h\bigg{[}\frac{1}{2}\tilde{A}^{\prime}\left(1+2Z_{1}\right)+2(1+\tilde{A})Z_{1}^{\prime}$
$\displaystyle-3X_{2}^{2}\left(\frac{S_{2}}{X_{2}}\right)^{\prime}\bigg{]}+\lambda-\frac{1}{4}U_{1}^{2}+2S_{2}U_{2}+U_{1}(S_{1}+\beta
X_{1})-S_{1}^{\prime}X_{1}^{\prime}-2S_{2}^{\prime}X_{2}^{\prime}-X_{2}S_{2}^{\prime\prime}+S_{2}X_{2}^{\prime\prime}=0,$
(47) $\displaystyle 2$
$\displaystyle(\tilde{A}-2)\left[h^{\prime}Z_{2}+3h^{2}\left(S_{2}X_{1}-S_{1}X_{2}+\tilde{A}Z_{2}\right)\right]+2h\left[\tilde{A}^{\prime}Z_{2}+\tilde{A}\left(4Z_{2}^{\prime}-6(S_{1}X_{2}^{\prime}+X_{1}S_{2}^{\prime})\right)+Z_{2}-6(X_{1}S_{2}^{\prime}+X_{2}S_{1}^{\prime})\right]$
$\displaystyle-4\beta U_{2}X_{1}+U_{1}(2U_{2}-4\beta
X_{2})+2(S_{1}X_{2}^{\prime\prime}+S_{2}X_{1}^{\prime\prime}-X_{2}S_{1}^{\prime\prime}-X_{1}S_{2}^{\prime\prime})=0,$
(48) $\displaystyle-2$
$\displaystyle(\tilde{A}+1)\left[3h^{2}+2h^{\prime}\right]S_{2}X_{2}+h\left[5X_{2}S_{2}^{\prime}+S_{2}\left(2\tilde{A}^{\prime}X_{2}-X_{2}^{\prime}\right)+2\tilde{A}\left(S_{2}X_{2}\right)^{\prime}\right]+\frac{1}{2}U_{2}\left(U_{2}-4\beta
X_{2}\right)$ $\displaystyle-
X_{2}S_{2}^{\prime\prime}+S_{2}X_{2}^{\prime\prime}=0,$ (49) $\displaystyle
U_{1}^{\prime\prime}+3hU_{1}^{\prime}+U_{1}\left(4\beta-2h^{2}(1+\tilde{A})^{2}\right)=0,$
(50)

and

$\displaystyle
U_{2}^{\prime\prime}+3hU_{2}^{\prime}+U_{2}\left(4\beta-h^{2}(2-2\tilde{A}+5\tilde{A}^{2})\right)=0,$
(51)

where prime represents derivative with respect to $\tau$ and we have defined
$\tilde{A}=\sqrt{A/2}$ and we have defined

$\displaystyle S_{j}$ $\displaystyle=\beta X_{j}-\frac{1}{2}U_{j},$
$\displaystyle Z_{1}$ $\displaystyle=S_{1}X_{1}+S_{2}X_{2},$ $\displaystyle
Z_{2}$ $\displaystyle=S_{1}X_{2}+S_{2}X_{1}.$

Figure 1: Variation of the deceleration parameter $q$ as a function of $\tau$
for different values of the parameters $\beta$ and $\lambda$: $\beta=1$ and
$\lambda=1$ (solid curve), $\beta=0.4$ and $\lambda=0.8$ (dotted curve),
$\beta=1.1$ and $\lambda=1.9$ (dot-dashed curve), respectively. The initial
conditions are set to $h(0)=1$ and $A(0)=2$.

In figures (1) and (2) we have shown the numerical solution of the above
equations for deceleration parameter $q$ and anisotropy parameter $A$ in terms
of the dimensionless time parameter $\tau$ for three different values of
$\beta=1.1,1,0.4$ and $\lambda=1.9,1,0.8$ respectively.

Figure 2: Variation of the anisotropy parameter $A$ as a function of $\tau$
for different values of the parameters $\beta$ and $\lambda$: $\beta=1$ and
$\lambda=1$ (solid curve), $\beta=0.4$ and $\lambda=0.8$ (dotted curve),
$\beta=1.1$ and $\lambda=1.9$ (dot-dashed curve), respectively.

These figures show that the late time behavior of the theory is consistent
with the late time dynamics of the universe; the dynamic of universe begins
from a highly anisotropic state and at the late time the anisotropy parameter
goes to zero which dictated that universe becomes isotropic. Also, the study
of the deceleration parameter (1), shows that the universe enters to the
accelerated expansion phase. This shows that the non-linear BP interaction
affects the IR behavior of the space-time making the universe to accelerate at
late times.

##  V Conclusions

In this paper we have considered the non-local modification of the non-linear
Bopp-Podolski electrodynamics. The Bopp-Podolski term is a higher derivative
electromagnetic self-interaction term which can solve the problem of the
divergence of the electric potential at the source point. Because of the
higher derivative nature of the interaction, the theory suffers from
Ostrogradski ghost instability. The cosmological implications of this higher
derivative self-interaction term was investigated extensively in BPcos  where
the authors shows that the vector and scalar modes suffer from instability on
top of de Sitter background. In the present paper we have considered the
possibility that whether the promotion of Bopp-Podolski interaction to a non-
local term can cure the theory from instability. The non-local version of the
Bopp-Podolski interaction term in flat space time is reduced to the co-called
Proca action, describing a massive spin-1 field. However, in curved space-
times, the action differs from Einstein-Proca theory. In fact one can
transform the action to a local theory with an additional Lagrange multiplier,
and the resulting action reduces to Einstein-Maxwell theory with two
additional antisymmetric tensor fields one massive and one massless.

The theory admits a de Sitter solution on the FRW background trivially, since
the vector field can not contribute to an isotropic and homogeneous universe.
Cosmological perturbations on top of this de Sitter solution, reveals that the
tensor mode is identical to $\Lambda$CDM model. This shows that the
propagation speed of the gravitational waves in this theory is equal to the
speed of light, satisfying recent gravitational wave observations gw . The
vector mode however, differs from Bopp-Podolski theory since there exists an
additional Lagrange multiplier which kills the ghost degree of freedom in the
vector mode. As a result we have shown that promoting to non-local theory,
makes the vector sector of the Bopp-Podolski theory healthy. The scalar sector
however has a mode which behaves like the Pais-Uhlenbeck oscillator. Because
of higher derivative nature of the Pais-Uhlenbeck Lagrangian, the scalar
sector suffers from Ostrogradski instability. This instability can be cured by
reducing the dimensionality of the phase space, as is done tolle . It is worth
mentioning that in proof  the authors has proven that there is no local
healthy theory with $U(1)$ symmetry. However, our theory is non-local and such
a theorem does not apply to our case. As a result such an interaction can in
principle be constructed which will be the scope of our future works.

As we have discussed above, the Bopp-Podolski interaction term does not
contribute to the homogeneous and isotropic FRW universe. So, we have
investigated the dynamics of the universe in anisotropic Bianchi-I universe.
We have seen that the universe starts from a highly anisotropic and
decelerating phase and at late times the universe becomes isotropic and
accelerating. As a result the model satisfies observational data on the late
time behavior of the universe.

## References

  * (1) A. G. Riess et al., Astron. J. 116 (1998) 1009, arXiv:astro-ph/9805201; S. Perlmutter et al., Astrophys. J. 517(1999) 565, arXiv:astro-ph/9812133; Planck Collaboration: P. A. R. Ade et al., Astron. Astrophys. 594, A13 (2016).
  * (2) N. Frusciante and L. Perenon, arXiv:1907.03150 [astro-ph.CO]; D. Huterer and D. L Shafer, Rep. Prog. Phys. 81 (2018) 016901; M. J. Mortonson, D. H. Weinberg and M. White, arXiv:1401.0046 [astro-ph.CO]; S. Wang, Y. Wang and M. Li, Phys. Rep. 696 (2017) 1.
  * (3) V. Sahni and A. A. Starobinsky, Int. J. Mod. Phys. D 9, 373 (2000), [astro-ph/9904398]; P. J. E. Peebles and B. Ratra, Rev. Mod. Phys. 75, 559 (2003), [astro-ph/0207347]; S. Weinberg, Rev. Mod. Phys.61, 1 (1989); S. M. Carroll, Living Rev. Relativity 4, 1 (2001), [astro-ph/0004075].
  * (4) V. Faraoni, Cosmology in Scalar-Tensor Gravity, Dordrecht, The Netherlands: Kluwer Academic (2004); C. Armendariz-Picon, T. Damour and V. Mukhanov, Phys. Lett. B 458 (1999) 209; S. Tsujikawa, Class. Quant. Grav. 30 (2013) 214003; J. Gleyzes, D. Langlois, F. Piazza and F. Vernizzi, Phys. Rev. Lett. 114 (2015) 211101; T. P. Sotiriou and V. Faraoni, Rev. Mod. Phys. 82 (2010) 451; A. De Felice and S. Tsujikawa, Living Rev. Rel. 13 (2010) 3.
  * (5) T. Harko, F.S.N. Lobo, S. Nojiri, S.D. Odintsov, Phys. Rev. D 84 (2011) 024020; F. G. Alvarenga, M. J. S. Houndjo, A. V. Monwanou and J. B. C. Orou, J. Mod. Phys. 4 (2013) 130; F.G. Alvarenga, et. al, Phys. Rev. D 87 (2013) 103526; T. Harko and F. S. N. Lobo, Galaxies 2 (2014) 410.
  * (6) Z. Haghani, T. Harko, F.S.N. Lobo, H.R. Sepangi, S. Shahidi, Phys. Rev. D 88 (2013) 044023; S. D. Odintsov and D. Sáez-Gómez, Phys. Lett. B 725 (2013) 437; N. Katırcı and M. Kavuk, Eur. Phys. J. Plus 129 (2014) 163.
  * (7) O. Bertolami, C.G. Boehmer, T. Harko, F.S.N. Lobo, Phys. Rev. D 75 (2007) 104016; T. Harko and F. S. N. Lobo, Extensions of f(R) Gravity Curvature-Matter Couplings and Hybrid Metric-Palatini Theory, Cambridge University Press, UK, 2019.
  * (8) F. W. Hehl, P. von der Heyde, G. D. Kerlick, and J. M. Nester, Rev. Mod. Phys. 48 (1976) 393.
  * (9) E. Cartan, C. R. Acad. Sci. (Paris) 174, 593 (1922); F. W. Hehl, P. von der Heyde, G. D. Kerlick, and J. M. Nester, Rev. Mod. Phys. 48, 393 (1976); F. W. Hehl and B. K. Datta, J. Math. Phys. 12, 1334 (1971); A. Trautman, arXiv:gr-qc/0606062v1.
  * (10) Z. Haghani, N. Khosravi, and S. Shahidi, Classical Quantum Gravity 32 (2015) 215016; J. B. Jimenez and T. S. Koivisto, Classical Quantum Gravity 31 (2014) 135002; Y. Cai, S. Capozziello, M. De Laurentis and E. N. Saridakis, Rept. Prog. Phys. 79 (2016) 106901; Z. Haghani, T. Harko, H. R. Sepangi and S. Shahidi, Phys. Rev. D 88 (2013) 044024; Z. Haghani, T. Harko, H. R. Sepangi and S. Shahidi, JCAP 10 (2012) 061. R. Ferraro, AIP Conf. Proc. 1471 (2012)103.
  * (11) A. Golovnev, V. Mukhanov and V.Vanchurin, JCAP 0806 (2008) 009; A. Golovnev, Classical and Quantum Gravity, 28 (2011) 245018; M. R. Setare, V. Kamali, Phys. Lett. B 08 (2013) 076.
  * (12) S. D. Sadatian, Adv. High Energy Phys. 2013, 802685 (2013); K. Bamba, S. Nojiri and S. D. Odintsov, Phys. Rev. D 77, 123532 (2008); B. Vakili, Phys. Lett. B 738, 488 (2014); J. B. Hartle, S. W. Hawking and T. Hertog, JHEP, 11 (2013) 20; A. Tartaglia and N. Radicella, Phys. Rev. D 76, 083501 (2007); T. Koivisto and D. F Mota, JCAP 08 (2008) 021; R. C. G. Landim, EPJC 76 (2016) 480;
  * (13) M. Born and L. Infeld, Proc. R. Soc. A 144, 425 (1934).
  * (14) Z. Haghani, H. R. Sepangi, S. Shahidi, Phys. Rev. D 83, 064014 (2011); Novello, M., Bergliaffa, S. E. P. and Salim, J. Phys. Rev. D69, 127301 (2004), [arXiv:astro-ph/0312093]; Moniz, P. V. Phys. Rev. D 66, 103501 (2002); GarcaSalcedo, R. and Breton, N. Int. J. Mod. Phys. A 15, 4341 (2000),[arXiv:gr-qc/0004017]; Garca-Salcedo, R. and Breton, N. Class. Quantum Grav. 20, 5425 (2003), [arXiv:hep-th/0212130]; Garca-Salcedo, R. and Breton, N. Class. Quantum Grav. 22, 4783 (2005), [arXiv:grqc/0410142]; Dyadichev, V. V., Galtsov, D. V., Zorin, A. G. and Zotov, M. Yu. Phys. Rev. D 65, 084007 (2002), [arXiv:hep-th/0111099]; Vollick, D. N. Gen. Relativ. Grav. 35, 1511 (2003), [arXiv:hep-th/0102187].
  * (15) F. Bopp, Ann. Phys. (Leipzig) 38, 345 (1940); B. Podolsky, Phys. Rev. 62, 68 (1942).
  * (16) Z. Haghani, T. Harko, H. R. Sepangi and S. Shahidi, Eur. Phys. J. C 77, 137 (2017); S. I. Kruglov, J. Phys. A 43, 245403 (2010); R. B. B. Santos, Modern Physics Letters A 26, 1909 (2011); P. Gaete, International Journal of Modern Physics A 27, 1250061 (2012); A. E. Zayats, Ann. Phys. 342, 11 (2014); V. Perlick, arXiv:1411.0286 (2014); J. Gratus, V. Perlick, and R. W. Tucker, J. Phys. A: Math. Theor. 48, 435401 (2015).
  * (17) A.O. Barvinsky, Mod. Phys. Lett. A 30 (2015), 1540003 arXiv:1408.6112 [hep-th].
  * (18) G. Dvali, G.Gabadadze and M. Porrati, Phys. Lett. B 485, 208 (2000); A. Lue, Phys. Rept. 423, 1 (2006). bibitemnongra L. Modesto and S. Tsujikawa, Phys. Lett. B 727 (2013) 48; L. Modesto and L. Rachwał , Int. J. Mod. Phys. D 26 (2017) 1730020; A. J. Tolley, Lect. Notes Phys. 892 (2015) 203.
  * (19) H. Nersisyan, Y. Akrami, L. Amendola, T. S. Koivisto and J. Rubio, Phys.Rev. D 94 (2016) 043531; X. Zhang, et. al, JCAP 1607 (2016) 003; Y. Dirian et. al, JCAP 1406 (2014) 033; S. Nesserisand S. Tsujikawa, Phys. Rev. D 90 (2014) 024070; M. Maggiore and M. Mancarella, Phys. Rev. D 90 (2014) 023005.
  * (20) N. Arkani-Hamed, S. Dimopoulos, G. Dvali, and G. Gabadadze, hep-th/0209227; G. Dvali, S. Hofmann, and J. Khoury, Phys.Rev. D76 (2007) 084006, hep-th/0703027.
  * (21) M. Jaccard, M. Maggiore, and E. Mitsou, Phys. Rev. D 88, 044033 (2013), arXiv:1305.3034v2 [hep-th].
  * (22) L. Heisenberg, JCAP 05, 015 (2014); L. Heisenberg, arXiv:1801.01523 [gr-qc]; S. Nakamura, R. Kase and S. Tsujikawa, Phys. Rev. D 95, 104001 (2017); R. Kimura, A. Naruko and D. Yoshida, JCAP 1701, 002 (2017); L. Heisenberg, R. Kase and S. Tsujikawa, Phys. Lett. B 760, 617 (2016);J. B. Jimenez and L. Heisenberg, Phys. Lett. B 757, 405 (2016); M. Hull, K. Koyama and G.Tasinato, Phys. Rev. D 9, 064012 3 (2016); N. Khosravi, Phys. Rev. D 89, 124027 (2014); J. B. Jimenez, arXiv:1606.04361 [gr-qc]; Z. Haghani, T. Harko and S. Shahidi, Phys. Dark Univ. 21 (2018) 27, arXiv:1707.00939v2 [gr-qc].
  * (23) E. Allys, P. Peter and Y. Rodriguez, Phys. Rev. D 94, 084041 (2016), arXiv:1609.05870 [hep-th]; C. Armendariz-Picon, JCAP 0407, 007 (2004); R. C. G. Landim, Eur. Phys. J. C 76, 480 (2016).
  * (24) M. Nedaie and Z. Haghani, Int. J. Phys. Research 19, 463 (2019).
  * (25) M. Maggiore, “ A Modern Introduction to Quantum Field Theory”, Oxford University Press, (2005).
  * (26) Virgo, Fermi-GBM, INTEGRAL, LIGO Scientific Collaboration, B. P. Abbott et. al., Astrophys. J. 848 (2017) L13; P. Creminelli and F. Vernizzi, Phys. Rev. Lett. 119 (2017) 251302; J. Sakstein and B. Jain, Phys. Rev. Lett. 119 (2017) 251303; J. María Ezquiaga and M. Zumalacárregui, Phys. Rev. Lett. 119 (2017) 251304.
  * (27) A. Pais and G.E. Uhlenbeck, Phys. Rev. 79, 145 (1950).
  * (28) T. J. Chen, M. Fasiello, E. A. Lim and A. J. Tolley, JCAP 1302 (2013) 042.
  * (29) C. Deffayet, A. E. Gumrukcuoglu, S. Mukohyama and Y. Wang, JHEP 04 (2014) 082.