1 Introduction

Since it was proposed, general relativity has been tested and proved to be quite effective to describe phenomena in the solar system and beyond [1,2,3]. Some of the most famous results of general relativity are: the correction of precession of perihelion of mercury [3, 4], the bending of light due to the gravitational field [5] and the existence of gravitational waves [6,7,8,9]. Another important prediction is the existence of black holes; astrophysical objects whose gravitational interaction is so strong that even light can not escape from [2, 10].

Even though it is effective in describing certain phenomena, general relativity presents problems to describe some situations, such as the accelerating expansion of the universe and the galaxy rotation curves [11]. To solve these problems, we must consider modifications in the field equations. This can be done in two ways. The first is considering changes in the matter/energy sector, through the presence of dark matter and dark energy [12, 13]. The other alternative would be to modify the geometric sector of the equations, which results in the known alternative theories of gravity [14].

One way to obtain the field equations of general relativity is through the variational principle using the Einstein–Hilbert action [15]. In the case of alternative theories of gravity, we can make modifications in this action in such a way that we obtain new equations of motion [14]. In 1980, Starobinsky proposed a model that describes an inflation scenario, where he inserts a term of \(R^2\) in the Einstein–Hilbert action [16]. This kind of modification generates what we call the f(R) theory, where the curvature scalar, in the Einstein–Hilbert action, is replaced by an arbitrary function of this scalar. It is also possible to use f(R) models as substitutes for dark matter and dark energy [17]. It is possible to modify the action by inserting functions of the stress-energy tensor, \(T_{\mu \nu }\), as f(R, T) [18,19,20] and \( f(R,T,R_{\mu \nu }T^{\mu \nu })\) [21], where T is the trace of \(T_{\mu \nu }\) and \(R_{\mu \nu }\) is the Ricci tensor.

The curvature scalar is not the only curvature invariant that can be used to construct the Lagrangian density. There are other invariants like \(R^{\mu \nu }R_{\mu \nu }\) and \(R^{\mu \nu \alpha \beta }R_{\mu \nu \alpha \beta }\) and even combinations of them as \(G=R^2-4R^{\mu \nu }R_{\mu \nu }+R^{\mu \nu \alpha \beta }R_{\mu \nu \alpha \beta }\), Gauss–Bonnet invariant, that could be used. Something interesting about the scalar G is that it is a topological invariant in \(3 + 1\) dimensions or lower. So that, if we include the Gauss–Bonnet term in the Einstein–Hilbert action we will not have modifications in the equations of motion. Even if the linear term does not make changes in the equations of motion, it is possible to modify the field equations including nonlinear terms in G, which are the f(G) theories [21,22,23,24,25,26,27,28,29]. This theory has been used to study the late-time accelerating expansion of the universe [23].

Another possibility of an alternative theory of gravity is the Teleparallel Theory (TT), where we consider zero curvature and nonzero torsion of the spacetime [30,31,32]. In this theory, instead of the curvature scalar, we have the presence of the torsion scalar in Lagrangian density and the field equations obtained from the variational principle are equivalent to the Einstein equations. As it has been done in the case of the f(R) theory, it is also possible to generalize the TT theory by replacing the scalar torsion in the Lagrangian density for a general function of \(\mathcal {T}\) and with that we obtain the known \(f(\mathcal {T})\) theory [33,34,35,36,37,38], where we can explain the present cosmic accelerating expansion with no dark energy [34, 36, 37], it’s also possible to include the trace of the stress-energy tensor to build a \(f(\mathcal {T},T)\) function [39,40,41,42,43,44]. In the context of TT there is also an analogue of the Gauss–Bonnet term, \(T_G\), that can be used to construct the \(f(\mathcal {T},T_G)\) theory [45,46,47].

In addition to the need for dark matter and dark energy, a problem that arises in general relativity, in the study of black holes, is the presence of singularities [48,49,50,51], that are points or set of points where the geodesic is interrupted and the physical quantities diverge [52]. It is believed that the problem of singularity arises because the theory is classical and that in a quantum theory of gravity this problem would be solved [48].

As we do not yet possess this complete theory of quantum gravity, an alternative to avoiding the presence of singularities is the study of regular black holes, which has this name because there are no singularities within it [53]. The first regular solution was proposed by James Bardeen, in 1968 [54]. The Bardeen metric has a de Sitter core with an equation of state \(\rho =-p\), where \(\rho \) is the energy density and p is the pressure. This solution also satisfies the weak energy condition and violates de strong energy condition [55, 56]. Since this solution does not satisfy Einstein equations in vacuum, Ayón-Beato and Garcia have shown that the Bardeen solution can be obtained through Einstein equations with a nonlinear magnetic source [55]. It is also possible to show that this metric is a solution for the case with a nonlinear electrical source, but in this case, it is not possible to obtain a closed form for the electromagnetic Lagrangian [56].

After Bardeen’s proposal, several regular solutions were studied in the literature. Most of these were developed in general relativity [57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72], others in f(R) theory [73,74,75] and \(f(\mathcal {T})\) theory [76]. In Einstein’s theory it is also possible to find solutions with rotation [77,78,79,80,81,82,83] and even with several horizons [84], solutions with multihorizons are also studied in alternative theories of gravity. Regular black holes have already been study in Einstein–Gauss–Bonnet gravity, considering a linear term in the f(G) function, for \(4+1\) dimensions [85].

This work is organized as follows. In the Sect. 2 we use the equations of motion to construct the expressions to the gravity theory and to the electromagnetic, considering a spherically symmetric and static source, and we formulated a theorem for the energy conditions in f(G) gravity. In Sect. 3 we first impose the quasi-global coordinate and, through an Ansatz, we rewrite the expressions in terms of a mass function. After that, we show that it is possible to generalize black hole solutions from general relativity to the f(G) theory. In Sect. 4 we construct the regular models, where we generalized the solutions already known from general relativity to f(G) gravity. As some expressions are too much large or there is no way to find an analytical form, we analyze the results graphically. Our conclusions and perspectives are present in Sect. 5. We dedicate the Appendix A to show the analytical forms of some functions as f(G) and L(P).

Through this work we consider natural units, where \(c=\hbar =G=1\), and metric signature \((+,-,-,-)\).

2 Nonlinear electrodynamics coupled with f(G) gravity: equations of motion

The action that describes f(G) gravity coupled with a nonlinear electrodynamics is given by [22]

$$\begin{aligned} S=\int d^4x\sqrt{-g}\left[ R+f(G)+4L(F)\right] , \end{aligned}$$
(1)

where L(F) is the Lagrangian density of the nonlinear electrodynamics, f(G) represents a general function of the Gauss–Bonnet invariant, \(G=R^2-4R^{\mu \nu }R_{\mu \nu }+R^{\mu \nu \alpha \beta }R_{\mu \nu \alpha \beta }\), with R being the curvature scalar, \(R^{\mu \nu }\) is the Ricci tensor and \(R^{\mu \nu \alpha \beta }\) is the Riemann tensor and g is the determinant of the metric \(g_{\mu \nu }\).

To obtain the field equations, we may vary the action (1) over \(g_{\mu \nu }\) that’s results in

$$\begin{aligned}&G_{\mu \nu }+8\left[ R_{\mu \alpha \nu \beta }+R_{\alpha \nu }g_{\beta \mu }-R_{\alpha \beta }g_{\nu \mu }-R_{\mu \nu }g_{\beta \alpha }\right. \nonumber \\&\quad \left. +R_{\mu \beta }g_{\nu \alpha }\right] \nabla ^\alpha \nabla ^\beta f_G+\left( Gf_G-f\right) g_{\mu \nu }=\kappa ^2T_{\mu \nu }, \end{aligned}$$
(2)

