1 Introduction

Homogeneity and isotropy, together with matter being treated as a specific gas, are the basic ingredients of the Friedmann–Lemaître–Robertson–Walker (FLRW) model. Since the pioneering work of E. Hubble, who suggested a correlation between the observed redshifts and distances of 24 galaxies, this is a very successfull model, being currently considered a very good measure to describe our cosmos.

In the last decades a great effort was dedicated to understanding the local inhomogeneities that occur in the Universe. Such efforts suggest that these inhomogeneities may be related to the expansion of the Universe. Inhomogeneities, as an alternative to dark energy, were first discussed in Ref. [1] as a means of explaning the observational results of the expansion of the Universe [2,3,4,5,6,7,8] without the need for postulating dark energy. Inhomogeneities inside astronomical objects can define different instability ranges, an effect that can describe distinct features of their evolution and structures formation.

The current standard model of cosmology, the \(\Lambda \)-Cold Dark Matter (\(\Lambda \hbox {CDM}\)) model, is a spatially homogeneous solutions of the FLRW Einstein’s field equations with only six free parameters that successfully accounts for most cosmological data, especially the characteristics of the Cosmic Microwave Background (CMB) and the structure formation on large scales considered through the theory of cosmological perturbations in spatially homogeneous and isotropic background. However, in the last fifteen years “standard” inhomogeneous cosmological models that are generalizations of FLRW cosmologies have been the subject of growing interest in astrophysics community in order to investigate cosmological phenomena. Some authors have demonstrated that spatially inhomogeneous models with spherical symmetry and dust source can be fitted to supernovae Ia (SNIa) data, as well as the position of the first peak of the CMB. These models show that the apparent accelerated expansion of the universe may not be a consequence of the repulsive gravity due to dark energy, but rather the result of inhomogeneities in the distribution of matter. In this context one needs to mention that an important cosmological model to describe the inhomogeneous universe is the Lemaître–Tolman–Bondi (LTB) spacetime [11,12,13,14,15,16,17,18,19], which is a spatially inhomogeneous description of a spherically symmetric distribution of dust matter in the Universe.

In the last few years several studies analyzing the scenario of electrodynamics embedded in the isotropic and homogeneous FLRW gravitational background models [20] have been produced. Astrophysical effects were explored together with an anisotropic expansion of the Universe. The interesting result of polarized electromagnetic radiation occurs when it travels through local anisotropic regions. Concerning the inhomogeneities, in Ref. [21] the authors investigated that the inhomogeneity with electromagnetic field caused a new scale factor. The propagation of photons was also affected, which is important phenomenum since most information obtained about the Universe is by means of photons. Reference [22] studied the effects of Palatini f(R) gravity together with the so-called tilted observer on the dynamics of LTB spacetime embedded in an electromagnetic field. The Proca theory in the FLRW cosmology has also been studied in Refs. [23,24,25,26].

In this work we studied the effect of a high-derivative electromagnetic field embedded in the inhomogeneous LTB geometry. We analyzed the Proca and high-derivatives Podolsky electromagnetic models. The electrodynamic contribution was inserted separately through the energy-momentum tensor (EMT) of the respective models. For the Proca model, besides the matter density there also are the electromagnetic contributions obtained from the Lagrangian of the Proca model in curved space-time. An analysis of the Proca model in curved space-time was carried out for a particular case of interest by Bekenstein [9]. Here, we have analyzed the electrodynamics effects in LTB cosmological model and calculate the scale factor in LTB universe. We also computed the luminosity distance in the presence of electromagnetic field.

Concerning the Podolsky electrodynamics effects in LTB cosmological model, the study of Podolsky model in curved space-time was made in Ref. [10] together with the analysis of the Bopp–Podolsky black holes. Therefore, we have obtained the line element and the equations that define the LTB model with Podolsky contributions.

This paper is structured as follows. In Sect. 2 we review the main aspects of the LTB metric, where the Einstein tensor and the matter contribution for the EMT are obtained. In Sect. 3 we analyze the Proca model of electrodynamics in curved space-time, where the EMT and the Maxwell–Proca equations in curved space-time are derived. In Sect. 4 we solve the Einstein equations including the Proca contribution. We obtain the scale factor for this model and examine the inhomogeneities and luminosity distance. In Sect. 5 we analyze the Podolsky electrodynamics contributions for LTB model. Section 6 discusses the results and presents some final considerations.

2 The LTB cosmological model: a brief review

The LTB model depicts a self-gravitating spherically symmetric distribution of spatially inhomogeneous nondissipative dust cloud where the EMT can be written as \(T_{\mu \nu } = \epsilon (\tau ,\rho )u_\mu u_\nu \) and \(u_\mu = u_\mu (\tau ,\rho )\) is the dust particle’s four-velocity vector. The proper time is represented by \(\tau \) and the several shells are labeled by \(\rho \), which helps in the formation of the dust cloud. As mentioned above, the LTB metric is the most popular inhomogeneous cosmological model, since it is appropriate for both small and large scale inhomogeneities, and also is the simplest inhomogeneous solution of the Einstein equations.

In this section we will analyze the LTB model through the point of view of the field equations of general relativity, as well as the matter contribution concerning the EMT. In LTB cosmology, the line element is written as

$$\begin{aligned} dS^2=dt^2-e^Adr^2-R^2\Big (d\theta ^2+sin^2{\theta }d\phi ^2\Big ), \end{aligned}$$
(1)

where A and R, unlike FLRW model, depend also on the r coordinate, namely, \(A{\equiv }A(r,t)\) and \(R{\equiv }R(r,t)\). Hence, the metric element is given by

$$\begin{aligned} g_{\mu \nu }=\Big (1,-e^A,-R^2,-R^2sin^2\theta \Big ), \end{aligned}$$
(2)

where

$$\begin{aligned} \sqrt{-g}=e^{A/2}R^2sin\theta . \end{aligned}$$
(3)

The spatially homogeneous FLRW metric is a special case of the metric in Eq. (1), where

$$\begin{aligned} e^{A(r,t)} \,=\,\frac{a(t)}{\sqrt{1\,-\,kr^2}} \end{aligned}$$
(4)

and \(R(r,t) = a(t) r\), where a(t) is the scale factor.

The Einstein equations, where \(c=1\), are given by

$$\begin{aligned} G_{\mu \nu }=8{\pi }GT_{\mu \nu }, \end{aligned}$$
(5)

where \(G_{\mu \nu }\) is the Einstein tensor

$$\begin{aligned} G_{\mu \nu }=R_{\mu \nu }-\frac{1}{2}g_{\mu \nu }R, \end{aligned}$$
(6)

G is the gravitational constant and \(T_{\mu \nu }\) is the EMT. So, we can make use of the line element in Eq. (1) to obtain the non-zero components of the Einstein tensor,