where the subscript G denotes the derivation with respect to the Gauss–Bonnet invariant, \(\kappa ^2=8\pi \) and \(T_{\mu \nu }\) is the stress-energy tensor defined as [1]

$$\begin{aligned} T_{\mu \nu }=-\frac{4}{\kappa ^2\sqrt{-g}}\frac{\delta \left( \sqrt{-g}L\right) }{\delta \left( g^{\mu \nu }\right) }. \end{aligned}$$
(3)

For nonlinear electrodynamics, L(F) is a general function of the scalar \(F=F^{\mu \nu }F_{\mu \nu }/4\), where \(F_{\mu \nu }=\partial _\mu A_\nu -\partial _\nu A_\mu \) is the Maxwell–Faraday tensor and \(A_\mu \) the gauge potential. The stress-energy tensor for this theory is given by [74]

$$\begin{aligned} T_{\mu \nu }=\frac{1}{4\pi }\left[ g_{\mu \nu }L(F)-L_FF_{\mu }^{\ \alpha }F_{\nu \alpha }\right] , \end{aligned}$$
(4)

with \(L_F\equiv \partial L(F)/\partial F\). If we vary the action (1) with respect to \(A_\mu \), we obtain the modified Maxwell equations,

$$\begin{aligned} \nabla _\mu \left[ F^{\mu \nu }L_F\right] \equiv \partial _\mu \left[ \sqrt{-g}F^{\mu \nu }L_F\right] . \end{aligned}$$
(5)

We can always recover the Maxwell theory for \(L(F)=F\) and \(L_F=1\).

At least, the trace of (2) is given by

$$\begin{aligned}&R+8\left[ R_{\alpha \beta }+Rg_{\beta \alpha }\right] \nabla ^\alpha \nabla ^\beta f_G+4\left( f-Gf_G\right) \nonumber \\&\quad =2\left[ FL_F-4L\right] . \end{aligned}$$
(6)

To construct regular solutions, we consider spherically symmetric and static spacetime, which is described by the line element

$$\begin{aligned} ds^2=e^{a(r)}dt^2-e^{b(r)}dr^2-r^2\left( d\theta ^2+\sin ^2\theta d\phi ^2\right) , \end{aligned}$$
(7)

where \(a=a(r)\) and \(b=b(r)\) are arbitrary functions of the radial coordinate. In this work we are considering only electrically charged source, so that, we can integrate the modified Maxwell equations (5) for the line element (7) and we find that the only nonzero component of the Maxwell–Faraday tensor is given by

$$\begin{aligned} F^{10}(r)=\frac{q}{r^2}e^{-(a+b)/2}L_F^{-1}, \end{aligned}$$
(8)

with q being an integration constant and may be interpreted as the electric charge of the source. By the Eq. (8), the electric field will be determined since we have \(L_F\), which is found solving the field equations.

If we write (2) in the mixed form, the nonzero components are

$$\begin{aligned}&\frac{e^{-2 b}}{r^2} \bigg \{2 \left( f_G \left( 2 a''-3 a'b'+a'^2\right) +6 b' f_G'-4f_G''\right) \nonumber \\&\quad +e^{b} \left( 2 f_G \left( 2 a''-a' b'+a'^2\right) -b' \left( r-4f_G'\right) -8 f_G''+1\right) \nonumber \\&\quad +e^{2 b} \left( 1-r^2 f\right) \bigg \}=2\left[ L+\frac{q^2}{r^4}L_F^{-1}\right] , \end{aligned}$$
(9)
$$\begin{aligned}&\frac{e^{-2 b} }{r^2}\bigg \{\left( e^{b}-1\right) \left( e^{b}-4 f_G a''\right) +a' \left( 2\left( e^{b}-3\right) f_G b'\right. \nonumber \\&\ \ \ \qquad \left. +4 \left( e^{b}-3\right) f_G'-re^{b}\right) -2 \left( e^{b}-1\right) f_G a'^2-r^2 e^{2 b}f\bigg \}\nonumber \\&\qquad =2\left[ L+\frac{q^2}{r^4}L_F^{-1}\right] , \end{aligned}$$
(10)
$$\begin{aligned}&\quad -\frac{e^{-2 b}}{4 r^2} \bigg \{16 r a' f_G''-8 \left( f_G-r f_G'\right) \left( 2 a''-3a' b'+a'^2\right) \nonumber \\&\quad +e^{b} \left( 2 \left( 8 f_G+r^2\right) a'' +\left( a'-b'\right) \left( \left( 8 f_G+r^2\right) a'+2 r\right) \right) \nonumber \\&\quad +4 r^2 e^{2 b} f\bigg \}=2L. \end{aligned}$$
(11)