$$\begin{aligned} G^0_{\;0}&=\Big (-2\frac{R''}{R}{+}\frac{R'A'}{R}-\frac{R'^2}{R}\Big )e^{-A}{+} \Big (\frac{{\dot{R}}{\dot{A}}}{R}{+}\frac{{\dot{R}}^2}{R^2}+\frac{1}{R^2} \Big ),\nonumber \\ G^1_{\;1}&=\frac{{\dot{R}}^2}{R^2}+2\frac{{\ddot{R}}}{R}-e^{-A}\frac{R'^2}{R^2}+\frac{1}{R^2},\nonumber \\ G^2_{\;2}&=G^3_{\;3}=\Big (-\frac{R''}{R}+\frac{A'R'}{2R}\Big )e^{-A}+ \frac{{\ddot{R}}}{R}+\frac{{\ddot{A}}}{2}\nonumber \\&\quad +\frac{{\dot{A}}^2}{4}- \frac{{\dot{R}}{\dot{A}}R}{2},\nonumber \\ G^1_{\;0}&=e^{-A}\Big (2\frac{{\dot{R}}'}{R}-\frac{{\dot{A}}R'}{R}\Big ). \end{aligned}$$
(7)

In LTB cosmology, the EMT is the zero pressure diagonal perfect fluid, \((\rho ,0,0,0)\), where \(\rho \) is the mass density. So, the matter contribution for the EMT \(T_{00}^M\) is

$$\begin{aligned} T_{00}^M=\rho . \end{aligned}$$
(8)

The electromagnetic contribution for the LTB cosmology will be obtained by adding on the right side of the Einstein equations. So, we have for the EMT, for both the Proca model and for the Podolsky model, that

$$\begin{aligned} T_{\mu \nu }=T_{\mu \nu }^M+T_{\mu \nu }^{EM}, \end{aligned}$$
(9)

where \(T_{\mu \nu }^{EM}\) is the electromagnetic contribution for the EMT. Therefore, determining \(T_{\mu \nu }^{EM}\), then we have the Einstein equation.

3 The curved spacetime Proca model

The Proca model of electrodynamics is defined by the following Lagrangian

$$\begin{aligned} {\mathcal {L}}_{Proca}=-\frac{1}{4}F_{\mu \nu }F^{\mu \nu }+\frac{1}{2} m^2A_{\mu }A^{\mu }, \end{aligned}$$
(10)

where m is the mass of electromagnetic field, \(A_\mu \) is the four-potential and \(F_{\mu \nu }\) is the electromagnetic tensor that, in curved space-time, is defined by \(F_{\mu \nu }=\nabla _{\mu } A_\nu -\nabla _{\mu }A_\nu \), where \(\nabla \) is the covariant derivative given by

$$\begin{aligned} \nabla _{\gamma }F_{\alpha \beta }=\partial _{\gamma }F_{\alpha \beta }- \Gamma ^{\rho }_{\gamma \alpha }F_{\rho \beta }-\Gamma ^{\rho }_{\gamma \beta } F_{\alpha \rho }. \end{aligned}$$
(11)

Moreover, the electromagnetic field and its dual can be obtained from the following expressions,

$$\begin{aligned} F_{\mu \nu }&=u_{\mu }E_\nu -u_{\nu }E_\mu +\epsilon _{\mu \nu \alpha \beta } B^\alpha {u}^\beta \nonumber \\ ^*F^{\mu \nu }&=\frac{1}{2}\epsilon ^{\mu \nu \alpha \beta }F_{\alpha \beta } \end{aligned}$$
(12)

and the Maxwell–Proca equations in curved spacetime are

$$\begin{aligned} \nabla _{\nu }F^{\mu \nu }+m^2A^{\mu }&=4\,{\pi }J^\mu \nonumber \\ {\nabla _{\nu }}^*F^{\mu \nu }&=0, \end{aligned}$$
(13)

where \(J^\mu =(\rho ,\vec {J})\) is the electromagnetic four-current density and \(^*F^{\mu \nu }\) is the dual electromagnetic tensor.

From Eq. (10), the action for the Proca model is

$$\begin{aligned} S=\frac{1}{16\pi }{\int }d^4x\sqrt{-g}\Big (-R-4{\mathcal {L}}\Big ), \end{aligned}$$
(14)

namely,

$$\begin{aligned} S=\frac{1}{16\pi }{\int }d^4x\sqrt{-g}\Big (-R+F_{\mu \nu }F^{\mu \nu }-2m^2A_{\mu }A^{\mu }\Big ). \end{aligned}$$
(15)

and, from the EMT definition

$$\begin{aligned} T_{\mu \nu }=\frac{2}{\sqrt{-g}}\frac{\delta {S}}{\delta {g^{\mu \nu }}}, \end{aligned}$$
(16)

the Proca EMT in curved spacetime is

$$\begin{aligned} T_{\mu \nu }= & {} \frac{1}{4\pi }\Big [F_{\mu }^{\rho }F_{\nu \rho }+m^2A_{\mu } A_{\nu }+g_{\mu \nu }\Big (\frac{1}{4}F_{\alpha \beta }F^{\alpha \beta }\nonumber \\&- \frac{1}{2}m^2A_{\alpha }A^{\alpha }+J^\alpha {\tilde{A}}_\alpha \Big )\Big ], \end{aligned}$$
(17)

In this work, we neglect the coupling term \(J^\alpha {\tilde{A}}_\alpha \) because most of the matter is electrically neutral. Thus, any fluctuations that may appear can be ignored. The trace of the energy-momentum tensor is:

$$\begin{aligned} T=-\frac{5}{8\pi }m^2A_{\mu }A^{\mu }, \end{aligned}$$
(18)

i.e., unlike Maxwell’s case, the tensor trace is not zero.

Using relations in Eq. (17), we find that the non-zero components of the EMT are

$$\begin{aligned} T_{00}&=e^A(E^{(1)})^2+R^2(E^{(2)})^2+R^2sin^2(E^{(3)})^2\nonumber \\&\quad +\frac{1}{4}F^{\alpha \beta }F_{\alpha \beta }+\frac{m^2}{4\pi }A_0^2-\frac{1}{8\pi }m^2A_{\alpha }A^\alpha ,\nonumber \\ T_{11}&=-e^{2A}(E^{(1)})^2-e^AR^2sin^2(B^{(3)})^2+e^AR^2(B^{(2)})^2\nonumber \\&\quad +\frac{e^A}{4}F^{\alpha \beta }F_{\alpha \beta }+\frac{m^2}{4\pi }A_1^2-\frac{1}{8\pi }e^Am^2A_{\alpha }A^\alpha ,\nonumber \\ T_{22}&=-R^4E^{(2)}+R^4sin^2\theta (B^{(3)})^2+e^AR^2(B^{(1)})^2\nonumber \\&\quad -\frac{R^2}{4}F^{\alpha \beta }F_{\alpha \beta }+\frac{m^2}{4\pi }A_2^2 -\frac{1}{8\pi }g_{22}m^2A_{\alpha }A^\alpha \nonumber \\ T_{33}&=-R^4E^{(2)}+R^4sin^2\theta (B^{(3)})^2+e^AR^2(B^{(1)})^2 \nonumber \\&\quad -\frac{R^2}{4}F^{\alpha \beta }F_{\alpha \beta }+\frac{m^2}{4\pi }A_3^2-\frac{1}{8\pi }g_{33}m^2A_{\alpha } A^\alpha \nonumber \\ T_{10}&=e^{A/2}R^2sin\theta (E^{(3)}B^{2)}-E^{(2)}B^{(3)}) +\frac{m^2}{4}A_1A_0,\nonumber \\ T_{20}&=e^{A/2}R^2sin\theta (E^{(1)}B^{(3)}-E^{(3)}B^{(1)})+\frac{m^2}{4}A_2A_0,\nonumber \\ T_{30}&=e^{A/2}R^2sin\theta (E^{(2)}B^{(1)}-E^{(1)}B^{(2)})+\frac{m^2}{4}A_3A_0,\nonumber \\ T_{12}&=-e^{A}R^2(E^{(1)}E^{(2)}+B^{(1)}B^{(2)})+\frac{m^2}{4}A_1A_2,\nonumber \\ T_{13}&=-e^{A}R^2sin^2\theta (E^{(1)}E^{(3)}+B^{(1)}B^{(3)})+\frac{m^2}{4}A_1A_3,\nonumber \\ T_{23}&=-R^4sin^2\theta (E^{(2)}E^{(3)}-B^{(2)}B^{(3)})+\frac{m^2}{4}A_2A_3. \end{aligned}$$
(19)

The next step is introduce Eq. (19) into the Einstein equations in Eq. (5). We will do this in the next section where we will analyze the effects of the Proca model on the LTB model.

4 LTB scenario Proca electrodynamics

Now that we have the Proca’s contribution in hand, we must analyze the terms of the Einstein tensor to obtain Einstein’s equations. As we can see from Eq. (7), the only non-zero off-diagonal term in \(G_{\mu \nu }\) is \(G_{10}\). The only way to satisfy this condition is to impose that