The derivative with respect to the radial coordinate is represented by the prime \((')\). The field equation (2) are equivalent to the equations presented in the original work of f(G) gravity (Eq. (8) in [22]) since we have the same set of Eqs. (9)–(11). As we are working with f(G) theory, it’s also important to calculate the curvature invariants, that are given by

$$\begin{aligned} R= & {} e^{-b}\left[ a''+\left( a'-b'\right) \left( \frac{a'}{2}+\frac{2}{r}\right) +\frac{2}{r}\right] -\frac{2}{r^2}, \end{aligned}$$
(12)
$$\begin{aligned} K\equiv & {} R^{\mu \nu \alpha \beta }R_{\mu \nu \alpha \beta }\nonumber \\= & {} \frac{4}{r^4}-\frac{8e^{-b}}{r^4}+\frac{e^{-2 b}}{4 r^2}\bigg [4 r^2 a''^2-2 r^2 a'^3 b'+r^2 a'^4\nonumber \\&+a'^2 \left( r^2 \left( 4a''+b'^2\right) +8\right) -4 r^2 a' a'' b'+8 b'^2+\frac{16}{r^2}\bigg ].\nonumber \\ \end{aligned}$$
(13)

The Gauss–Bonnet invariant is

$$\begin{aligned} G(r)=\frac{2e^{-2 b} }{r^2}\left[ \left( 2a''+a'^2\right) \left( 1-e^{b}\right) +\left( e^{b}-3\right) a' b'\right] . \end{aligned}$$
(14)

Since we have the equations of motion, we need to solve them to find regular black holes solutions. In general relativity, from the condition \(G^{0}_{\ 0}-G^{1}_{\ 1}=0\) in the Einstein tensor, we have that \(a+b=0\). In alternatives theories of gravity, it’s not necessarily the truth, for example, in the f(R) theory it’s is just one of the possibilities [73,74,75]. If we subtract (9) from (10) we get

$$\begin{aligned}&\frac{e^{-2 b} }{r^2}\left\{ \left( 12f_G'+e^b\left( r-4f_G'\right) \right) \left( a'+b'\right) \right. \nonumber \\&\quad \left. +8\left( e^b-1\right) f_G''\right\} =0, \end{aligned}$$
(15)

that is a second-order differential equations in \(f_G\). The solution of this equation is given by

$$\begin{aligned} f_G(r)= & {} c_0\nonumber \\&+\int _1^r \left( c_1 e^{\int _1^{k_3} -\frac{\left( 12-4 e^{b(k_1)}\right) \left( a'(k_1)+b'(k_1)\right) }{8 \left( e^{b(k_1)}-1\right) } \, dk_1}\right. -e^{\int _1^{k_3} \frac{\left( 4 e^{b(k_1)}-12\right) \left( a'(k_1)+b'(k_1)\right) }{8 \left( e^{b(k_1)}-1\right) } \, dk_1} \nonumber \\&\times \left. \int _1^{k_3} \frac{k_2\left( a'(k_2)+b'(k_2)\right) }{8 \left( e^{b(k_2)}-1\right) }e^{b(k_2)+\int _1^{k_2} \frac{\left( 12-4 e^{b(k_1)}\right) \left( a'(k_1)+b'(k_1)\right) }{8 \left( e^{b(k_1)}-1\right) } \, dk_1} \, dk_2\right) \, dk_3. \end{aligned}$$
(16)

We may solve (15) since we know the relation between the functions a and b. In the next section, we will show a way to generalize solutions constructed in general relativity to f(G) theory.

Since we are interested in studying regular black holes, it is very important to analyze the energy conditions associated with these solutions. In order to analyze the energy conditions, we rewrite (8) as the Einstein equations with an effective stress-energy tensor.

$$\begin{aligned} G_{\mu \nu }= & {} \kappa ^2T_{\mu \nu }-8\left[ R_{\mu \alpha \nu \beta }-R_{\alpha \nu }g_{\beta \mu }-R_{\alpha \beta }g_{\nu \mu }\right. \nonumber \\&\times \left. -R_{\mu \nu }g_{\beta \alpha }+R_{\mu \beta }g_{\nu \alpha }\right] \nabla ^\alpha \nabla ^\beta f_G-\left( Gf_G-f\right) g_{\mu \nu }\nonumber \\= & {} \kappa ^2 T^{eff}_{\mu \nu }, \end{aligned}$$
(17)

where \(T_{\mu \nu }^{eff}\) has the components

$$\begin{aligned} T_{\mu \nu }^{eff}=diag\left( \rho ^{eff}(r),-p^{eff}_r(r),-p^{eff}_t(r),-p^{eff}_t(r)\right) , \end{aligned}$$
(18)

where \(\rho ^{eff}(r)\), \(p^{eff}_r(r)\) and \(p^{eff}_t(r)\) are the effective energy density, radial and tangential pressure, respectively. The energy conditions are given by [29]

$$\begin{aligned} SEC(r)= & {} \rho ^{eff}+p^{eff}_r+2p^{eff}_t \ge 0, \end{aligned}$$
(19)
$$\begin{aligned} WEC_{1,2}(r)= & {} NEC_{1,2}(r)=\rho ^{eff}+p^{eff}_{r,t}\ge 0, \end{aligned}$$
(20)
$$\begin{aligned} WEC_{3}(r)= & {} DEC_{1}(r)=\rho ^{eff} \ge 0, \end{aligned}$$
(21)
$$\begin{aligned} DEC_{2,3}(r)= & {} \rho ^{eff}-p^{eff}_{r,t} \ge 0. \end{aligned}$$
(22)

From (7), (17) and (18) we get

$$\begin{aligned} \rho ^{eff}(r)= & {} \frac{e^{-2 b} }{\kappa ^2 r^2}\bigg \{-4 f_G a''+6 b' \left( f_G a'-2 f_G'\right) \nonumber \\&-2f_G a'^2+2 e^{b} \left( f_G \left( 2 a''-a' b'+a'^2\right) \right. \nonumber \\&\times \left. +2 b'f_G'-4 f_G''\right) +8 f_G''\nonumber \\&+2 r^2 \left( F^{10}\right) ^2 L_F e^{a+3 b}+r^2 e^{2b} f+2 r^2 e^{2 b} L\bigg \}, \end{aligned}$$
(23)
$$\begin{aligned} p^{eff}_r(r)= & {} \frac{e^{-2 b} }{\kappa ^2 r^2}\bigg \{4 f_G a''+2 a' \left( f_G \left( a'-3 b'\right) -6f_G'\right) \nonumber \\&-2 e^{b} \left( 2 f_G a''+a' \left( f_G\left( a'-b'\right) -2 f_G'\right) \right) \nonumber \\&-2 r^2 \left( F^{10}\right) ^2 L_Fe^{a+3 b}-r^2e^{2 b} f-2 r^2 e^{2 b} L\bigg \}, \end{aligned}$$
(24)
$$\begin{aligned} p^{eff}_t(r)= & {} \frac{e^{-2 b} }{\kappa ^2 r^2}\bigg \{2 \left( f_G-r f_G'\right) \left( 2a''-3 a' b'+a'^2\right) -4 r a' f_G''\nonumber \\&-2 e^{b} f_G \left( 2 a''-a' b'+a'^2\right) -r^2e^{2 b} \left( f+2 L\right) \bigg \}. \end{aligned}$$
(25)

In the way that the equations of motion are written, the effective energy density and the effective pressures are equal to the components of Einstein tensor. In this sense, if we consider the same metric, the energy conditions in general relativity and f(G) gravity will be the same. In [75] there is a theorem that says that the energy conditions in f(R) are equal to Einstein theory. We can generalize this theorem to f(G) gravity.

Theorem

Given a solution of (9)–(11), described by the functions \(S_1=\left\{ a(r), b(r), f(G), L(F),F^{10}(r)\right\} \), if we have a solution in general relativity described by \(S_2=\left\{ a(r), b(r), \bar{L}(F),\bar{F}^{10}(r)\right\} \), then the energy conditions are identical for \(S_1\) and \(S_2\) since \(T_{\mu \nu }^{eff}\) in (17) is equal to \(T_{\mu \nu }\) in general relativity.

3 New black hole solutions

Imposing the symmetry that we have from the general relativity, \(b=-a\), (15) becomes

$$\begin{aligned} \frac{8 e^{a} \left( e^{a}-1\right) f_G''}{r^2}=0. \end{aligned}$$
(26)

In general, \(e^a\ne 0\) and \(e^a\ne 1\). In this sense, the solution of (26) is given by

$$\begin{aligned} f_G(r)=c_0+c_1r, \end{aligned}$$
(27)

with \(c_0\) and \(c_1\) being integration constants. This type of behavior is already known in f(R), where we have the No-Go theorem. To construct the function f(G), we need calculate G(r) from (14) and invert the function to get \(r=r(G)\) to replace in (27). Integrating (27) in G we get that

$$\begin{aligned} f(G)=c_0 G+c_1\int r(G)dG. \end{aligned}$$
(28)

A way to construct regular black holes solutions is introducing a mass function, M(r), through the component \(g_{00}\) by

$$\begin{aligned} e^a=1-\frac{2M(r)}{r}, \end{aligned}$$
(29)

where M(r) satisfies the conditions \(\lim \nolimits _{r\rightarrow 0}\left[ M(r)/r\right] \rightarrow 0\) and \(\lim \nolimits _{r\rightarrow \infty }M(r)\rightarrow m\), with m being the ADM mass. In terms of M(r), the Gauss–Bonnet invariant is

$$\begin{aligned} G(r)=\frac{16}{r^6} \left[ r^2 M'^2+r M \left( r M''-4 M'\right) +3 M^2\right] . \end{aligned}$$
(30)

Through (27) and (14), is possible write f(G) in terms of the radial coordinate by

$$\begin{aligned} f(r)=\int f_G(r)\frac{dG(r)}{dr}dr. \end{aligned}$$
(31)

We may use mass functions that are already known from general relativity and verify if the solutions are black holes or regular black holes. From the equations of motion (9)–(11), the quantities related to the electromagnetic theory are given by

$$\begin{aligned} L= & {} \frac{1}{2 r^6}\Bigg [r^2 M'' \Bigg (16 c_0 M+r^2 (8 c_1+r)\Bigg )-r^6f \nonumber \\&+16 \Bigg (M-r M'\Bigg ) \Bigg (r \Bigg (c_1 r-c_0 M'\Bigg )+3c_0 M\Bigg )\Bigg ], \end{aligned}$$
(32)
$$\begin{aligned} L_F= & {} 2 q^2 r\Bigg \{r^2 \Bigg (2 (16 c_1+r) M'+16c_1 M'^2-r (8 c_1+r) M''\Bigg )\nonumber \\&+16 c_1 r M \Bigg (r M''-7 M'-2\Bigg )+96 c_1M^2\Bigg \}^{-1}, \end{aligned}$$
(33)
$$\begin{aligned} F^{10}= & {} \frac{1}{2q r^3}\Bigg \{r^2 \Bigg (-r (8 c_1+r) M''+16 c_1 M'^2+2 (16 c_1+r) M'\Bigg )\nonumber \\&+16 c_1 r M \Bigg (r M''-7 M'-2\Bigg )+96 c_1 M^2\Bigg \}. \end{aligned}$$
(34)

Is also necessary to verify the consistency between the functions L and \(L_F\). This is made by

$$\begin{aligned} L_F=\frac{\partial L}{\partial F}=\frac{\partial L}{\partial r}\left( \frac{\partial F}{\partial r}\right) ^{-1}, \end{aligned}$$
(35)

with \(F=-e^{a+b}\left( F^{10}\right) ^2/2\). To electric sources, it’s very difficult to construct an analytical form to L(F), so that is more convenient use another function, that comes from the Legendre transformation [57, 58]

$$\begin{aligned} \mathcal {H}=2 F L_F-L. \end{aligned}$$
(36)

If we define a new tensor \(P^{\mu \nu }\equiv L_F F^{\mu \nu }\), is possible show that \(\mathcal {H}\) is a function of the scalar \(P=P^{\mu \nu }P_{\mu \nu }/4=\left( L_ F\right) ^2F\). Using (8), the scalar P is given by

$$\begin{aligned} P(r)=-\frac{q^2}{2r^4}, \end{aligned}$$
(37)

and, in terms of the mass function, \(\mathcal {H}\) is

$$\begin{aligned} \mathcal {H}= & {} \frac{1}{2 r^6}\Bigg [16 r M \Bigg (r \Bigg (c_1-(c_0+c_1 r) M''\Bigg )+(4 c_0 +7c_1 r) M'\Bigg )\nonumber \\&-2 r^2 M' \Bigg (8 (c_0+c_1 r) M'+r (8 c_1+r)\Bigg )\nonumber \\&-48 M^2 (c_0+2 c_1 r)+r^6 f\Bigg ]. \end{aligned}$$
(38)

From (37) is possible invert the function to r(P) and then we substitute in (38) to find H(P). Using the auxiliary field \(P^{\mu \nu }\), we can write L(P) as

$$\begin{aligned} L(P)=2P\mathcal {H}_P-\mathcal {H}, \end{aligned}$$
(39)

with \(\mathcal {H}_P=\partial \mathcal {H}/\partial P\).

As the expressions for f(G), L are written in terms of the mass function, each different mass function will generate a different nonlinear electrodynamics and f(G) theory. In [86] the authors proposed a different way to work with the nonlinear electrodynamics in terms of an auxiliary field \(X=q\sqrt{-2F}\) and with that, they could find a closed form to L(X). The methods applied in [86] and here are equivalent since we must obtain the same energy conditions.

3.1 Schwarzschild case

For the Schwarzschild case, the mass function is the ADM mass, \(M(r)=m\), and (29) is

$$\begin{aligned} e^a=e^{-b}=1-\frac{2m}{r}. \end{aligned}$$
(40)

The curvature scalar and the Ricci tensor are zero for Schwarzschild. So that, the Gauss–Bonnet invariant is equal to the Kretschmann scalar.

$$\begin{aligned} G(r)=\frac{48m^2}{r^6}. \end{aligned}$$
(41)

Since we have G(r), we are able to find the r(G) and construct the function f(G). Then, substituting r(G) in (28) we get

$$\begin{aligned} f(G)=c_0 G+\frac{6c_1}{5}\left( 48m^2 G^5\right) ^{1/6}. \end{aligned}$$
(42)

Some important to emphasize is that, in the Einstein’s theory, the Schwarzschild is a vacuum solution, however, in the f(G) theory we have a nonzero stress-energy tensor. If we analyze the components of \(T_{\mu \nu }\) we find some kind of anisotropic matter/field, that is very similar to the nonlinear electrodynamics. If we considered the stress-energy tensor for the nonlinear electrodynamics, despite the fact that we don’t have charge in the geometry, we will find expressions to the electric field and the nonlinear Lagrangian,

$$\begin{aligned} L=\frac{16c_1m\left( 5r-18m\right) }{10r^5}, \end{aligned}$$
(43)
$$\begin{aligned} F^{10}=\frac{32c_1m\left( 3m-r\right) }{2qr^3}. \end{aligned}$$
(44)

In the limit of general relativity, \(c_1\rightarrow 0\), these expressions are zero and then the Schwarzschild solution becomes again a vacuum solution.

3.2 Reissner–Nordström-anti-de Sitter case

To work with de Sitter/anti-de Sitter type solutions, we insert a term with the cosmological constant, \(\Lambda g_{\mu \nu }\), in the left side of (2). As we are interested in solutions with charge, we will consider the Reissner–Nordström-anti-de Sitter (RNAdS) solution, that is described by the mass function

$$\begin{aligned} M(r)=m-\frac{q^2}{2r}-\frac{\Lambda r^3}{6}, \end{aligned}$$
(45)

where \(\Lambda \) is the cosmological constant. The component \(g_{00}\) is

$$\begin{aligned} e^{a(r)}=1-\frac{2m}{r}+\frac{q^2}{r^2}+\frac{\Lambda r^2}{3}. \end{aligned}$$
(46)

The curvature and Gauss–Bonnet invariants are

$$\begin{aligned}&R(r)=4\Lambda , \end{aligned}$$
(47)
$$\begin{aligned}&K(r)=\frac{8\left( 21q^4-36mq^2r+18m^2r^2+r^8\Lambda ^2\right) }{3r^8}, \end{aligned}$$
(48)
$$\begin{aligned}&G(r)=\frac{8\left( 15q^4-36mq^2r+18m^2r^2+r^8\Lambda ^2\right) }{3r^8}. \end{aligned}$$
(49)

We can see that this solution presents a singularity at the center of the black hole. From (49) is not possible to construct an analytical form to r(G) and then we cannot find a closed form to f(G). An alternative procedure is calculated f(G) in terms of the radial coordinate, that is

$$\begin{aligned} f(G)= & {} \frac{8}{35 r^8} \Bigg (35c_0 \Bigg (6 m^2 r^2-12 m q^2 r+5 q^4\Bigg )\nonumber \\&+2 c_1 r \Bigg (126 m^2 r^2-245 m q^2 r+100 q^4\Bigg )\Bigg ), \end{aligned}$$
(50)

and with that we can make a parametric plot showing the behavior \(f_G(G)\times G\) and \(f(G)\times G\). In Fig. 1 we show the difference for the linear and nonlinear case of the theory in G.

Fig. 1
figure 1

Difference between the behavior of f(G) for \(mc_1=0\) and \(mc_1=1\), with \(q=0.1\,m\), \(m^2\Lambda =0.02\) and \(c_0=1\)

Full size image

From (32)–(34), the electromagnetic quantities are

$$\begin{aligned} L(r)= & {} \frac{1}{210 r^7}\Bigg [35 \Bigg (8 c_0 \Lambda ^2 r^7-3 q^2 r^3\Bigg )\nonumber \\&-8 c_1 \Bigg (35 \Lambda r^6+600 q^4+42 m r^2 (18 m-5 r)\nonumber \\&+105 q^2 r (3 r-14 m)\Bigg )\Bigg ], \end{aligned}$$
(51)
$$\begin{aligned} L_F(r)= & {} 6 q^2 r^3\Bigg [6 q^2 r^3+8 c_1 \Bigg (3 q^2 r (5 r-21 m)\nonumber \\&-r^2 (r-3 m) \Bigg (12 m+\Lambda r^3\Bigg )+24 q^4\Bigg )\Bigg ]^{-1}, \end{aligned}$$
(52)
$$\begin{aligned} F^{10}(r)= & {} \frac{1}{6 q r^5}\Bigg [6 q^2 r^3+8 c_1 \Bigg (3 q^2 r (5 r-21 m)\nonumber \\&-r^2 (r-3 m) \Bigg (12 m+\Lambda r^3\Bigg )+24 q^4\Bigg )\Bigg ]. \end{aligned}$$
(53)

So that, the coupled between the electromagnetic and gravity theory makes corrections in these functions. If we consider \(c_1=0\) and \(\Lambda =0\) we recovered the Reissner–Nordström solution in general relativity where the electromagnetic theory is linear with

$$\begin{aligned} L(F)=F=-\frac{q^2}{2r^4}\ \ \text{ and }\ \ L_F=1. \end{aligned}$$
(54)

Something interesting to comments is that, if we expand the electric field for points far from the black hole we get

$$\begin{aligned} F^{10}\approx -\frac{4c_1\Lambda r}{3q}+\frac{4c_1m\Lambda }{q}+\frac{1}{r^2}\left( q-\frac{16c_1m}{q}\right) , \end{aligned}$$
(55)

which clearly diverges due to the presence of the cosmological constant. If we have \(\Lambda =0\) the electric field is well behaved. In f(R) theory, due to the coupled between the gravity and the electromagnetic theory, the electric field diverges at the infinity of the radial coordinate [74, 75], that not necessary happen in the f(G) theory.

As is not possible write an analytical form to L(F), we will show numerically the behavior \(L(F)\times F\) and \(L(P)\times P\) (the analytical form of L(P) is written in the Appendix A). In Fig. 2 we show the nonlinear behavior of the electromagnetic theory.

Fig. 2
figure 2

Behavior of \(L(F)\times F\) and \(L(P)\times P\) for \(c_0=1\), \(mc_1=1\), \(m^2\Lambda =0.02\) and \(q=0.4\,m\)

Full size image

4 Regular black hole solutions

4.1 First regular solution

The first regular solution proposed by Bardeen [54] may be generated by the mass function

$$\begin{aligned} M(r)=\frac{mr^3}{\left( r^2+q^2\right) ^{3/2}}. \end{aligned}$$
(56)

Some generalizations that have Bardeen as a particular case have already been proposed [64, 67]. Let’s consider the mass function

$$\begin{aligned} M(r)=\frac{\alpha m r^{n_1}}{\left( r^{n_2} + q^{n_2}\right) ^{n_3}}, \end{aligned}$$
(57)

where the constant \(\alpha \) is a parameter whose unit depends on the constants \(n_1\), \(n_2\) and \(n_3\). The Bardeen solution is recover for \(n_1=3\), \(n_2=2\) and \(n_3=3/2\). To analyze the regular solutions, we will consider the case \(n_1=3\), \(n_2=2\), \(n_3=3/2\) and \(\alpha =1\). Than find the following functions

$$\begin{aligned} e^{a}=e^{-b}=1-\frac{2mr^2}{\sqrt{q^6+r^6}}, \end{aligned}$$
(58)
$$\begin{aligned}&K(r)=\frac{12 m^2}{\Bigg (q^6+r^6\Bigg )^5}\nonumber \\&\quad \times \Bigg (8 q^{24}-20 q^{18} r^6+183 q^{12} r^{12} -28 q^6 r^{18}+4 r^{24}\Bigg ), \end{aligned}$$
(59)
$$\begin{aligned} G(r)= & {} \frac{48 m^2 \Bigg (2 q^{12}-9 q^6 r^6+r^{12}\Bigg )}{\Bigg (q^6+r^6\Bigg )^3}, \end{aligned}$$
(60)
$$\begin{aligned} f(G)= & {} \frac{48 c_0m^2 \Bigg (2 q^{12}-9 q^6 r^6+r^{12}\Bigg )}{\Bigg (q^6+r^6\Bigg )^3}\nonumber \\&+4c_1m^2 \left\{ \frac{12 r \Bigg (q^{12}-10 q^6 r^6+r^{12}\Bigg )}{\Bigg (q^6+r^6\Bigg )^3}-\frac{1}{q^5}\Bigg [4 \tan ^{-1}\Bigg (\frac{r}{q}\Bigg )\right. \nonumber \\&+\sqrt{3} \ln \Bigg (\frac{q^2+\sqrt{3} q r+r^2}{q^2-\sqrt{3} q r+r^2}\Bigg )\nonumber \\&\left. -2 \tan ^{-1}\Bigg (\sqrt{3}-\frac{2 r}{q}\Bigg )+2 \tan ^{-1}\Bigg (\frac{2 r}{q}+\sqrt{3}\Bigg )\Bigg ]\right\} . \end{aligned}$$
(61)

The analytical forms of \(f_G(G)\) and f(G) in terms of the Gauss–Bonnet invariant are given by (A3) and A4. With that, we can see that the solution is regular in all points of the spacetime and asymptotically flat. We have the limits \(\lim \nolimits _{r\rightarrow \infty }\left\{ e^{a(r)},K(r)\right\} =\left\{ 1,0\right\} \) and \(\lim \nolimits _{r\rightarrow 0}\left\{ e^{a(r)},K(r)\right\} =\left\{ 1,96~m^{2}/q^6\right\} \). It clearly shows that the solution is regular.

From (60) we can write r(G) and with (27) and (61) we construct the functions f(G) and \(f_G(G)\), whose nonlinear behavior is shown in the Fig. 3. We see that \(f_G\) diverges for small values of G(r), which corresponds to \(r\rightarrow \infty \) in (27). The gravity theory clearly is not general relativity, however, can always be recover to \(c_1=0\).

Fig. 3
figure 3

Graphical representation of the functions f(G) and \(f_G(G)\) with \(c_0=1\), \(mc_1=1\) and \(q=0.4\,m\)

Full size image
Fig. 4
figure 4

Graphical representation of L(F) and L(P) with \(mc_1=1\) and \(q=0.4\,m\)

Full size image

The next step is calculate the electromagnetic quantities, that are written

$$\begin{aligned} F^{10}(r)= & {} \frac{m r }{2 q\Bigg (q^6+r^6\Bigg )^3}\Bigg \{96 c_1 m r^8 \Bigg (r^6-5 q^6\Bigg )\nonumber \\&+8 c_1 \sqrt{q^6+r^6} \Bigg (2 q^{12}+25 q^6 r^6-4 r^{12}\Bigg )\nonumber \\&+27 q^6 r^7 \sqrt{q^6+r^6}\Bigg \}, \end{aligned}$$
(62)
$$\begin{aligned} L(F)= & {} -\frac{24 c_1 m^2 r \Bigg (q^{12}-10 q^6 r^6+r^{12}\Bigg )}{\Bigg (q^6+r^6\Bigg )^3}\nonumber \\&+\frac{2 c_1 m^2}{q^5} \Bigg (\sqrt{3}\ln \Bigg [\frac{q^2+\sqrt{3} qr+r^2}{q^2-\sqrt{3} qr+r^2} \Bigg ]+4 \tan ^{-1}\Bigg (\frac{r}{q}\Bigg )\nonumber \\&-2 \tan ^{-1}\Bigg (\sqrt{3}-\frac{2 r}{q}\Bigg )+2 \tan ^{-1}\Bigg (\frac{2r}{q}+\sqrt{3}\Bigg )\Bigg )\nonumber \\&+\frac{4m c_1 \Bigg (2 q^{12}-23 q^6 r^6+2 r^{12}\Bigg )}{r \Bigg (q^6+r^6\Bigg )^{5/2}}+\frac{3m q^6 \Bigg (2 q^6-7r^6\Bigg )}{2\Bigg (q^6+r^6\Bigg )^{5/2}}.\nonumber \\ \end{aligned}$$
(63)

Since we have (37), (62) and (63), we are able to find the behavior of L(F) and L(P). However, is not possible construct L(F) in an analytical form. So that, we show the nonlinearity of the electromagnetic in Fig. 4, which clearly is not Maxwell [the analytical form of L(P) is written in (A2)]. From 5 we can see some new aspects in relation to general relativity. In the Einstein theory, the electric field and the scalar F of a regular black hole go to zero in the origin and at the infinity of the radial coordinate and the \(F^{10}\) has a maximum at the same point where F has a minimum. Here, due to the coupled with the f(G) gravity, the function \(F^{10}\) presents more than one extremum, where which extremum corresponds to a minimum in F. Expand \(F^{10}\) at the infinity we get

$$\begin{aligned} F^{10}\approx -\frac{16c_1m}{qr^2}+\frac{48c_1m^2}{qr^3}+O\left( \frac{1}{r^4}\right) . \end{aligned}$$
(64)

Therefore, the electric field goes to zero at the infinity as in general relativity and different of the f(R) gravity [74, 75].

Fig. 5
figure 5

Electric field and the scalar F in terms of the radial coordinate for \(mc_1=1\), \(\kappa ^2=8\pi \) and \(q=0.7\,m\)

Full size image

Now we will check if the solution satisfies the energy conditions. To do that, we need first calculate the effective energy density and the effective pressures, that are given by

$$\begin{aligned} \rho ^{eff}(r)= & {} \frac{6 m q^6}{\kappa ^2 \left( q^6+r^6\right) ^{3/2}}, \end{aligned}$$
(65)
$$\begin{aligned} p_r^{eff}(r)= & {} -\frac{6 m q^6}{\kappa ^2 \left( q^6+r^6\right) ^{3/2}}, \end{aligned}$$
(66)
$$\begin{aligned} p_t^{eff}(r)= & {} \frac{3 m q^6 \left( 7 r^6-2 q^6\right) }{\kappa ^2 \left( q^6+r^6\right) ^{5/2}}. \end{aligned}$$
(67)

From the fluid quantities, it is possible to notice that the energy density is always positive, that obeys a equation of state of the type \(p_r^{eff}=-\rho ^{eff}\) and an anisotropic behavior \(p_r^{eff}\ne p_t^{eff}\). From the Fig. 6 we can see that close to the center of the black hole we have the behavior of a isotropic fluid \(p_r^{eff}\approx -p_t^{eff}\) with a de Sitter equation of state \(p^{eff}=-\rho ^{eff}\).

Fig. 6
figure 6

Graphical representation of the fluid quantities for \(mc_1=1\) and \(q=0.4\,m\)

Full size image
Fig. 7
figure 7

Kretschmann scalar and Gauss–Bonnet invariant in terms of the radial coordinate for \(q=0.7\,m\)

Full size image
Fig. 8
figure 8

Graphical representation for the functions \(f_G(G)\) and f(G) in terms of the Gauss–Bonnet invariant with \(q=0.7\,m\), \(c_0=1\), \(mc_1=0,2\) for f(G) and \(mc_1=2\) for \(f_G(G)\)

Full size image

Since we have the fluid quantities is easily find the energy conditions. From (19)–(22), the energy conditions are given by

$$\begin{aligned}&WEC_1(r)=0, \end{aligned}$$
(68)
$$\begin{aligned}&WEC_2(r)=\frac{27 m q^6 r^6}{\kappa ^2 \left( q^6+r^6\right) ^{5/2}}, \end{aligned}$$
(69)
$$\begin{aligned}&DEC_2(r)=2DEC_1(r)=\frac{12 m q^6}{\kappa ^2 \left( q^6+r^6\right) ^{3/2}}, \end{aligned}$$
(70)
$$\begin{aligned}&DEC_3(r)=\frac{3 m q^6 \left( 4 q^6-5 r^6\right) }{\kappa ^2 \left( q^6+r^6\right) ^{5/2}}, \end{aligned}$$
(71)
$$\begin{aligned}&SEC(r)=\frac{6 m q^6 \left( 7 r^6-2 q^6\right) }{\kappa ^2 \left( q^6+r^6\right) ^{5/2}}. \end{aligned}$$
(72)

The strong energy condition is violated for \(r<\root 6 \of {2/7}\left| q\right| \). As SEC is related with the fact that the gravitational interaction is attractive, inside the regular black hole, we have a surface of zero-gravity, \(r=\root 6 \of {2/7}\left| q\right| \), where inside this surface the gravitational interaction is repulsive. Outside the black hole \(DEC_3\) is violated for \(r>\root 6 \of {4/5}\left| q\right| \). This type of behavior is already known from the Bardeen solution [56].

4.2 Hayward-type solution

An important regular black hole solution was proposed by Sean Hayward in [69]. If we consider the mass function (57) for \(n_1=3\), \(n_2=3\), \(n_3=1\) and \(\alpha =1\) we get a mass function that generates a Hayward-type solution. As we did with the first solution, we get

$$\begin{aligned} e^{a}= & {} e^{-b}=1-\frac{2 m r^2}{q^3+r^3}, \end{aligned}$$
(73)
$$\begin{aligned} K(r)= & {} \frac{48 m^2 }{\Bigg (q^3+r^3\Bigg )^6}\Bigg (2 q^{12}-2 q^9 r^3+18 q^6 r^6-4 q^3 r^9+r^{12}\Bigg ),\nonumber \\ \end{aligned}$$
(74)
$$\begin{aligned} G(r)= & {} \frac{48 m^2 \Bigg (2 q^6-6 q^3 r^3+r^6\Bigg )}{\Bigg (q^3+r^3\Bigg )^4}, \end{aligned}$$
(75)
$$\begin{aligned} f(G)= & {} \frac{48 c_0 m^2 \Bigg (2 q^6-6 q^3 r^3+r^6\Bigg )}{\Bigg (q^3+r^3\Bigg )^4}\nonumber \\&-\frac{16 c_1 m^2}{3 q^5} \Bigg (\frac{3 q^2 r \Big (r^9+24 q^6 r^3-2q^9\Bigg )}{\Bigg (q^3+r^3\Bigg )^4}\nonumber \\&-\ln \Bigg (\frac{q^2-q r+r^2}{q^2+2qr+r^2}\Bigg )+2 \sqrt{3} \tan ^{-1}\Bigg (\frac{2 r-q}{\sqrt{3} q}\Bigg )\Bigg ). \end{aligned}$$
(76)
Fig. 9
figure 9

Graphical representation of the electric field and the electromagnetic scalar for \(mc_1=1\) and \(q=0.7\,m\)

Full size image
Fig. 10
figure 10

Behavior of the electromagnetic Lagrangian in terms of the electromagnetic scalar and the auxiliary field P for \(mc_1=1\) and \(q=0.7\,m\)

Full size image

In the Fig. 7, by the behavior of K(r) and G(r), we can see the regularity of the spacetime. For points close to the origin these functions tend to a constant and tend to zero in the infinity of the radial coordinate. From (75) together with (27) and (76), we can prove that the gravity theory is not general relativity. In Fig. 8 we compare the f(G) function for the linear case (\(mc_1=0\)) with the case that generates the Hayward solution in f(G) gravity (with \(mc_1=2\)). From \(f_G(G)\), it’s clearly the nonlinearity of the theory since from the linearity case \(f_G(G)\) must be a constant. The analytical expressions are represented by (A8) and (A10).

Fig. 11
figure 11

Effective energy density and effective pressures for the Hayward solution for \(\kappa =\sqrt{8\pi }\) and \(q=0.7\,m\)

Full size image

From (32) and (34) we have

$$\begin{aligned} L= & {} \frac{2 \root 3 \of {2m} c_1}{3 q^{10/3}} \Bigg (\ln \Bigg (\frac{4\root 3 \of {m2q^4}+4\root 3 \of {mq^2}r +2\root 3 \of {2} r^2}{2\root 3 \of {m^2q^4} - \root 3 \of {4q^2m} r+\root 3 \of {2} r^2}\Bigg )\nonumber \\&+2 \sqrt{3} \tan ^{-1}\Bigg (\root 3 \of {\frac{4}{mq^2}}\frac{r}{\sqrt{3}}-\frac{1}{\sqrt{3}}\Bigg )\Bigg )\nonumber \\&-\frac{m }{3 q^5 r }\Bigg (8 c_1 m r \Bigg (\ln \Bigg (\frac{q^2-q r+r^2}{(q+r)^2}\Bigg )\nonumber \\&+2 \sqrt{3} \tan ^{-1}\Bigg (\frac{q-2r}{\sqrt{3} q}\Bigg )\Bigg )\nonumber \\&- \frac{24 c_1q^2}{\Bigg (q^3+r^3\Bigg )^4}\Bigg (6 q^6 r^5 (4 m-r)-2 q^9 r^2 (m+3 r)\nonumber \\&+m r^{11}+q^{12}+q^3 r^9\Bigg )-9 q^8 r \Bigg (q^6-q^3 r^3-2r^6\Bigg )\Bigg ), \end{aligned}$$
(77)
$$\begin{aligned} F^{10}= & {} \frac{m r }{q \Bigg (q^3+r^3\Bigg )^4}\Bigg (8 c_1 \Bigg (3 q^3 r^5 (2 r-7 m)+2 r^8 (3 m-r)\nonumber \\&+q^9+9 q^6 r^3\Bigg )+9 q^3 r^4 \Bigg (q^3+r^3\Bigg )\Bigg ). \end{aligned}$$
(78)

With \(F^{10}\) is simple obtain F and then we show the behavior of these functions in Fig. 9. We can see that the electric field and the electromagnetic scalar are always regular and have zero value in the origin and in the infinity of the radial coordinate. Since we have F and P we may construct L(F) and L(P), whose the nonlinear behaviors are shown in Fig. 10. Finally, the analytical form of L(P) is given by the Eq. (A11).

At least, in order to analyze the energy conditions, we calculate the effective fluid quantities

$$\begin{aligned} \rho ^{eff}(r)= & {} \frac{6 m q^3}{\kappa ^2 \left( q^3+r^3\right) ^2}, \end{aligned}$$
(79)
$$\begin{aligned} p_r^{eff}(r)= & {} -\frac{6 m q^3}{\kappa ^2 \left( q^3+r^3\right) ^2}, \end{aligned}$$
(80)
$$\begin{aligned} p_t^{eff}(r)= & {} -\frac{6 m q^3 \left( q^3-2 r^3\right) }{\kappa ^2 \left( q^3+r^3\right) ^3}. \end{aligned}$$
(81)

It’s clear that we have the behavior of an anisotropic fluid with \(\rho =-p_r\). In Fig. 11 we can see how the fluid quantities, given by the Eqs. (79)–(81), behave in terms of the radial coordinate. The effective energy density is always positive and it’s interesting to notice that all these quantities tend zero in the infinity and zero in the origin of the radial coordinate. It is also important to note that near the center of the black hole we have an isotropic behavior, \(p_r\approx p_t\). The energy conditions are

$$\begin{aligned}&WEC_1(r)=0, \end{aligned}$$
(82)
$$\begin{aligned}&WEC_2(r)=\frac{18 m q^3 r^3}{\kappa ^2 \left( q^3+r^3\right) ^3}, \end{aligned}$$
(83)
$$\begin{aligned}&DEC_2(r)=2DEC_1(r)=\frac{12 m q^3}{\kappa ^2 \left( q^3+r^3\right) ^2}, \end{aligned}$$
(84)
$$\begin{aligned}&DEC_3(r)=\frac{6 m q^3 \left( 2 q^3-r^3\right) }{\kappa ^2 \left( q^3+r^3\right) ^3}, \end{aligned}$$
(85)
$$\begin{aligned}&SEC(r)=-\frac{12 m q^3 \left( q^3-2 r^3\right) }{\kappa ^2 \left( q^3+r^3\right) ^3}. \end{aligned}$$
(86)

The energy conditions do not depend on the parameters of f(G) gravity, being equal to the energy condition on general relativity in agreement with the theorem present in Sect. 2.

4.3 Culetu solution

In [70, 71] Culetu proposed a regular solution that is described by the mass function

$$\begin{aligned} M(r)=me^{-q^2/(2mr)}. \end{aligned}$$
(87)

This solution behaves asymptotically as Reissner–Nordström for regions far from the black hole and as de Sitter close to the center of the black hole. The Culetu solution has already been generalized for f(R) gravity in [74]. From the curvature scalar, is possible construct an analytical expression for r and then construct the f(R) function to the Culetu solution. We will see that it is not so simple in f(G) gravity. The Kretschmann scalar and the Gauss–Bonnet invariant are given by

$$\begin{aligned} K(r)= & {} \frac{e^{-\frac{q^2}{m r}}}{4 m^2 r^{10}} \Bigg (192 m^4 r^4-192 m^3 q^2 r^3\nonumber \\&+96 m^2 q^4 r^2-16 m q^6 r+q^8\Bigg ), \end{aligned}$$
(88)
$$\begin{aligned} G(r)= & {} \frac{8 e^{-\frac{q^2}{m r}} \Bigg (6 m^2 r^2-6 m q^2 r+q^4\Bigg )}{r^8}. \end{aligned}$$
(89)

The regularity of these functions is shown in the Fig. 12. We can see that the Kretschmann scalar zero far from the black hole and at the center and has a maximum. The Gauss–Bonnet has a minimum and a maximum value and goes to zero in the limits \(r\rightarrow 0\) and \(r\rightarrow \infty \). Different of the curvature scalar, it’s not possible to invert G(r) to obtain an analytical form of r(G). So that, to show the nonlinear behave of the gravitational theory, we use (89) with (27) and f(G) in terms of the radial coordinate, that is given by

Fig. 12
figure 12

Kretschmann scalar and Gauss–Bonnet invariant for the Culetu solution for \(q=0.7\,m\)

Full size image
$$\begin{aligned} f(G)= & {} \frac{e^{-\frac{q^2}{m r}}}{q^{10} r^8}, \Bigg (8 c_0 q^{10} \Bigg (6 m^2 r^2-6 m q^2 r+q^4\Bigg )\nonumber \\&-8 c_1 r \Bigg (144 m^7 r^7+144 m^6 q^2 r^6+72 m^5 q^4 r^5+24 m^4 q^6 r^4\nonumber \\&+6 m^3 q^8r^3-6 m^2 q^{10} r^2+7 m q^{12} r-q^{14}\Bigg )\Bigg ), \end{aligned}$$
(90)

to build a parametric plot, Fig. 13, showing that the gravity theory is not the general relativity.

Fig. 13
figure 13

Comparison between the f(G) no the linear case, \(mc_1=0\), and the nonlinear, \(mc_1=2\), generates by the Culetu solution and the \(f_G(G)\) function for \(c_0=1\) and \(q=0.7\,m\)

Full size image

Substituting the mass function that generates the Culetu solution in (34) and (32), we find that the only nonzero component of the Faraday–Maxwell tensor and the electromagnetic Lagrangian, in terms of the radial coordinate, are

$$\begin{aligned} F^{10}= & {} \frac{e^{-\frac{q^2}{m r}}}{8 m q r^5} \Bigg (32 c_1 m \Bigg (12 m^2 r^2-9 m q^2 r+q^4\Bigg )\nonumber \\&+r e^{\frac{q^2}{2 m r}} \Bigg (8 m r^2 \Bigg (q^2-16 c_1 m\Bigg )-q^2 r \Bigg (q^2-96c_1 m\Bigg )8 c_1 q^4\Bigg )\Bigg ), \end{aligned}$$
(91)
$$\begin{aligned} L= & {} \frac{e^{-\frac{q^2}{m r}}}{8 r^7} \Bigg (\frac{r e^{\frac{q^2}{2 m r}}}{m} \Bigg (8 c_1 \Bigg (8 m^2 r^2-8 m q^2 r+q^4\Bigg )+q^2 r \Bigg (q^2-4 m r\Bigg )\Bigg )\nonumber \\&+\frac{32 c_1}{q^{10}}\Bigg (144 m^7 r^7+144 m^6 q^2 r^6+72 m^5 q^4 r^5+24 m^4 q^6 r^4+6 m^3 q^8 r^3\nonumber \\&-6 m^2 q^{10} r^2+7 m q^{12} r-q^{14}\Bigg )\Bigg ). \end{aligned}$$
(92)

If we expand the electric field for larges values of the radial coordinate we find

$$\begin{aligned} F^{10}\approx \frac{4c_1}{qr^2} \left( \frac{12 m^2}{r}-4 m+\frac{5 q^2}{r}\right) -\frac{5 q^3}{8 m r^3}+\frac{q}{r^2}+O\left( \frac{1}{r^4}\right) , \end{aligned}$$
(93)

which clearly is regular for \(r\rightarrow \infty \). It’s interesting since in the f(R) gravity, due to the coupled with the gravitational theory, the electric field diverges for \(r\rightarrow \infty \) (see equation A5 in [74]). From (92), we show the nonlinear behavior of the electromagnetic theory in Fig. 14 where we analyze L in terms of the scalars F and P. The analytical expression for L(P) is given by (A12).

Fig. 14
figure 14

Graphical representation of the electromagnetic Lagrangian in terms of the electromagnetic scalar and the auxiliary field for \(mc_1=1\) and \(q=0.7\,m\)

Full size image

At least, in order to find the energy conditions, we calculate the fluid quantities, that are

$$\begin{aligned}&\rho ^{eff}(r)=\frac{q^2 e^{-\frac{q^2}{2 m r}}}{\kappa ^2 r^4}, \end{aligned}$$
(94)
$$\begin{aligned}&p_r^{eff}(r)=-\frac{q^2 e^{-\frac{q^2}{2 m r}}}{\kappa ^2 r^4}, \end{aligned}$$
(95)
$$\begin{aligned}&p_t^{eff}(r)=\frac{q^2 e^{-\frac{q^2}{2 m r}} \left( 4 m r-q^2\right) }{4 \kappa ^2 m r^5}. \end{aligned}$$
(96)

With that, the energy conditions become

$$\begin{aligned}&WEC_1(r)=0, \end{aligned}$$
(97)
$$\begin{aligned}&WEC_2(r)=\frac{q^2 e^{-\frac{q^2}{2 m r}} \left( 8 m r-q^2\right) }{4 \kappa ^2 m r^5}, \end{aligned}$$
(98)
$$\begin{aligned}&DEC_2(r)=2DEC_1(r)=\frac{2 q^2 e^{-\frac{q^2}{2 m r}}}{\kappa ^2 r^4}, \end{aligned}$$
(99)
$$\begin{aligned}&DEC_3(r)=\frac{q^4 e^{-\frac{q^2}{2 m r}}}{4 \kappa ^2 m r^5}, \end{aligned}$$
(100)
$$\begin{aligned}&SEC(r)=\frac{q^2 e^{-\frac{q^2}{2 m r}} \left( 4 m r-q^2\right) }{2 \kappa ^2 m r^5}. \end{aligned}$$
(101)

So, as in general relativity and in f(R) gravity, SEC and \(NEC_2\) are violated inside the event horizon while the other conditions are always satisfied.

5 Conclusion

In this work, we have developed a method that generalizes solutions of regular black holes, already known from general relativity, to the f(G) theory. The method consists of writing the gravitational and electromagnetic quantities in terms of a mass function, in such a way that each solution generates a different electromagnetic and gravitational theory. Through the equations of motion and considering the symmetry \(a=-b\), we show that the \(f_G(G)\) function has a linear dependence on the radial coordinate, which clearly diverge in the limit \(r\rightarrow \infty \). The divergence in \(f_G(G)\) does not imply that there will also be divergences in f(G). For the models of regular black hole present here f(G) is always regular.

As examples of black holes, we construct the generalization of the Schwarzschild and RNAdS solutions and we showed that in the f(G) gravity the Schwarzschild solution could not be interpreted as a vacuum solution. In the Schwarzschild case, we obtained an analytical form for f(G) and a numerical form to RNAdS. The linear term with the Gauss invariant usually does not modify the equations of motion, however, for RNAdS this term is coupled with the cosmological constant such that to recover the results of general relativity it is necessary to make both \(c_1=0\) and \(c_0 =0\), which is when we recover the linearity of the electromagnetic theory.

For the regular models, we choose a mass function that has the Bardeen and the Hayward solution as specific cases. We also chose the mass function that generates the Culetu solution, so that we can compare this result with those already known from general relativity and the f(R) theory. For the case of coupling with general relativity, the electric field associated with the source is always regular, tends to zero at the origin and at infinity and has a maximum point. For the f(R) gravity, the electric field diverge in infinity and therefore it is necessary to analyze the electric induction tensor. This type of divergence does not appear in the case of coupling with the f(G) theory, but now the electric field has maximum and minimum points.

Since it is not possible to find an analytical form to L(F), we obtain the nonlinear behavior of these functions numerically. As we are working with electrical sources, it becomes useful to use the auxiliary field P and with it, we can find a closed form for L(P).

The energy conditions are the same as those obtained in general relativity and in the f(R). For all solutions, the strong energy condition is violated inside the black hole, which implies that within the event horizon we have regions in which the gravitational interaction is repulsive. For the Culetu solution, we have that the null energy condition is violated outside the event horizon and for the other solutions the dominant energy condition is violated outside the black hole.

As continuations of this work, we can try to verify the stability and the possibility of application of the f(G) models developed here to compare with cosmological data [87]. To do that, since in the Newtonian limit we do not have corrections from the perturbations in f(G) gravity in relation to the general relativity (compare (22) for \(f(R,G)=R\) with (51) in [88]), it’s necessary to use the Post Newtonian and Parameterized Post Newtonian limit to get new corrections. With that, it’s possible to find a regime of validity for \(c_1\), as was done in [89] using the deflection of light, Cassini experiment, perihelion shift retardation of light, gravitational redshift and the equivalence principle.