$$\begin{aligned} T_{02}=T_{03}=T_{12}=T_{13}=T_{23}, \end{aligned}$$
(20)

and, consequently

$$\begin{aligned} {\left\{ \begin{array}{ll}A_2=A_3=0\\ F_{02}=F_{03}=F_{12}=F_{13}=F_{23}=0.\end{array}\right. } \end{aligned}$$
(21)

Taking into account the last conditions, the electromagnetic fields are

$$\begin{aligned} E^\mu&=(0,E(t,r),0,0),\nonumber \\ B^\mu&=(0,0,0,0), \end{aligned}$$
(22)

as we can see, as consequence of conditions in Eq. (21), we have that the \(B_\mu \) field is zero. The only contribution of electromagnetic field is by \(E_\mu \) field. Now, imposing these same conditions in Eq. (19), the energy-momentum components are

$$\begin{aligned} T_{00}&=\frac{1}{2}e^{A}E^2+\frac{m^2}{8\pi }A_0^2-\frac{m^2}{8\pi }g^{11}A_1^2,\nonumber \\ T_{11}&=\frac{1}{2}e^{2A}E^2-\frac{m^2e^A}{8\pi }\Big (A_0^2+A_1^2\Big ),\nonumber \\ T_{22}&=T_{33}=-\frac{R^2e^A}{2}E^2-g_{33} \frac{1}{8\pi }m^2A_{0}^2,\nonumber \\ T_{01}&=\frac{1}{8\pi }\Big (m^2A_0A_1\Big ). \end{aligned}$$
(23)

In order to write A in function of R, we impose that \(A_1=0\), that is \(A_\mu =(A_0,0,0,0)\). So, the Eq. (19) turns

$$\begin{aligned} T_{00}&=\frac{1}{2}e^{A}E^2+\frac{m^2}{8\pi }A_0^2,\nonumber \\ T_{11}&=\frac{1}{2}e^{2A}E^2-\frac{m^2e^A}{8\pi }\Big (A_0^2\Big ),\nonumber \\ T_{22}&=T_{33}=-\frac{R^2e^AE^2}{2}-g_{33}\frac{1}{8\pi }m^2A_{0}^2,\nonumber \\ T_{01}&=0. \end{aligned}$$
(24)

and the Einstein equations with matter and electromagnetic contribution in EMT are

$$\begin{aligned} G^0_0&=\Big (-2\frac{R''}{R}+\frac{R'A'}{R}-\frac{R'^2}{R}\Big )e^{-A}\nonumber \\&\quad +\Big (\frac{{\dot{R}}{\dot{A}}}{R}+\frac{{\dot{R}}^2}{R^2}+\frac{1}{R^2}\Big )\nonumber \\&=8{\pi }G\Big [\rho _m+\frac{1}{2}e^{A}E^2+\frac{m^2}{8\pi }A_0^2\Big ],\nonumber \\ G^1_1&=-\frac{{\dot{R}}^2}{R^2}+2\frac{{\ddot{R}}}{R}-e^{-A}\frac{R'^2}{R^2}\nonumber \\&\quad +\frac{1}{R^2}=-8{\pi }G\Big [\frac{1}{2}e^{A}E^2-\frac{m^2}{8\pi }\Big (A_0^2\Big )\Big ],\nonumber \\ G^2_2&=G^3_3=\Big (-\frac{R''}{R}+\frac{A'R'}{2R}\Big )e^{-A}\nonumber \\&\quad +\frac{{\ddot{R}}}{R}+\frac{{\ddot{A}}}{2}+\frac{{\dot{A}}^2}{4}-\frac{{\dot{R}}{\dot{A}}R}{2}\nonumber \\&=-g^{33}\frac{R^2e^AE^2}{2}-\frac{1}{8\pi }m^2A_{0}^2,\nonumber \\ G^1_0&=e^{-A}\Big (2\frac{{\dot{R}}'}{R}-\frac{{\dot{A}}R'}{R}\Big )=0. \end{aligned}$$
(25)

The last expression of Eqs. (25) allows us as to write A in terms of \(R'\) such as

$$\begin{aligned} e^A=\frac{R'^2}{1-k(r)} \end{aligned}$$
(26)

where k(r) is an arbitrary function. So, the line element is

$$\begin{aligned} dS^2=dt^2-\frac{R'^2}{1-k(r)}dr^2-R^2\Big (d\theta ^2+sin^2{\theta }d\phi ^2\Big )\,\,. \end{aligned}$$
(27)

Moreover, we can rewrite the Einstein equations for the (00) and (11) components like

$$\begin{aligned}&\frac{{\dot{R}}^2+k}{R^2} +\frac{2{\dot{R}}{\dot{R}}'+k'}{RR'}\nonumber \\&\quad = 8{\pi }G\, \bigg [\rho _m+\frac{R'^2E^2}{(1-k)}+m^2A_0^2-J^{\alpha }A_\alpha \bigg ], \end{aligned}$$
(28)
$$\begin{aligned}&\frac{{\dot{R}}^2+2R{\ddot{R}}+k}{R^2}\nonumber \\&\quad = 8{\pi }G\bigg [\frac{R'^2E^2}{(1-k)}+m^2A_0^2-J^{\alpha }A_\alpha \bigg ]. \end{aligned}$$
(29)

On the other hand, the non-zero components of the Maxwell–Proca equations in Eq. (12) are

$$\begin{aligned}&\partial _1\Big (e^{A/2}R^2F^{01}+m^2e^{A/2}R^2A^0\Big )=4{\pi }e^{A/2}R^2J^0, \end{aligned}$$
(30)
$$\begin{aligned}&\partial _0\Big (e^{A/2}R^2F^{10}\Big )=4{\pi }e^{A/2}R^2J^1, \end{aligned}$$
(31)
$$\begin{aligned}&\partial _2E^{(1)}=0, \end{aligned}$$
(32)
$$\begin{aligned}&\partial _3E^{(1)}=0. \end{aligned}$$
(33)

So, we can write the electric field E(t, r) using the integration of Eq. (30), as

$$\begin{aligned} E(t,r)=\frac{\epsilon (t)+\xi (t,r)}{e^{A/2}R^2}, \end{aligned}$$
(34)

where

$$\begin{aligned} \xi (t,r)\equiv 4\pi \int _0^re^{A/2}R^2J^0d{\bar{r}}-m^2e^{A/2}R^2A^0, \end{aligned}$$
(35)

and \(\epsilon (t)\) is the constant of integration. Now, from Eq.  (31), using the fact that the charged matter belongs to the comoving matter, namely, the four-current appears as \(J^\alpha =(J^0,\vec {0})\)

$$\begin{aligned} {\dot{\epsilon }}+{\dot{\xi }}=0, \end{aligned}$$
(36)

which leads us to

$$\begin{aligned} E(t,r)=\frac{\epsilon _0+\xi _0(r)}{e^{A}R^4}. \end{aligned}$$
(37)

Now, defining

$$\begin{aligned} \sigma (r)\,{\equiv }\,4{\pi }Ge^AR^4(E^2+2e^{-A}m^2A^2_0), \end{aligned}$$
(38)

multiplying the Eq. (29) by \({\dot{R}}\), we obtain that

$$\begin{aligned} {\dot{R}}^3+2A{\dot{R}}{\ddot{R}}=-k{\dot{R}}+\sigma \frac{{\dot{R}}}{R^2}, \end{aligned}$$
(39)

or yet, using the fact that \(\partial _0(R{\dot{R}})=2A{\dot{R}}{\ddot{R}}\), we have

$$\begin{aligned} R{\dot{R}}&=-k(r)R+\alpha (r)+\int \frac{\sigma (r)}{R^2}{\dot{R}}dt,\nonumber \\ {\dot{R}}^2&={-k(r)}+\frac{\alpha (r)}{R}-\frac{\sigma (r)}{R^2}, \end{aligned}$$
(40)

and using Eq. (28) we can write

$$\begin{aligned}&8\pi {G}\Big (\rho +\frac{R'E^2}{(1-k)}+m^2A_0^2\Big )\nonumber \\&\quad =\frac{k'}{RR'}+\frac{\alpha (r)}{R^3} -\frac{\sigma (r)}{R^4}+\frac{2}{RR'}\Big [\Big (-k(r)\nonumber \\&\qquad -\frac{\alpha (r)}{R}-\frac{\sigma (r)}{R^2}\Big )\Big (\frac{\alpha {R'}}{R^2}+2\frac{\sigma (r)R'}{R^3}\Big )\Big ] \end{aligned}$$
(41)

Equations (27), (40) and (41) define the LTB model with the Proca electrodynamics contributions. The singularities arising from \(R=0\), \(R'=0\) and \(k=1\). The \(R=0\) singularity we interpreted as the Big Bang singularity, and the \(R'=0\) singularity as the shell cross singularity. The last singularity, \(k=1\), that occurs if \(R'\ne {0}\), come from the Proca contribution.

Now, defining

$$\begin{aligned} -k(r)\,{\equiv }\,H_0^2(r)\Omega _k(r)R_0^2(r),\nonumber \\ \alpha (r)\,{\equiv }\,H_0^2(r)\Omega _m(r)R_0^3(r),\nonumber \\ -\sigma (r)\,{\equiv }\,H_0^2(r)\Omega _\sigma (r)R_0^4(r), \end{aligned}$$
(42)

where the Hubble constant is defined as

$$\begin{aligned} H(t,r)=\frac{{\dot{R}}}{R}. \end{aligned}$$
(43)

We have imposed the boundary values at \(t_0\) through \(A_0(r){\equiv }A(t_0,r)\), \(H_0(r){\equiv }H(t_0,r)\). Moreover, \(\Omega _k\), \(\Omega _m\) and \(\Omega _\sigma \) are subjected to the constraint

$$\begin{aligned} \Omega _k(r)+\Omega _\sigma (r)+\Omega _m(r)=1. \end{aligned}$$
(44)

Therefore, we can write Eq. (40) as

$$\begin{aligned} \frac{{\dot{R}}}{R}=H_0(r)\Big [\Omega _k(r)\Big (\frac{R_0}{R}\Big )^2+\Omega _m(r)\Big (\frac{R_0}{R}\Big )^3+\Omega _\sigma (r)\Big (\frac{R_0}{R}\Big )^4\Big ]^{1/2} \end{aligned}$$
(45)

and the Hubble constant is computed by

$$\begin{aligned} -H_0(r)t=-\int _{\frac{R}{R_0}}^1\frac{dx}{\sqrt{\Omega _k(r)+\Omega _\sigma (r)x^{-1}+\Omega _m(r)x^{-2}}} \end{aligned}$$
(46)

and, considering the solution using that \({\dot{R}}>0\), we have

$$\begin{aligned} R(t,r)=R_0(r)\bigg [\frac{\Omega _\sigma (r)}{\Omega _m(r)}+\frac{\Omega _m(r)\Omega _\sigma ^2(r)}{M(t,r)}+\frac{M(t,r)}{\Omega _m^3(r)}\bigg ]\,\,, \end{aligned}$$
(47)

where

$$\begin{aligned}&M(t,r)=\Omega _m^2(r)\bigg (\frac{N+2\Omega _\sigma ^3(r)+\sqrt{N^2+4N\Omega _\sigma ^6(r)}}{2}\bigg )^{1/3}\,\,, \end{aligned}$$
(48)
$$\begin{aligned}&N(t,r)=\bigg [\frac{3\Omega _m(r)H_0(r)t}{2}+\Omega _m(r)-2\Omega _\sigma (r)\bigg ]^2{-}4\Omega _\sigma ^3(r), \end{aligned}$$
(49)

where \(R_0(r)\) corresponds to the current shape of the scale factor, \(H_0(r)\) is the current value of the Hubble constant in each point and \(\Omega _\sigma (r)\) is the density of the electromagnetic field, that, because of Proca’s contributions, is purely electric.

Remember that the density of the electromagnetic field is given by

$$\begin{aligned} \Omega _\sigma =-\frac{\sigma (r)}{H_0^2R_0^4} \end{aligned}$$
(50)

where

$$\begin{aligned} \sigma (r)\,{\equiv }\,4{\pi }Ge^AR^4\Big (E^2+2e^{-A}m^2A^2_0\Big )\,\,, \end{aligned}$$
(51)

and the Proca contribution for the density of the electromagnetic field is the mass term \(-2e^{-A}m^2A^2_0\).

Now, let us examine the inhomogeneities and luminosity distance for the Proca electrodynamics perspective. Therefore, the radial geodesic requires that \(d\theta =d\phi =0\). Moreover, since light always travels along null geodesics, we have \(dS^2=0\). So, using Eq.  (27), we have

$$\begin{aligned} dt^2=\frac{R'^2}{1+2k(r)}dr^2 \end{aligned}$$
(52)

so that

$$\begin{aligned} \frac{dt}{du}=\frac{R'}{\sqrt{1+2k(r)}}\frac{dr}{du}\,\,. \end{aligned}$$
(53)

Consider two light rays with solutions of Eq. (53) given by \(t_1 = t(u)\) and \(t_2 = t(u) + \lambda (u)\). Substituting these two into Eq. (53) we obtain

$$\begin{aligned} \frac{dt_1}{du}&=-\frac{dr}{du}\frac{R'}{\sqrt{1+2k(r)}},\nonumber \\ \frac{dt_2}{du}&=-\frac{dr}{du}\frac{R'+{\dot{R}}\lambda }{\sqrt{1+2k(r)}},\nonumber \\ \frac{d\lambda }{du}&=-\frac{dr}{du}\frac{{\dot{R}}'\lambda }{\sqrt{1+2k(r)}}. \end{aligned}$$
(54)

Differentiating the definition of the redshift, \(z\,{\equiv }\,[\lambda (0)-\lambda (u)]/\lambda (u)\), we have,

$$\begin{aligned} \frac{dz}{du}=\frac{dr}{du}\frac{(1+z){\dot{R}}'}{\sqrt{1+2k(r)}} \end{aligned}$$
(55)

and, consequently,

$$\begin{aligned} \frac{dt}{dz}&=-\frac{R'}{(1+z){\dot{R}}'}\,\,,\nonumber \\ \frac{dr}{dz}&=\sqrt{\frac{1-k(r)}{(1+z){\dot{R}}'}}\,\,, \end{aligned}$$
(56)

which determine the relation between the coordinates and the observable redshift.

For the last equation, using the expression for k(r)

$$\begin{aligned} k(r)\,{\equiv }\,H_0^2\Omega _kA_0^2 \end{aligned}$$
(57)

we find

$$\begin{aligned} \frac{dr}{dz}=\sqrt{\frac{1+H_0^2(r)(1-\Omega _m(r)+\Omega \sigma (r))R_0^2(r)}{(1+z){\dot{R}}'(t,r)}}\,\,. \end{aligned}$$
(58)

Moreover, the relation between the redshift and the energy flux F, defined as \(d_L{\equiv }\sqrt{L/(a\pi {F})}\), where L is the total power distance radiated by the source, is

$$\begin{aligned} d_L(z)=(1-z)^2R(r(z),t(z)) \end{aligned}$$
(59)

and the angular distance diameter is given by

$$\begin{aligned} d_A(z)=R(r(z),t(z)) \end{aligned}$$
(60)

which is a direct relation to the scale factor as a function of the redshift.

5 LTB geometry Podolsky electrodynamics

The Podolsky electrodynamics in curved space-time is given by

$$\begin{aligned} {\mathcal {L}}_{Pod}=-\frac{1}{4}F^{\alpha \beta }F_{\alpha \beta }+\frac{a^2}{2}\nabla _{\beta }F^{\alpha \beta }\nabla _{\gamma }F_\alpha ^{\;\gamma }+\frac{b^2}{2}\nabla _{\beta }F^{\alpha \gamma }\nabla ^{\beta }F_{\alpha \gamma } \end{aligned}$$
(61)

or,

$$\begin{aligned} {\mathcal {L}}_{Pod}&=-\frac{1}{4}F^{\alpha \beta }F_{\alpha \beta }+\frac{a^2+2b^2}{2}\nabla _{\beta }F^{\alpha \beta }\nabla _{\gamma }F_\alpha ^{\;\gamma }\nonumber \\&\quad +\frac{b^2}{2}\Big (R_{\sigma \beta }F^{\sigma \alpha }F_\alpha ^{\;\beta }+R_{\alpha \sigma \beta \gamma }F^{\sigma \gamma }F^{\alpha \beta }\Big ) \end{aligned}$$
(62)

where \(R_{\sigma \beta }\) is the Ricci tensor.

The Einstein–Podolsky action is

$$\begin{aligned} S=\frac{1}{16\pi }{\int }d^4x\sqrt{-g}\Big [-R+4{\mathcal {L}}_{Pod}\Big ] \end{aligned}$$
(63)

which leads us, from the variation with respect to \(A_\mu \), to the Einstein–Podolsky equation

$$\begin{aligned} \nabla _\nu \Big [F^{\mu \nu }-(a^2+2b^2)H^{\mu \nu }+2b^2S^{\mu \nu }\Big ]=0 \end{aligned}$$
(64)

where

$$\begin{aligned} H^{\mu \nu }&\,{\equiv }\,\nabla ^{\mu }K^\nu -\nabla ^{\nu }K^\mu \nonumber \\ S^{\mu \nu }&\,{\equiv }\,F^{\mu \sigma }R_\sigma ^\nu -F^{\nu \sigma }R_\sigma ^\mu +2R^{\mu \;\nu }_{\;\sigma \;\beta }{F^{\beta \sigma }} \end{aligned}$$
(65)

and

$$\begin{aligned} K^\mu \,{\equiv }\,\nabla _{\nu }F^{\mu \nu }. \end{aligned}$$
(66)

The EMT for the Podolsky electrodynamics in curved space-time is given by

$$\begin{aligned} T_{\mu \nu }&=\frac{1}{4\pi }\Big [F_{\mu \sigma }F^\sigma _{\;\nu }+\frac{1}{4}g_{\mu \nu }F^{\alpha \beta }F_{\alpha \beta }\Big ]\nonumber \\&\quad +\frac{a^2}{4\pi }\Big [g_{\mu \nu }F_\beta ^{\;\gamma }\nabla _{\gamma }K^\beta +\frac{1}{2}g_{\mu \nu }K^{\beta }K_{\beta }\nonumber \\&\quad +2F_{(\mu }^\alpha \nabla _{\nu )}K_\alpha -2F_{(\mu }^\alpha \nabla _{\alpha }K_{\nu )}-K_\mu {K_\nu }\Big ]\nonumber \\&\quad +\frac{b^2}{2\pi }\Big [\frac{1}{4}g_{\mu \nu }\nabla ^{\beta }F^{\alpha \gamma } \nabla _{\beta }F_{\alpha \gamma }+F^\gamma _{\;(\mu }\nabla ^\beta \nabla _{\beta }F_{\nu )\gamma }\nonumber \\&\quad +F^{\beta }_{(\mu }\nabla _\beta \nabla _{\nu )}-\nabla _{\beta }\Big (F_\gamma ^{\;\beta }\nabla _{(\mu }F_{\nu )}^{\;\gamma }\Big )\Big ] \end{aligned}$$
(67)

and the trace of the EMT is

$$\begin{aligned} T=\frac{a^2+2b^2}{4\pi }K^{\mu }K_\mu +\frac{b^2}{4\pi }\Big [2\nabla _\beta \Big (F^{\gamma \alpha }\nabla ^{\beta }F_{\alpha \gamma }\Big )+F^{\nu \mu }S_{\mu \nu }\Big ] \end{aligned}$$
(68)

Now, using the metric element Eq. (1), the non-null elements for the energy-momentum tensor are

$$\begin{aligned} T^0_{\;0}&=e^A(E^{(1)})^2+R^2(E^{(2)})^2+R^2sin^2(E^{(3)})^2\nonumber \\&\quad +\frac{1}{4}F^{\alpha \beta }F_{\alpha \beta } +\frac{a^2}{4\pi }\Big [F_\beta ^{\;\gamma }\nabla _{\gamma } K^\beta \nonumber \\&\quad -\frac{K_0K^0}{2}+\frac{1}{2}K^iK_i+2F^{0i}\nabla _0K_i-2F^{0i} \nabla _iK_0\Big ]\nonumber \\&\quad +\frac{b^2}{2\pi }\Big [-\frac{1}{4}\nabla ^{\beta }F^{\alpha \gamma } \nabla _{\beta }F_{\alpha \gamma }\nonumber \\&\quad +F^{i0}\nabla ^\beta \nabla _{\beta }F_{0i}+F_{0i} \nabla _\beta \nabla ^0{F}^{\beta {i}}-\nabla _\beta \Big (F_i^\beta \nabla ^0{F}_0^i\Big )\Big ] \end{aligned}$$
(69)
$$\begin{aligned} T^1_{\;1}&=e^A(E^{(1)})^2-R^2sin^2(B^{(3)})^2-R^2(B^{(2)})^2 \nonumber \\&\quad +\frac{1}{4}F^{\alpha \beta }F_{\alpha \beta } +\frac{a^2}{4\pi }\Big [F_\beta ^{\;\gamma } \nabla _{\gamma }K^\beta -\frac{1}{2}K^{\beta }K_{\beta }\nonumber \\&\quad -2e^{-A} F_1^\alpha \nabla _{1}K_{\alpha }+2e^{-A}F_1^\alpha \nabla _{\alpha }K_1+K^1K_1\Big ]\nonumber \\&\quad +\frac{b^2}{2\pi }\Big [-\frac{1}{4}\nabla ^{\beta } F^{\alpha \gamma }\nabla _{\beta }F_{\alpha \gamma }+F^{\gamma {1}}\nabla ^\beta \nabla _{\beta }F_{1\gamma }\nonumber \\&\quad +F_\gamma ^{1}\nabla _\beta \nabla _1{F}^{\beta \gamma } -\nabla _\beta \Big (F_{\gamma }^\beta \nabla _1{F}_1^\gamma \Big )\Big ] \end{aligned}$$
(70)
$$\begin{aligned} T^2_{\;2}&= R^2E^{(2)}-R^2sin^2\theta (B^{(3)})^2-e^A(B^{(1)})^2\nonumber \\&\quad +\frac{1}{4}F_{\alpha \beta }F^{\alpha \beta } +\frac{a^2}{4\pi }\Big [F_\beta ^{\;\gamma }\nabla _{\gamma }K^\beta \nonumber \\&\quad +\frac{1}{2}K^\beta _\beta -F^{2\alpha }\nabla _2K_\alpha -2F^{2\gamma }\nabla _{\gamma }K_2-K^2K_2\Big ]\nonumber \\&\quad +\frac{b^2}{2\pi }\Big [-\frac{1}{4}\nabla ^{\beta } F^{\alpha \gamma }\nabla _{\beta }F_{\alpha \gamma }\nonumber \\&\quad +F^\gamma _2\nabla ^\beta \nabla _\beta {F}^2_\gamma +F_{\gamma {2}}\nabla _\beta \nabla ^2F^{\beta \gamma } -\nabla _\beta \Big (F_\gamma ^\beta \nabla _{2}F_{2}^\gamma \Big )\Big ] \end{aligned}$$
(71)
$$\begin{aligned} T^3_{\;3}&=\frac{1}{4}g^{33}\Big [-R^4E^{(2)}+R^4sin^2\theta (B^{(3)})^2\nonumber \\&\quad +e^AR^2(B^{(1)})^2-\frac{R^2}{4}F^{\alpha \beta }F_{\alpha \beta }\Big ]\nonumber \\&\quad +\frac{a^2}{4\pi }\Big [g_{33}F_\beta ^{\;\gamma }\nabla _{\gamma }K^\beta \nonumber \\&\quad +\frac{1}{2}g_{33}K^\beta _\beta -F^{3\alpha }\nabla _3K_\alpha -2F^{3\gamma }\nabla _{\gamma }K_3-K^3K_3\Big ]\nonumber \\&\quad +\frac{b^2}{2\pi }\Big [-\frac{1}{4}\nabla ^{\beta }F^{\alpha \gamma } \nabla _{\beta }F_{\alpha \gamma }\nonumber \\&\quad +F^\gamma _3\nabla ^\beta \nabla _\beta {F}^3_\gamma +F_{\gamma {3}}\nabla _\beta \nabla ^3F^{\beta \gamma }-\nabla _\beta \Big (F_\gamma ^\beta \nabla _{3}F_{3}^\gamma \Big )\Big ] \end{aligned}$$
(72)
$$\begin{aligned} T^1_{\;0}&=e^{A/2}R^2sin\theta \Big (E^{(3)}B^{2)}-E^{(2)}B^{(3)}\Big )\nonumber \\&\quad +\frac{a^2}{4\pi }\Big [2F_{(1}^\alpha \nabla _{0)}K_\alpha -2F_{(1}^\alpha \nabla _{\alpha }K_{0)}\Big ]\nonumber \\&\quad +\frac{b^2}{2\pi }\big [F^\gamma _{(1}\nabla _ \alpha \nabla _ {\beta }F_{0)\gamma }F^{\beta \alpha }\nonumber \\&\quad -\nabla _{\beta }F_\gamma ^\beta \nabla _{(1}F_{0)}^\gamma +\nabla _\beta {F}_\gamma ^\beta \nabla _{(1}F_{0)}^\gamma -F_\gamma ^\beta \nabla _\beta \nabla _{(1}F_{0)}^\gamma \Big ] \end{aligned}$$
(73)
$$\begin{aligned} T^2_{\;0}&=g^{22}\Big \{e^{A}R^2sin\theta \Big (E^{(3)}B^{(2)} -E^{(2)}B^{(3)}\Big )\nonumber \\&\quad +\frac{a^2}{4\pi }\Big [2F_{(2}^\alpha \nabla _{0)} K_\alpha -2F_{(2}^\alpha \nabla _{\alpha }K_{0)}\Big ]\nonumber \\&\quad +\frac{b^2}{2\pi }\big [F^\gamma _{(2}\nabla _\alpha \nabla _{\beta } F_{0)\gamma }F^{\beta \alpha }\nonumber \\&\quad -\nabla _{\beta }F_\gamma ^\beta \nabla _{(2}F_{0)}^\gamma +\nabla _\beta {F}_\gamma ^\beta \nabla _{(2}F_{0)}^\gamma -F_\gamma ^\beta \nabla _\beta \nabla _{(2}F_{0)}^\gamma \Big ]\Big \} \end{aligned}$$
(74)
$$\begin{aligned} T^3_{\;0}&=g^{33}\Big \{e^{A}R^2sin\theta \Big (E^{(3)}B^{(2)} -E^{(2)}B^{(3)}\Big )\nonumber \\&\quad +\frac{a^2}{4\pi }\Big [2F_{(3}^\alpha \nabla _{0)} K_\alpha -2F_{(3}^\alpha \nabla _{\alpha }K_{0)}\Big ]\nonumber \\&\quad +\frac{b^2}{2\pi }\big [F^\gamma _{(3}\nabla _\alpha \nabla _{\beta }F_{0)\gamma } F^{\beta \alpha }\nonumber \\&\quad -\nabla _{\beta }F_\gamma ^\beta \nabla _{(3}F_{0)}^\gamma +\nabla _\beta {F}_\gamma ^\beta \nabla _{(3}F_{0)}^\gamma -F_\gamma ^\beta \nabla _\beta \nabla _{(3}F_{0)}^\gamma \Big ]\Big \} \end{aligned}$$
(75)
$$\begin{aligned} T^1_{\;2}&=g^{11}\Big \{e^{A}R^2sin\theta \Big (E^{(3)}B^{(2)}-E^{(2)}B^{(3)} \Big )\nonumber \\&\quad +\frac{a^2}{4\pi }\Big [2F_{(1}^\alpha \nabla _{2)}K_\alpha -2F_{(1}^\alpha \nabla _{\alpha }K_{2)}\Big ]\nonumber \\&\quad +\frac{b^2}{2\pi }\big [F^\gamma _{(1}\nabla _\alpha \nabla _{\beta }F_{2) \gamma }F^{\beta \alpha }\nonumber \\&\quad -\nabla _{\beta }F_\gamma ^\beta \nabla _{(1}F_{2)}^\gamma +\nabla _\beta {F}_\gamma ^\beta \nabla _{(1}F_{2)}^\gamma -F_\gamma ^\beta \nabla _\beta \nabla _{(1}F_{2)}^\gamma \Big ]\Big \} \end{aligned}$$
(76)
$$\begin{aligned} T^1_{\;3}&=g^{11}\Big \{e^{A}R^2sin\theta \Big (E^{(3)}B^{(2)}-E^{(2)}B^{(3)}\Big )\nonumber \\&\quad +\frac{a^2}{4\pi }\Big [2F_{(1}^\alpha \nabla _{3)}K_\alpha -2F_{(1}^\alpha \nabla _{\alpha }K_{3)}\Big ]\nonumber \\&\quad +\frac{b^2}{2\pi }\big [F^\gamma _{(1}\nabla _\alpha \nabla _{\beta } F_{3)\gamma }F^{\beta \alpha }\nonumber \\&\quad -\nabla _{\beta }F_\gamma ^\beta \nabla _{(1}F_{3)} ^\gamma +\nabla _\beta {F}_\gamma ^\beta \nabla _{(1}F_{3)}^\gamma -F_\gamma ^\beta \nabla _\beta \nabla _{(1}F_{3)}^\gamma \Big ]\Big \} \end{aligned}$$
(77)
$$\begin{aligned} T^2_{\;3}&=g^{22}\Big \{e^{A}R^2sin\theta \Big (E^{(3)}B^{(2)} -E^{(2)}B^{(3)}\Big )\nonumber \\&\quad +\frac{a^2}{4\pi }\Big [2F_{(2}^\alpha \nabla _{3)} K_\alpha -2F_{(2}^\alpha \nabla _{\alpha }K_{3)}\Big ]\nonumber \\&\quad +\frac{b^2}{2\pi }\big [F^\gamma _{(2}\nabla _\alpha \nabla _{\beta }F_{3) \gamma }F^{\beta \alpha }\nonumber \\&\quad -\nabla _{\beta }F_\gamma ^\beta \nabla _{(2}F_{3)}^\gamma +\nabla _\beta {F}_\gamma ^\beta \nabla _{(2}F_{3)}^\gamma -F_\gamma ^\beta \nabla _\beta \nabla _{(2}F_{3)}^\gamma \Big ]\Big \} \end{aligned}$$
(78)

However, by Eq. (7), the only non-null off diagonal component of Einstein tensor is \(G^1_0\). Therefore, in order to have

$$\begin{aligned} T_{02}=T_{03}=T_{12}=T_{13}=T_{23}=0 \end{aligned}$$
(79)

we need to impose the conditions

$$\begin{aligned} F_{02}=F_{03}=F_{12}=F_{13}=F_{23}=0. \end{aligned}$$
(80)

After imposing these conditions into Eqs.  (69)–(78) we have that

$$\begin{aligned} T^0_{\;0}&=\frac{1}{2}e^A(E^{(1)})^2+\frac{a^2}{4\pi } \Big [F^{01}H_{01}-\frac{1}{2}K^0K_0+\frac{1}{2}K^1K_1\Big ]\nonumber \\&\quad +\frac{b^2}{2\pi }\Big [-\frac{1}{4}\nabla ^{\beta }F^{\alpha \gamma } \nabla _{\beta }F_{\alpha \gamma }+F^{01}\nabla ^0\nabla _0F_{01}\nonumber \\&\quad +F^{01}\nabla ^1\nabla _1F_{01}-\nabla ^0F_{10}\nabla _0F^{01})\Big ] \end{aligned}$$
(81)
$$\begin{aligned} T^1_{\;1}&=\frac{1}{2}e^A(E^{(1)})^2+\frac{a^2}{4\pi } \Big [F^{01}H_{01} +\frac{1}{2}K^0K_0-\frac{1}{2}K^1K_1\Big ]\nonumber \\&\quad +\frac{b^2}{2\pi }\Big [-\frac{1}{4}\nabla ^{\beta } F^{\alpha \gamma }\nabla _{\beta }F_{\alpha \gamma }+F^{01}\nabla ^0\nabla _0F_{01}\nonumber \\&\quad +F^{01}\nabla ^1\nabla _1F_{01}-\nabla ^0F_{10}\nabla _0F^{01})\Big ] \end{aligned}$$
(82)
$$\begin{aligned} T^2_{\;2}=T^3_{\;3}&=\Big [\frac{1}{2}F_{01}F^{01}\Big ] +\frac{a^2}{4\pi }\Big [F_0^{\;1}\nabla _1K^0\nonumber \\&\quad +F_1^{\;0}\nabla _0K^1+K^0K_0+K^1K_1\Big ]\nonumber \\&\quad +\frac{b^2}{2\pi }\Big [-\frac{1}{4}\nabla ^{\beta }F^{\alpha \gamma } \nabla _{\beta }F_{\alpha \gamma }\Big ] \end{aligned}$$
(83)
$$\begin{aligned} T_{10}&=-\Big (\frac{a^2+2b^2}{4\pi }\Big )K_1K_0 \end{aligned}$$
(84)

So, together with Eq. (7), we have

$$\begin{aligned} G^0_{\;0}&=\Big (-2\frac{R''}{R}+\frac{R'A'}{R}-\frac{R'^2}{R}\Big )e^{-A}+\Big (\frac{{\dot{R}}{\dot{A}}}{R}+\frac{{\dot{R}}^2}{R^2}+\frac{1}{R^2}\Big )\nonumber \\&=8{\pi }G\Big \{{\kappa \rho }+\frac{1}{2}e^A(E^{(1)})^2\nonumber \\&\quad +\frac{a^2}{4\pi }\Big [F^{01}H_{01}-\frac{1}{2}K^0K_0+\frac{1}{2}K^1K_1\Big ]\nonumber \\&\quad +\frac{b^2}{2\pi }\Big [-\frac{1}{4}\nabla ^{\beta } F^{\alpha \gamma }\nabla _{\beta }F_{\alpha \gamma }\nonumber \\&\quad +F^{01}\nabla ^0\nabla _0F_{01}+F^{01}\nabla ^1\nabla _1F_{01}-\nabla ^0F_{10}\nabla _0F^{01})\Big ]\Big \}\nonumber \\ \end{aligned}$$
(85)
$$\begin{aligned} G^1_{\;1}&=-\frac{{\dot{R}}^2}{R^2}+2\frac{{\ddot{R}}}{R}-e^{-A}\frac{R'^2}{R^2}+\frac{1}{R^2},\nonumber \\&=8{\pi }G\Big \{+\frac{1}{2}e^A(E^{(1)})^2+\frac{a^2}{4\pi }\Big [F^{01}H_{01}\nonumber \\&\quad +\frac{1}{2}K^0K_0-\frac{1}{2}K^1K_1\Big ] +\frac{b^2}{2\pi }\Big [-\frac{1}{4}\nabla ^{\beta } F^{\alpha \gamma }\nabla _{\beta }F_{\alpha \gamma }\nonumber \\&\quad +F^{01}\nabla ^0\nabla _0F_{01}+F^{01}\nabla ^1\nabla _1F_{01}-\nabla ^0F_{10}\nabla _0F^{01})\Big ]\Big \} \end{aligned}$$
(86)
$$\begin{aligned} G^2_{\;2}&= G^3_{\;3}=\Big (-\frac{R''}{R}+\frac{A'R'}{2R}\Big )e^{-A}+\frac{{\ddot{R}}}{R}\nonumber \\&\quad +\frac{{\ddot{A}}}{2}+\frac{{\dot{A}}^2}{4}-\frac{{\dot{R}}{\dot{A}}R}{2}\nonumber \\&=8{\pi }G\Big [+\frac{1}{2}F_{01}F^{01}\Big ] +\frac{a^2}{4\pi }\Big [F_0^{\;1}\nabla _1K^0\nonumber \\&\quad +F_1^{\;0}\nabla _0K^1+K^0K_0+K^1K_1\Big ]\nonumber \\&\quad +\frac{b^2}{2\pi }\Big [-\frac{1}{4}\nabla ^{\beta }F^{\alpha \gamma }\nabla _{\beta }F_{\alpha \gamma }\Big ] \end{aligned}$$
(87)
$$\begin{aligned} G^1_{\;0}&=e^{-A}\Big (2\frac{{\dot{R}}'}{R}-\frac{{\dot{A}}R'}{R}\Big )=-2G\Big (\frac{a^2+2b^2}{4\pi }\Big )K^1K_0. \end{aligned}$$
(88)

Hence, in order to obtain a more concise notation, let us define

$$\begin{aligned} \Big (\frac{a^2+2b^2}{4\pi }\Big ){\equiv }M^2\,\,. \end{aligned}$$
(89)

Hence, that Eq. (88) can be written as

$$\begin{aligned} e^{-A}\Big (2\frac{{\dot{R}}'}{R}-\frac{{\dot{A}}R'}{R}\Big )=-2GM^2K^1K_0. \end{aligned}$$
(90)

where, together with the definition of covariant derivative in Eq.  (11) we have that

$$\begin{aligned} K_0&=\partial _1\Big (\frac{1}{2}A(t,r)-E(t,r)\Big )\,\,,\nonumber \\ K_1&=\partial _0\Big (\frac{1}{2}A(t,r)-E(t,r)\Big )\,\,. \end{aligned}$$
(91)

Integrating Eq. (90), we obtain

$$\begin{aligned} R'=e^{A/2}\chi (t,r), \end{aligned}$$
(92)

where

$$\begin{aligned} \chi (t,r)=-{\int }2GM^2K^1K_0dt+\alpha (r), \end{aligned}$$
(93)

where \(\alpha \) is the constant of integration and, as we can see, if the M is zero, we obtain the Maxwell case discussed in [21]. Therefore, defining \(\chi \,{\equiv }\,1-k(t,r)\), we obtain

$$\begin{aligned} e^A=\frac{R'^2}{(1-k(t,r))} \end{aligned}$$
(94)

and the line element in Eq. (1) can be write as

$$\begin{aligned} dS^2=dt^2-\frac{R'^2}{(1-k(t,r))}dr^2-R^2\Big (d\theta ^2+sin^2{\theta }d\phi ^2\Big )\,\,. \end{aligned}$$
(95)

The new result here is that the k-term depends on both r and t coordinate.

For the Podolsky–Maxwell Eq. (64), with the conditions in Eqs. (69)–(78), the non-zero components are

$$\begin{aligned} F_{01}-(a^2+2b^2)H_{01}+2b^2S_{01}&=C(t)\frac{e^{A/2}}{R^2}\nonumber \\ F_{10}-(a^2+2b^2)H_{10}+2b^2S_{10}&=D(r)\frac{e^{A/2}}{R^2} \end{aligned}$$
(96)

where clearly,

$$\begin{aligned} C(t)=-D(r) \end{aligned}$$
(97)

and

$$\begin{aligned} F_{01}-(a^2+2b^2)H_{01}+2b^2S_{01}&=C\frac{e^{A/2}}{R^2} \end{aligned}$$
(98)

The line element in Eq. (95) together with Eq. (65) lead us to

$$\begin{aligned} S_{10}&=-F_{01}\bigg [2\frac{{\ddot{R}}}{R}-\frac{{\dot{R}}{\dot{A}}}{R}-\Big (\frac{R'A'}{R}-2\frac{R''}{R}\Big )e^{-A}\bigg ] \end{aligned}$$
(99)

Therefore, using Eq. (95), we rewrite the Einstein equations for the (00) and (11) components as

$$\begin{aligned} G^0_{\;0}&=\frac{{\dot{R}}^2}{R^2}+2\frac{{\dot{R}}\dot{R'}}{RR'}-\frac{{\dot{R}}{\dot{k}}}{R(1-k)}+\frac{k}{R^2}+\frac{kk'R'}{R^3(1-k)}\nonumber \\&=8{\pi }G\Big \{\kappa \rho +\frac{R'^2}{2(1+2k(t,r))}(E^{(1)})^2\nonumber \\&\quad +\frac{a^2}{4\pi }\Big [F^{01}H_{01}-\frac{1}{2}K^0K_0+\frac{1}{2}K^1K_1\Big ]\nonumber \\&\quad +\frac{b^2}{2\pi }\Big [-\frac{1}{4}\nabla ^{\beta } F^{\alpha \gamma }\nabla _{\beta }F_{\alpha \gamma }\nonumber \\&\quad +F^{01}\nabla ^0\nabla _0F_{01}+F^{01}\nabla ^1\nabla _1F_{01}-\nabla ^0F_{10}\nabla _0F^{01})\Big ]\Big \}\nonumber \\ \end{aligned}$$
(100)
$$\begin{aligned} G^1_{\;1}&=-\frac{{\dot{R}}^2}{R^2}+2\frac{{\ddot{R}}}{R}-(1+2k(t,r)){R^2}+\frac{1}{R^2}\nonumber \\&= 8{\pi }G\Big \{\frac{R'^2}{2(1-k(t,r))}(E^{(1)})^2 +\frac{a^2}{4\pi }\Big [F^{01}H_{01}\nonumber \\&\quad +\frac{1}{2}K^0K_0-\frac{1}{2}K^1K_1\Big ]\nonumber \\&\quad \times \frac{b^2}{2\pi }\Big [-\frac{1}{4}\nabla ^{\beta }F^{\alpha \gamma } \nabla _{\beta }F_{\alpha \gamma }+F^{01}\nabla ^0\nabla _0F_{01}\nonumber \\&\quad +F^{01}\nabla ^1\nabla _1F_{01}-\nabla ^0F_{10}\nabla _0F^{01})\Big ]\Big \} \end{aligned}$$
(101)

Defining,

$$\begin{aligned} \epsilon (t,r)\,{\equiv }\,\frac{1}{2}K^1K_1-\frac{1}{2}K^0K_0+F^{01}\nabla ^0\nabla _0F_{01}-F^{01}\nabla ^1\nabla _1F_{01} \end{aligned}$$
(102)

and

$$\begin{aligned} \epsilon '(t,r)\,{\equiv }\,\frac{1}{2}K^1K_1-\frac{1}{2}K^0K_0-F^{01}\nabla ^0\nabla _0F_{01}+F^{01}\nabla ^1\nabla _1F_{01}\,\,, \end{aligned}$$
(103)

we can rewrite these Einstein equations as

$$\begin{aligned}&\frac{{\dot{R}}^2}{R^2}+2\frac{{\dot{R}}\dot{R'}}{RR'}-\frac{{\dot{R}}{\dot{k}}}{R(1-k)}+\frac{k}{R^2}+\frac{kk'R'}{R^3(1-k)}\nonumber \\&\quad =8{\pi }G\bigg ({\kappa \rho }+\frac{R'^2}{2(1-k(t,r))}(E^{(1)})^2+M^2\epsilon '(t,r)\bigg ) \end{aligned}$$
(104)

and

$$\begin{aligned}&\frac{{\dot{R}}^2}{R^2}+2\frac{{\ddot{R}}}{R}+k(t,r){R^2}\nonumber \\&\quad =8{\pi }G\bigg \{\frac{R'^2}{2(1-k(t,r))}(E^{(1)})^2+M^2\epsilon (t,r)\bigg \} \end{aligned}$$
(105)

which lead us to

$$\begin{aligned} {\dot{R}}^3+2R{\dot{R}}{\ddot{R}}=k(t,r){\dot{R}}+\sigma (t,r)\frac{{\dot{R}}}{R^2} \end{aligned}$$
(106)

where

$$\begin{aligned} \sigma (t,r)\,{\equiv }\,4{\pi }Ge^AR^4\Big (E^2+2M^2\epsilon (t,r){e^A}\Big ) \end{aligned}$$
(107)

So, we write \({\dot{R}}\) as

$$\begin{aligned} {\dot{R}}=\frac{1}{\sqrt{R}}{\int }K(t,r)dt+\frac{1}{\sqrt{R}}{\int }\sigma (t,r)dt+\frac{\alpha (r)}{\sqrt{R}} \end{aligned}$$
(108)

substituting this result into Eq. (104) we obtain

$$\begin{aligned}&8{\pi }G\Big \{{\kappa \rho +\frac{R'^2}{2(1-k(t,r))}(E^{(1)})^2-M^2\epsilon '(t,r) }\Big \}\nonumber \\&\quad =\,\frac{2}{RR'}\bigg (\Sigma (t,r)+\Omega (t,r)+\frac{\alpha }{\sqrt{R}}\bigg )\nonumber \\&\ \qquad \bigg [\Sigma '(t,r)+\Omega '(t,r)+\partial _r\bigg (\frac{\alpha }{\sqrt{R}}\bigg )\bigg ]\nonumber \\&\quad +\frac{1}{R^2}\bigg (\Sigma (t,r)+\Omega (t,r) +\frac{\alpha }{\sqrt{R}}\bigg )^2\nonumber \\&\quad +\frac{1}{R}\bigg (\Sigma (t,r)+\Omega (t,r) +\frac{\alpha }{\sqrt{R}}\bigg )\frac{{\dot{k}}}{R(1-k)}\nonumber \\&\quad +\frac{k}{R^2}+\frac{kk'R'}{R^3(1-k)} \end{aligned}$$
(109)

Equations (95), (108) and (109) define the LTB model with the Podolsky electrodynamics contributions. The singularities arising from \(R=0\), \(R'=0\) and \(k=1\). The \(R=0\) singularity can be understood as the Big Bang singularity, and the \(R'=0\) singularity as the shell cross singularity. The last singularity, \(k=1\), comes from the Podolsky contribution.

6 Conclusion

The analysis of the electromagnetic field in an LTB background was explored by other authors in Ref. [21], however they studied only the Maxwell case, which means that the effects of a mass term or the complication of high-derivative electromagnetic term were not investigated until now. The objective of this paper is to fill this gap and to analyze the possible theoretical effects of these rich electrodynamic scenarios in a more realistic inhomogeneous cosmological background.

In this work, we have analyzed the Proca and the higher-derivative Podolsky models embedded in an inhomogeneous LTB background. In the case of Proca model we found a new singularity at \(k(r)=0\). Moreover, the magnetic field must be zero to satisfy the Einstein’s equations. Considering the Proca model, we have determined the scale factor and provided an analysis of the luminosity distance.

In the Podolsky case, we analyzed the function k, which appears in the line element that defines the LTB model where the Podolsky electrodynamics is dependent of both t and r coordinate. This is different from Maxwell and Proca cases.

As a perspective, a more precise analysis of the Podolsky contributions for the cosmological LTB model can be carried out, as the analysis of black holes to generate the singularity at \(R'=0\), for example.