1 Introduction

General relativity has passed all tests successfully. The latest one was the detection of gravitational wave [1], almost one hundred years after Einstein predicted it. Despite these achievements, it is unavoidable to modify general relativity when spacetime curvature becomes extremely large, say, near a singularity. The most natural modification is to take into account the higher-order curvature terms. The well-known higher-order Lovelock terms provide this kind of modification while they respect the constraints of early version of general relativity [2, 3]. However, these terms have no contribution in four dimensions. Recently, a cubic order curvature model with contribution in four dimensions called Einsteinian cubic gravity (ECG) has been proposed [4]. This model has attracted a lot of attention [5,6,7,8,9,10,11,12,13,14,15,16,17,18,19]. ECG respects some of the constraints that Lovelock theories have. For example, on a maximally symmetric background, it just propagates a transverse and massless graviton and in all dimensions it has the same relative coefficients of the different curvature invariants involved. The nontrivial contribution in four dimensions along with other features have made this model intriguing and important. The latter property allows us to see the effects of higher-order curvature modifications on (\(2+1\))-dimensional holographic duals of gravity theory solutions.

The solutions in the context of ECG have been explored from different points of view. In [5], by constructing perturbative five-dimensional black hole solutions of ECG, the holographic entanglement Rényi entropy has been computed in the dual field theory. The first examples of black hole solutions in ECG have been obtained in [6], and thermal behaviors of them have been explored. In [7], the static and spherically symmetric generalizations of four-dimensional linearly charged and uncharged black hole solutions in ECG have been constructed and their thermodynamics has been studied. The most general theory of gravity to cubic order in curvature called Generalized Quasi-Topological Gravity (GQTG) whose static spherically symmetric vacuum solutions are fully described by a single field equation has been constructed in [8]. In this theory, the ECG as well as Lovelock and quasi-topological gravities have been recovered in four dimensions as special cases. General results corresponding to static and spherically symmetric black hole solutions of general higher-derivative gravities including GQTG have also been established [9]. It has been proved as well that the four-dimensional black hole solutions corresponding to an infinite family of ghost-free higher-order theories are universally stable below a certain mass [10]. In [11], by employing the continued fraction approximation, some interesting properties of ECG black hole solutions such as the innermost stable circular orbit of massive test bodies near a black hole and the shadow of a black hole have been obtained. Some properties of a nonsupersymmetric conformal field theory in three dimensions which could be a holographic dual to four-dimensional ECG model has been explored in [12]. Euclidean AdS-Taub-NUT and bolt solutions with various base spaces in four- and six-dimensions constructed respectively in the context of ECG and GQTG have been studied as well and their thermodynamics features have been explored [13].

In this paper, we study four-dimensional topological black hole solutions of ECG in the presence of nonlinear Born–Infeld (BI) electrodynamics and a bare cosmological constant. To the best of our knowledge, this is the first consideration of nonlinearly charged topological solutions in ECG. The importance of considering nonlinear BI electrodynamics is at least two fold. On the one hand, it resolves the singularity problem of Maxwell electrodynamics at the place of point charge [20]. On the other hand, it comes from the low energy limit of open superstring theory [21,22,23]. In addition, photon-photon interaction experiments have suggested that there is a nonlinear theory of electrodynamics in vacuum [24,25,26,27,28]. Black holes with different horizon’s topologies show drastically different thermodynamical properties as well. For instance, whereas Schwarzschild black holes with spherical horizon are not thermally stable, it has been argued that Schwarzschild-AdS black holes whose horizons have either planar or hyperbolic topologies are thermally stable and do not underlie Hawking-Page phase transition [29]. From holographic point of view, topological solutions are described as duals to thermal states of conformal field theories as well [30].

This paper proceeds as follows. In the next section, we will introduce the action of theory and its field equations. Then, we will calculate conserved quantities in Sect. 3. In Sect. 4, we will compute thermodynamical quantities and check the satisfaction of thermodynamics first law. We will explore thermal stability of our solutions in Sect. 5. Last section is devoted to summary and concluding remarks.

2 Action and field equations

The Einsteinian cubic gravity (ECG), which is the most general dimension-independent gravity theory that consists of metric and Riemann tensor contractions for which linearized spectrum coincides with Einstein gravity one, in the presence of nonlinear electrodynamics could be written as [4, 7]

$$\begin{aligned} {\mathcal {S}}=\frac{1}{16\pi }\int _{{\mathcal {M}}}d^{4}x\sqrt{-g}\bigg ( R+\sum _{i=2}^{3}\alpha _{i}L_{i} -\,2\Lambda _{0}-\lambda {\mathcal {P}}+\mathcal {L }\left( F\right) \bigg ),\nonumber \\ \end{aligned}$$
(1)

up to cubic order in curvature. We write action (1) in Planck units where we set \(G=c=1\). In the action above, \(\Lambda _{0}\), \(\lambda \) and \( \alpha _{i}\)’s are bare cosmological constant, cubic coupling constant and Lovelock coefficients, respectively. We will assume \(\lambda \ge 0\) throughout this paper. Also, \(L_{i}\)’s stand for ith-order Lovelock terms [2, 3]. Note that \(L_2\) is topological and \(L_3\) vanishes identically in four dimensions. The additional cubic contribution \({\mathcal {P}}\) is defined as [4]

$$\begin{aligned} {\mathcal {P}}&= 12R_{a\,b}^{\ c\,d}R_{c\ d}^{\ e\ f}R_{e\ f}^{\ a\, b}+R_{ab}^{cd}R_{cd}^{ef}R_{ef}^{ab} \nonumber \\&\quad -12R_{abcd}R^{ac}R^{bd}+8R_{a}^{b}R_{b}^{c}R_{c}^{a}\,. \end{aligned}$$
(2)

In this paper, we intend to consider Born–Infeld (BI) nonlinear electrodynamics for which

$$\begin{aligned} {\mathcal {L}}\left( F\right) =b^{2}\left( 1-\sqrt{1+\frac{F}{2b^{2}}}\right) , \end{aligned}$$
(3)

where b is a nonlinear parameter and \(F=F_{ab}F^{ab}\) in which \( F_{ab}=2\partial _{[a}A_{b]}\) and \(A_{a}\) is the electromagnetic potential. As b tends to infinity, \({\mathcal {L}}\) reduces to the linear Maxwell case i.e. \(-F/4\). Whereas the second- and third-order Lovelock terms are respectively topological and trivial in four-dimensions, the new cubic term \({\mathcal {P}}\) has contribution in field equations [7].

We make the following ansatz for four-dimensional topological nonlinearly charged black solutions

$$\begin{aligned} ds^{2}= & {} -N^{2}(r)f(r)dt^{2}+\frac{dr^{2}}{f(r)}+r^{2}d\Omega _{k}^{2}\,, \end{aligned}$$
(4)
$$\begin{aligned} A= & {} h\left( r\right) dt, \end{aligned}$$
(5)

where

$$\begin{aligned} d\Omega _{k}^{2}=\left\{ \begin{array}{ll} d\theta ^{2}+\sin ^{2}(\theta )d\phi ^{2}\ \ \ \ \ &{} k=1 \\ d\theta ^{2}+d\phi ^{2} &{} k=0 \\ d\theta ^{2}+\sinh ^{2}(\theta )d\phi ^{2} &{} k=-1 \end{array} \right. , \end{aligned}$$
(6)

represents a 2-dimensional hypersurface with constant curvature 2k and area \(A_{k}\). \(k=1\), 0 and \(-1\) represent two-sphere \(S^{2}\), plane \( {\mathbb {R}} ^{2}\) and hyperbolic \(H^{2}\) topologies for the event horizon, respectively. Varying the action (1) with respect to N(r), f(r) and h(r), we are left with three field equations. One of these field equations which arises from variation with respect to f(r), is satisfied by \(N(r)=const\). Then, we have two field equations

$$\begin{aligned} 0&=-r^{4}{\mathcal {H}}-2r^{2}\left( k-\Lambda _{0}r^{2}-rf^{\prime }-f\right) \nonumber \\&\quad -\,\frac{12\lambda }{r}(r^{3}ff^{\prime \prime 2}+2kr^{2}ff^{\prime \prime \prime }-4krff^{\prime \prime } \nonumber \\&\quad +\,r^{3}ff^{\prime }f^{\prime \prime }f^{\prime \prime \prime 2}+4kff^{\prime }+4rff^{\prime 2}-krf^{\prime 2} \nonumber \\&\quad - \, 4f^{2}f^{\prime }-2rf^{2}\left( rf^{\prime \prime \prime }-2f^{\prime \prime }\right) ), \end{aligned}$$
(7)
$$\begin{aligned} 0&=2h^{\prime 3}-b^{2}(2h^{\prime }+rh^{\prime \prime }), \end{aligned}$$
(8)

where \({\mathcal {H}}=b^{2}-b^{2}\left( 1-h^{\prime 2}/b^{2}\right) ^{-1/2}\), which is reduced to \(-h^{\prime 2}/2\) for Maxwell case. Note that prime denotes the derivative with respect to r. One could immediately solve the electrodynamic field equation (8) as

$$\begin{aligned} h(r)=-\frac{q}{r}\,{\mathbf {F}}\left( \frac{1}{2},\frac{1}{4},\frac{5}{4},- \frac{q^{2}}{b^{2}r^{4}}\right) , \end{aligned}$$
(9)

where \(\mathbf {F}\left( x,y,z,w\right) \) is the Gauss hypergeometric function and q is an integration constant related to total electric charge of black hole. Expanding h(r) about infinity b, the potential of Maxwell electrodynamics can be reproduced as \(h(r)=-q/r+{\mathcal {O}}\left( b^{-1}\right) \). Substituting h(r) from Eq. (9) to Einstein field Eq. (7) and then dividing expressions by \(2r^{2}\), one can integrate Einstein field equation once. After some manipulation, the result could be simplified as follows

$$\begin{aligned}&kr-m-\frac{r^{3}\Lambda _{0}}{3}-rf+\frac{1}{6}b^{2}r^{3} \nonumber \\&\quad \times \left[ 1-\,\mathbf {F}\left( -\frac{1}{2},-\frac{3}{4},\frac{1}{4},- \frac{q^{2}}{b^{2}r^{4}}\right) \right] \nonumber \\&\quad +\frac{\lambda }{r^{2}}\Big [6rff^{\prime \prime }\left( 2k+rf^{\prime }-2f\right) \nonumber \\&\quad -2rf^{\prime 2}\left( 3k+rf^{\prime }\right) -12f^{\prime }f\left( k-f\right) \Big ]=0, \end{aligned}$$
(10)

where m is an integration constant related to total mass of black hole. For BI-(A)dS gravity in four dimensions (\(\lambda =0\)), f(r) could be obtained from above equation as [31]

$$\begin{aligned} f(r)= & {} k-\frac{m}{r}-\frac{r^{2}\Lambda _{0}}{3}\nonumber \\&+ \,\frac{1}{6}b^{2}r^{2}\left[ 1-\,\mathbf {F}\left( -\frac{1}{2},-\frac{3}{4},\frac{1}{4},-\frac{q^{2}}{ b^{2}r^{4}}\right) \right] . \end{aligned}$$
(11)

In this case, the total mass per unit area \(A_{k}\) is [31,32,33,34,35]

$$\begin{aligned} M_{E}=\frac{m}{8\pi }. \end{aligned}$$
(12)

In the next section, we turn to calculate the conserved quantities related to our solutions, namely total mass and total charge.

3 Conserved quantities

In the present section, we intend to compute the total mass and the total charge of our black solutions. In order to find the total mass, we follow the Abbott–Deser–Tekin method [36,37,38] of finding conserved quantities for higher-order curvature gravities [39,40,41]. According to which, we have to find an equivalent quadratic curvature action of the form

$$\begin{aligned} {\mathcal {S}}_{\text {EQCA}}= & {} \int d^{4}x\sqrt{-g}\, [ {\tilde{\kappa }} ^{-1} ( R-2{\tilde{\Lambda }}_{0}) +{\tilde{\alpha }}R^{2} \nonumber \\&+{\tilde{\beta }}R_{ab}R^{ab}+{\tilde{\gamma }}L_{2}] , \end{aligned}$$
(13)

where \(L_{2}=R_{abcd}R^{abcd}-4R_{ab}R^{ab}+R^{2}\) is the second-order Lovelock term known as the Gauss-Bonnet term. Note that equivalent quadratic curvature action has the same vacuum solution and the same linearized field equations as the corresponding higher-order gravity theory. We can write the gravity part of our action (1) which is cubic-order as

$$\begin{aligned}&{\mathcal {S}}_{G}=\int _{{\mathcal {M}}}d^{4}x\sqrt{-g}f\left( R_{cd}^{ab}\right) , \nonumber \\&f\left( R_{cd}^{ab}\right) =\kappa ^{-1}\left[ R+\alpha _{2}L_{2}-2{\Lambda _{0}}-\lambda {\mathcal {P}}\right] , \end{aligned}$$
(14)

in which additional cubic contribution \({\mathcal {P}}\) defined in Eq. (2) can be re-written as

$$\begin{aligned} {\mathcal {P}}&= 12R_{cd}^{ab}R_{af}^{ce}R_{be}^{df}+4R_{cd}^{ab}R_{ab}^{ef}R_{ef}^{cd} \nonumber \\&\quad -12R_{cd}^{ab}R_{a}^{c}R_{b}^{d}+8R_{b}^{a}R_{a}^{c}R_{c}^{b}\,. \end{aligned}$$
(15)

In Appendix A, we show how one can obtain Eq. (15). We suppose that action (14) has an asymptotically (A)dS solution around \( {\bar{g}}_{ab}\) for which

$$\begin{aligned} {\bar{R}}_{cd}^{ab}=\frac{\Lambda }{3}\left( \delta _{c}^{a}\delta _{d}^{b}-\delta _{d}^{a}\delta _{c}^{b}\right) , \end{aligned}$$
(16)

where \(\Lambda \) is the effective cosmological constant. Indeed, for \( \Lambda =0\), we have an asymptotically flat solution. Comparing Eq. (14), with the most general cubic curvature action constructed from Riemann tensor contractions [42]

$$\begin{aligned} {\mathcal {S}}_{\text {GCCA}}= & {} \int d^{4}x\sqrt{-g}\,\left[ \kappa ^{-1}\left( R-2\Lambda _{0}\right) +\alpha R^{2}+\beta R_{ab}R^{ab}\right. \nonumber \\&\left. +\gamma L_{2}+F\left( R_{cd}^{ab}\right) \right] , \nonumber \\ \end{aligned}$$
(17)

in which

$$\begin{aligned} F\left( R_{cd}^{ab}\right)\equiv & {} a_{1}R_{cd}^{ab}R_{af}^{ce}R_{be}^{df}+a_{2}R_{cd}^{ab}R_{ab}^{ef}R_{ef}^{cd}+a_{3}R_{b}^{a}R_{ea}^{cd}R_{cd}^{eb}\\&+a_{4}RR_{cd}^{ab}R_{ab}^{cd}+a_{5}R_{b}^{a}R_{d}^{c}R_{ac}^{bd}+a_{6}R_{b}^{a}R_{a}^{c}R_{c}^{b}\\&+a_{7}RR_{b}^{a}R_{a}^{b}+a_{8}R^{3}, \end{aligned}$$

one can find that

$$\begin{aligned} \gamma= & {} \frac{\alpha _{2}}{\kappa },\,\,a_{1}=-\frac{12\lambda }{\kappa },\nonumber \\ a_{2}= & {} -\frac{4\lambda }{\kappa },\,a_{5}=\frac{12\lambda }{\kappa },\,\,a_{6}=-\frac{8\lambda }{\kappa }, \nonumber \\ \alpha= & {} \beta =a_{3}=a_{4}=a_{7}=a_{8}=0. \end{aligned}$$
(18)

In addition, parameters in actions (13) and (17) are related to each other, so that, in four dimensions [39,40,41]

$$\begin{aligned} \frac{1}{{\tilde{\kappa }}}&\equiv \frac{1}{\kappa }-\frac{\Lambda ^{2}}{3} \left[ a_{1}+4a_{2}+6\left( a_{3}+4a_{4}\right) \right. \nonumber \\&\quad \left. +9\left( a_{5}+a_{6}+4a_{7}+16a_{8}\right) \right] , \end{aligned}$$
(19)
$$\begin{aligned} {\tilde{\Lambda }}_{0}&\equiv \frac{{\tilde{\kappa }}}{\kappa }\Lambda _{0}+\frac{ 2\Lambda }{3}\left( 1-\frac{{\tilde{\kappa }}}{\kappa }\right) , \end{aligned}$$
(20)
$$\begin{aligned} {\tilde{\alpha }}&\equiv \alpha +\frac{\Lambda }{3}\left[ 3a_{1}-6a_{2}-8a_{4}+a_{5}\right. \nonumber \\&\quad \left. +3\left( -a_{3}+2a_{7}+12a_{8}\right) \right] , \end{aligned}$$
(21)
$$\begin{aligned} {\tilde{\beta }}&\equiv \beta +\frac{\Lambda }{3}\left[ -9a_{1}+24a_{2}+16a_{3}+5a_{5}\right. \nonumber \\&\quad \left. +3\left( 16a_{4}+3a_{6}+4a_{7}\right) \right] , \end{aligned}$$
(22)
$$\begin{aligned} {\tilde{\gamma }}&\equiv \gamma +\frac{\Lambda }{3}\left[ -3a_{1}+6a_{2}+3 \left( a_{3}+4a_{4}\right) \right] . \end{aligned}$$
(23)

Using Eq. (18), one obtains

$$\begin{aligned} \frac{\kappa }{{\tilde{\kappa }}}= & {} 1-\frac{8\lambda \Lambda ^{2}}{3}, \end{aligned}$$
(24)
$$\begin{aligned} \Lambda _{0}= & {} {\tilde{\Lambda }}_{0}\left( 1-\frac{8\lambda \Lambda ^{2}}{3} \right) +\frac{16\lambda \Lambda ^{3}}{9}, \end{aligned}$$
(25)
$$\begin{aligned} \kappa {\tilde{\gamma }}= & {} \left( \alpha _{2}+4\lambda \Lambda \right) , \end{aligned}$$
(26)
$$\begin{aligned} {\tilde{\alpha }}= & {} {\tilde{\beta }}=0. \end{aligned}$$
(27)

Moreover, in four dimensions, we have \(\Lambda ={\tilde{\Lambda }}_{0}\) [39,40,41]. Therefore, Eq. (25) leads to

$$\begin{aligned} \frac{8\lambda }{9}\Lambda ^{3}-\Lambda +\Lambda _{0}=0, \end{aligned}$$
(28)

whose solution is the effective cosmological constant \(\Lambda \). As Eq. (28) shows, \(\Lambda =\Lambda _{0}\) for BI-(A)dS gravity (\(\lambda =0 \)). One also needs to calculate the equivalent effective Newton’s constant as

$$\begin{aligned} \frac{1}{{\tilde{\kappa }}_{e}}\equiv \frac{1}{{\tilde{\kappa }}}+8\Lambda \tilde{ \alpha }+\frac{4\Lambda }{3}{\tilde{\beta }}. \end{aligned}$$
(29)

Putting \({\tilde{\alpha }}={\tilde{\beta }}=0\) from Eq. (27) to above relation, we find

$$\begin{aligned} {\tilde{\kappa }}_{e}={\tilde{\kappa }}. \end{aligned}$$
(30)

The total mass per unit area \(A_{k}\) is given by [39,40,41]

$$\begin{aligned} M=\frac{\kappa }{{\tilde{\kappa }}_{e}}M_{E}, \end{aligned}$$
(31)

where \(M_{E}\) is the total mass in Einstein gravity which can be found as \( M_{E}=m/8\pi \) in four dimensions (see Eq. (12)). Thus, using Eqs. (24) and (30), we have

$$\begin{aligned} \frac{\kappa }{{\tilde{\kappa }}_{e}}=\frac{\kappa }{{\tilde{\kappa }}}=1-\frac{ 8\lambda \Lambda ^{2}}{3}, \end{aligned}$$
(32)

which leads to

$$\begin{aligned} M=\left( 1-\frac{8\lambda \Lambda ^{2}}{3}\right) \frac{m}{8\pi }. \end{aligned}$$
(33)

For model to be unitary and free of ghost, it is necessary to have \(\kappa / {\tilde{\kappa }}_{e}>0\) [39,40,41]. Here, it is remarkable to discuss about the restricts imposed by Eq. (28) as well. This equation is cubic in \(\Lambda \) and its discriminant \(\Delta \) is

$$\begin{aligned} \Delta =\frac{32\lambda }{9}\left( 1-6\lambda \Lambda _{0}^{2}\right) . \end{aligned}$$

So, it has three real roots if

$$\begin{aligned} \Delta \ge 0\rightarrow \lambda \left( 1-6\lambda \Lambda _{0}^{2}\right) \ge 0. \end{aligned}$$

One of the possibilities is \(\Lambda _{0}=0\) for which

$$\begin{aligned} \Lambda =0\text { and }\Lambda =\pm \frac{3}{2\sqrt{2\lambda }}. \end{aligned}$$

Note that we have assumed \(\lambda \) to be positive throughout the paper. For \(\Lambda =\pm 3/2\sqrt{2\lambda }\), \(\kappa /{\tilde{\kappa }}_{e}\) given by Eq. (32) is negative. As a result, for vanishing bare cosmological constant \(\Lambda _{0}\), we have just a unitary model with asymptotically flat solution (\(\Lambda =0\)). Another possibility is \(\Lambda _{0}^{2}=1/6\lambda \) for which \(\Delta =0\). In this case, we have

$$\begin{aligned} \text {For }\Lambda _{0}= & {} +\frac{1}{\sqrt{6\lambda }}:\,\Lambda =- \sqrt{\frac{3}{2\lambda }}\text { and }\Lambda =+\frac{1}{2}\sqrt{\frac{3}{ 2\lambda }}, \\ \text {For }\Lambda _{0}= & {} -\frac{1}{\sqrt{6\lambda }}:\,\Lambda =+ \sqrt{\frac{3}{2\lambda }}\text { and }\Lambda =-\frac{1}{2}\sqrt{\frac{3}{ 2\lambda }}, \end{aligned}$$

where the second one in both cases is a double root. One could check that for all these roots \(\kappa /{\tilde{\kappa }}_{e}\le 0\) and therefore the model is not unitary. For some other values of parameters including ones for which \(\Delta <0\), one can find some vacuum solutions within a unitary model as well. Note that for \(\Delta <0\) (\(1-6\lambda \Lambda _{0}^{2}<0\)), the cubic equation (28) has two imaginary complex conjugate roots and just one real root for \(\Lambda \). In our thermodynamical studies which will be presented in next sections, we focus on vanishing bare cosmological constant case (\(\Lambda _{0}=0\)). As discussed above, in this case, an asymptotically flat solution (\(\Lambda =0\)) is just allowed. So, the mass per unit area reads

$$\begin{aligned} M=\frac{m}{8\pi }, \end{aligned}$$
(34)

according to Eq. (33).

Here, we turn to calculate the total charge of our black solutions via Gauss law. Using (9), one can show that \(F_{rt}=h^{\prime }(r)\) is

$$\begin{aligned} F_{rt}=\frac{q}{r^{2}\sqrt{1+q^{2}/b^{2}r^{4}}}. \end{aligned}$$
(35)

Then, using Gauss law, the total charge is given by

$$\begin{aligned} Q=\frac{\,{1}}{4\pi }\int r^{2}{\mathcal {L}}_{F}F_{\mu \nu }n^{\mu }u^{\nu }d{ \Omega }_{k}, \end{aligned}$$
(36)

where \({\mathcal {L}}_{F}=\partial {\mathcal {L}}/\partial F\) and \({\mathcal {L}}\) is the Lagrangian of BI nonlinear electrodynamics expressed in Eq. (3). Also, \(n^{\mu }\) and \(u^{\nu }\) are the unit spacelike and timelike normals to the hypersurface of radius r given as \(n^{\mu }=\left( \sqrt{-g_{tt}} \right) ^{-1}dt=\left( \sqrt{f(r)}\right) ^{-1}dt\) and \(u^{\nu }=\left( \sqrt{g_{rr}}\right) ^{-1}dr=\sqrt{f(r)}dr\). Using Eqs. (35) and (36), we can obtain the total charge of black solutions per unit area as

$$\begin{aligned} Q=\frac{q}{16\pi }. \end{aligned}$$
(37)

In the next section, we will obtain thermodynamical quantities corresponding to our solutions and check the first law of thermodynamics for them.

4 Thermodynamical quantities and Thermodynamics first law

In this section, we intend to calculate thermodynamical quantities and check the first law of thermodynamics. For this purpose, we first have to compute field equations near the black hole horizon. Taylor expansion of metric function \(f\left( r\right) \) near the black hole horizon \(r_{h}\) is

$$\begin{aligned} f(r)=\sum _{n=0}^{\infty }a_{n}(r-r_{h})^{n}\,, \end{aligned}$$
(38)

in which \(a_{n}=f^{(n)}(r_{h})/n!\). Note that \(a_{0}=f(r_{h})=0\) and \( a_{1}=f^{\prime }(r_{h})=2\kappa _{g}\) where \(\kappa _{g}\) is surface gravity on the horizon and \(f^{\prime }(r_{h})\ge 0\). Plugging above expansion into field equation (10), one receives the following equation up to quadratic order of \(r-r_{h}\):

$$\begin{aligned}&kr_{h}-m-8\lambda \kappa _{g}^{2}\left( 2\kappa _{g}+\frac{3k}{r_{h}}\right) -\frac{ \Lambda _{0}}{3}r_{h}^{3} \nonumber \\&\quad +\frac{1}{6}r_{h}^{3}b^{2}\left[ 1-\mathbf {F}\left( -\frac{1}{2},-\frac{3}{ 4},\frac{1}{4},-\frac{q^{2}}{b^{2}r_{h}^{4}}\right) \right] \,+\left[ k-\Lambda _{0}r_{h}^{2}\right. \nonumber \\&\quad \left. -2\kappa _{g}r_{h}-24\lambda k\frac{\kappa _{g}^{2}}{r_{h}^{2}}+\frac{1}{2} r_{h}^{2}b^{2}\left( 1-\sqrt{1+\frac{q^{2}}{b^{2}r_{h}^{4}}}\right) \right] (r-r_{h}) \nonumber \\&\quad +\left[ 72a_{3}\lambda \kappa _{g}(\kappa _{g}+\frac{k}{r_{h}} )+24a_{2}^{2}\lambda \kappa _{g}-a_{2}r_{h}-2\kappa _{g}-\Lambda _{0}r_{h}\right. \nonumber \\&\quad -\frac{24a_{2}\lambda \kappa _{g}}{r_{h}}\left( 4\kappa _{g}+\frac{3k}{ r_{h}}\right) +\frac{1}{2}r_{h}b^{2}\left[ 1-\left( 1+\frac{q^{2}}{ b^{2}r_{h}^{4}}\right) ^{-\frac{1}{2}}\right] \nonumber \\&\quad \left. +\frac{72\lambda \kappa _{g}^{2}}{r_{h}^{3}}k+\frac{96\lambda \kappa _{g}^{3}}{r_{h}^{2}}\right] (r-r_{h})^{2}+{\mathcal {O}}((r-r_{h})^{3})=0. \end{aligned}$$
(39)

Solving the above equation order by order in terms of \(r-r_{h}\) powers, up to second term, we get:

$$\begin{aligned}&kr_{h}-m-8\lambda \kappa _{g}^{2}\left( 2\kappa _{g}+\frac{3k}{r_{h}} \right) -\frac{\Lambda _{0}}{3}r_{h}^{3} \nonumber \\&\quad +\frac{1}{6}b^{2}r_{h}^{3}\left[ 1-\,\mathbf {F}\left( -\frac{1}{2},-\frac{3 }{4},\frac{1}{4},-\frac{q^{2}}{b^{2}r_{h}^{4}}\right) \right] =0, \end{aligned}$$
(40)
$$\begin{aligned}&k+\frac{1}{2}r_{h}^{2}b^{2}\left( 1-\sqrt{1+\frac{q^{2}}{b^{2}r_{h}^{4}}} \right) -\Lambda _{0}r_{h}^{2} \nonumber \\&\quad -2\kappa _{g}r_{h}-24\lambda \frac{\kappa _{g}^{2}}{r_{h}^{2}}k=0. \end{aligned}$$
(41)

As expected, these relations reproduce the following ECG-Maxwell results if \( b\rightarrow \infty \) [7]

$$\begin{aligned}&kr_{h}-m-8\lambda \kappa _{g}^{2}\left( 2\kappa _{g}+\frac{3k}{r_{h}}\right) -\frac{\Lambda _{0}}{3}r_{h}^{3}+\frac{q^{2}}{4r_{h}} =0\,, \end{aligned}$$
(42)
$$\begin{aligned}&k-\frac{q^{2}}{4r_{h}^{2}}-\Lambda _{0}r_{h}^{2}-2\kappa _{g}r_{h}-24\lambda \frac{\kappa _{g}^{2}}{r_{h}^{2}}k =0. \end{aligned}$$
(43)

It is notable to mention that since the third term in (39) is linear in terms of \(a_{3}\), we can compute that in terms of \(a_{2}\). Also, other higher order terms are linear with respect to other coefficients and therefore, all of them could finally be determined in terms of \(a_{2}\). Consequently, the family of solutions have only one free parameter, \(a_{2}\). This fact shows that the present model allows black solutions which have regular horizons where the regularity condition reduces the number of solutions from a two-parameter family to a one-parameter one.

Fig. 1
figure 1

The behavior of T versus \(r_{h}\) for spherical asymptotically flat solutions with various values of \(\lambda \), b and q. Note that the solid (dashed) curves show the temperature in the presence (absence) of cubic terms. In the right panel, dot-dashed curves show the temperature for Einstein–Maxwell regime

Full size image

From Eqs. (40) and (41), we could determine m and surface gravity \( \kappa _{g}\) as functions of model parameters:

$$\begin{aligned} \kappa _{g}&= \frac{r_{h}^{3}}{24\lambda k}\left[ 1+\frac{12k\lambda b^{2}}{ r_{h}^{2}}\left( 1-\sqrt{1+\frac{q^{2}}{b^{2}r_{h}^{4}}}\right) \right. \nonumber \\&\quad \left. -\frac{24k\lambda \Lambda _{0}}{r_{h}^{2}}+\frac{24k^{2}\lambda }{r_{h}^{4} }\right] ^{\frac{1}{2}}-\frac{r_{h}^{3}}{24\lambda k}, \end{aligned}$$
(44)
$$\begin{aligned} \frac{m}{r_{h}}&=k-\frac{24\lambda \kappa _{g}^{2}}{r_{h}^{2}}k-\frac{ \Lambda _{0}r_{h}^{2}}{3}-\frac{16\lambda \kappa _{g}^{3}}{r_{h}} \nonumber \\&\quad +\frac{b^{2}r_{h}^{2}}{6}\left[ 1-\,\mathbf {F}\left( -\frac{1}{2},-\frac{3 }{4},\frac{1}{4},-\frac{q^{2}}{b^{2}r_{h}^{4}}\right) \right] . \end{aligned}$$
(45)

It is notable to point out that m is related to total mass M via Eq. ( 33). If \(\lambda \) and b tends to zero and infinity, respectively, the above equations reduce to RN-(A)ds black hole ones:

$$\begin{aligned} \kappa _{g}&=\frac{k}{2r_{h}}-\frac{q^{2}}{8r_{h}^{3}}-\frac{\Lambda _{0}r_{h}}{2}, \end{aligned}$$
(46)
$$\begin{aligned} \frac{m}{r_{h}}&=k-\frac{\Lambda _{0}r_{h}^{2}}{3}+\frac{q^{2}}{4r_{h}^{2}}. \end{aligned}$$
(47)

The Hawking temperature of our solution can be written in terms of surface gravity [43]

$$\begin{aligned} T=\frac{\kappa _{g}}{2\pi }, \end{aligned}$$
(48)

where \(\kappa _{g}\) has been given in Eq. (44). Using (44) and (48), one could find the temperature as

$$\begin{aligned} T&=\frac{r_{h}^{3}}{48\pi \lambda k}\left[ 1+\frac{12k\lambda b^{2}}{ r_{h}^{2}}\left( 1-\sqrt{1+\frac{q^{2}}{b^{2}r_{h}^{4}}}\right) + \frac{24k^{2}\lambda }{r_{h}^{4}}\right] ^{\frac{1}{2}} \nonumber \\&\quad -\frac{r_{h}^{3}}{48\pi \lambda k}, \end{aligned}$$
(49)

where we set \(\Lambda _{0}=0\). Note that, according to discussions presented in Sect. 3, we study thermodynamics of asymptotically flat solutions (\(\Lambda =0\)) with vanishing bare cosmological constant \( \Lambda _{0}\).

Here, it is worthwhile to discuss about the allowed values of k for asymptotically flat solutions. We know that for Einstein–Maxwell gravity case, the metric function of asymptotically flat topological solution is \(f(r)= k - m/r + q^{2}/r^{2}\), which obviously does not present viable black holes if \(k=0\) and \(k=-1\) since that would imply a wrong sign for the metric for large values of radial coordinate r where \(-m/r\) or \(-1\) would dominate. It is exactly the case for bare ECG because the contributions from the cubic piece are very small, asymptotically (for large r) as well. For the BI charged black holes in ECG, this reasoning holds since the term related to nonlinear electrodynamics in metric function (as appears in Eq. (11)) is asymptotically proportional to \(r^{-2}\). So, we focus on spherically symmetric black hole solutions (\(k=1\)) because those are the only ones which could exist in asymptotically flat space. One also could check that in hyperbolic and planar cases, the temperature given by Eq. (49) is negative and from this point of view, as well, these cases are not physically meaningful. Therefore, we study the thermodynamics of asymptotically flat solutions with spherical topology on horizon.

Now, we are going to discuss the effects of each model parameter \(\lambda \) , b and q on the temperature of spherical asymptotically flat solutions. From Eq. (49), one can see that for \(k=1\), the first term of temperature formula is positive while the second term is negative. Increasing parameters b and q enhances the magnitude of

$$\begin{aligned} \frac{12k\lambda b^{2}}{r_{h}^{2}}\left( 1-\sqrt{1+\frac{q^{2}}{ b^{2}r_{h}^{4}}}\right) , \end{aligned}$$

in Eq. (49). This term is negative for \(k=1\), so, as b or q increases, the temperature T of spherical asymptotically flat solutions decreases. In order to exhibit these effects, we have plotted T versus \(r_{h}\) for different values of b, q and \(\lambda \) in Fig. 1. This figure shows that as each one of parameters b, q or \(\lambda \) grows, the temperature value becomes lower. This means that as the effect of the cubic term grows (or equivalently, as the nonlinearity of electrodynamics weakens), the temperature of spherical black holes decreases. Figure 1 shows that for spherical asymptotically flat solutions, there is a maximum temperature \(T_{\text {max}} \). The behavior of \(T_{\text {max}}\) in terms of model parameters is exhibited in Table 1 as well. One can find that the value of \(T_{ \text {max}}\) decreases as each of the parameters \(\lambda \), b and q increases.

Table 1 The behavior of maximum temperature \(T_{\text {max}}\) for various values of model parameters
Full size table

Let us now calculate the entropy of our black holes. The entropy of black solutions in higher-order gravities could be computed by Wald’s formula given by [44,45,46]

$$\begin{aligned} S=-2\pi \int _{H}d^{2}x\sqrt{h}\frac{\delta {\mathcal {L}}_{G}}{\delta R_{abcd}} \epsilon _{ab}\epsilon _{cd}, \end{aligned}$$
(50)

where h is the determinant of the induced metric on the horizon, \(\delta {\mathcal {L}}_{G}/\delta R_{abcd}\) is the Euler–Lagrange derivative of gravitational Lagrangian and \(\epsilon _{ab}\) is the binormal of the horizon normalized as \(\epsilon _{ab}\epsilon ^{ab}=-2\). Applying the above formula on our theory introduced by action (1), we receive

$$\begin{aligned} S&= \frac{1}{4}\int _{H}d^{2}x\sqrt{h} [1+2\alpha _{2}R_{(2)}+\lambda ( 36R_{b\ d}^{\ e\ f}R_{aecf} \nonumber \\&\quad +3R_{ab}^{\ \ ef}R_{cdef}-12R_{ac}R_{db}-24R^{ef}R_{ebfc}g_{bd} \nonumber \\&\quad +24g_{bd}R_{ce}R_{\ a}^{e})\epsilon ^{ab}\epsilon ^{cd}], \end{aligned}$$
(51)

in which \(R_{(2)}\) is the Ricci scalar of the induced metric on the horizon. This term comes from the Gauss–Bonnet term \( L_{2}=R^{2}-4R_{ab}R^{ab}+R_{abcd}R^{abcd}\) in the action (1). As we have pointed out before, this term is topological and hence has no effect in our calculations so far. However, it contributes to the entropy. Using the metric (4) with \(N=1\), one determines the entropy per unit area as

$$\begin{aligned} S=\frac{r_{h}^{2}}{4}\left[ 1-24\lambda \frac{\kappa _{g}^{2}}{r_{h}^{2}} \left( \frac{2k}{\kappa _{g}r_{h}}+1\right) \right] +k\alpha _{2}, \end{aligned}$$
(52)

in which \(\kappa _{g}\) could be replaced by Eq. (44). It is remarkable to mention that, the above relation reduces to \(r_{h}^{2}/4\) [31] if the higher-order terms disappear (\(\lambda =\alpha _{2}=0\)).

Now, we turn to calculate the electric potential. The electric potential U , measured at infinity with respect to horizon is defined by

$$\begin{aligned} U=A_{\mu }\chi ^{\mu }\left| _{r\rightarrow \infty }-A_{\mu }\chi ^{\mu }\right| _{r=r_{h}}, \end{aligned}$$
(53)

where \(\chi =\partial _{t}\) is the null generator of the horizon. Using (9) and (53), one finds

$$\begin{aligned} U=\frac{q}{r_{h}}\,\mathbf {F}\left( \frac{1}{2},\frac{1}{4},\frac{5}{4},- \frac{q^{2}}{b^{2}r_{h}^{4}}\right) . \end{aligned}$$
(54)
Fig. 2
figure 2

The behaviors of \((\partial ^{2}M/\partial S^{2})_{Q}\) and T (dashed) versus b for different values of \(\lambda \) with \(q=1\) and \(r_{h}=0.5\)

Full size image
Fig. 3
figure 3

The behavior of \(\mathbf {H}_{S,Q}^{M}\) and T (dashed) versus b for different values of \(\lambda \) with \(q=1\) and \(r_{h}=0.5\)

Full size image
Fig. 4
figure 4

The behaviors of \((\partial ^{2}M/\partial S^{2})_{Q}\), \(\mathbf {H} _{S,Q}^{M}\) and T (dashed) versus b for different values of \(\lambda \) with \(q=1\) and \(r_{h}=1\)

Full size image

In order to check the first law of thermodynamics, we first write a Smarr-type formula. Using Eqs. (34), (37) and (45), the Smarr-type formula \(M\left( r_{h},Q\right) \) can be written as

$$\begin{aligned}&M\left( r_{h},Q\right) =\frac{1}{8\pi }\Bigg [kr_{h}-\frac{4\lambda \kappa _{g}^{2}}{r_{h}}\left( 6k+4\kappa _{g}r_{h}\right) \nonumber \\&\quad +\frac{1}{6}b^{2}r_{h}^{3}\left( 1-\,\mathbf {F}\left( -\frac{1}{2},-\frac{3 }{4},\frac{1}{4},-\frac{\left( 16\pi Q\right) ^{2}}{b^{2}A_{k}^{2}r_{h}^{4}} \right) \right) \Bigg ]. \end{aligned}$$
(55)

According to Eq. (52), \(r_{h}=r_{h}(S,Q)\) and in general \( M=M(S,Q) \). We can then consider S and Q as a complete set of extensive quantities for mass. Therefore, temperature T and electric potential U are defined as conjugate intensive quantities for S and Q, respectively and

$$\begin{aligned} T=\left( \frac{\partial M}{\partial S}\right) _{Q},\quad U=\left( \frac{\partial M}{\partial Q}\right) _{S}. \end{aligned}$$

These intensive quantities can be calculated using Eqs. (37), (52) and (55) whereFootnote 1

$$\begin{aligned} T= & {} \left( \frac{\partial M}{\partial S}\right) _{Q}=\left( \frac{\partial M }{\partial r_{h}}\right) _{Q}\left( \frac{\partial S}{\partial r_{h}}\right) _{Q}^{-1}, \end{aligned}$$
(56)
$$\begin{aligned} U= & {} \left( \frac{\partial M}{\partial Q}\right) _{S}=\left( \frac{\partial M }{\partial Q}\right) _{r_{h}}+\left( \frac{\partial M}{\partial r_{h}} \right) _{Q}\left( \frac{\partial r_{h}}{\partial Q}\right) _{S} \nonumber \\= & {} \left( \frac{\partial M}{\partial Q}\right) _{r_{h}}-\frac{\left( \frac{ \partial S}{\partial Q}\right) _{r_{h}}\left( \frac{\partial M}{\partial r_{h}}\right) _{Q}}{\left( \frac{\partial S}{\partial r_{h}}\right) _{Q}} \nonumber \\= & {} \left( \frac{\partial M}{\partial Q}\right) _{r_{h}}-T\left( \frac{ \partial S}{\partial Q}\right) _{r_{h}}. \end{aligned}$$
(57)

Our calculations show that T and U computed by Eqs. (56) and (57) coincide with ones obtained from Eqs. (49) and (54) for nonlinearly charged spherical (\(k=1\)) asymptotically flat solutions. Thus, the first law of thermodynamics

$$\begin{aligned} dM=TdS+UdQ, \end{aligned}$$
(58)

is satisfied by these quantities.

In next section, we will study the thermal stability of spherical solutions in both canonical and grand canonical ensembles.

5 Thermal stability

In this section, we intend to study thermal stability of the four-dimensional nonlinearly charged asymptotically flat solutions of ECG with spherical horizon in both canonical and grand canonical ensembles. In canonical ensemble where the charge is a fixed parameter, the positivity of heat capacity \(C=T/(\partial T/\partial S)_{Q}=T/(\partial ^{2}M/\partial S^{2})_{Q}\) guarantees the local stability [47,48,49]. Therefore, it is sufficient to check that \((\partial ^{2}M/\partial S^{2})_{Q}\) is positive in order to explore the stability of solutions in the ranges where temperature is positive as well. Charge Q is no longer fixed in the grand canonical ensemble. In this ensemble, the system is locally stable provided that the determinant of Hessian matrix \(\mathbf {H}_{S,Q}^{M}=\left[ \partial ^{2}M/\partial S\partial Q\right] \) is positive, in positive temperature ranges. Since the explicit forms of \((\partial ^{2}M/\partial S^{2})_{Q}\) and \(\mathbf {H}_{S,Q}^{M}\) are complicated, we avoid writing them here. However, the main results are depicted in Figs. 2, 3 and 4.

Let us now turn to study the stability of black hole solutions in canonical and grand canonical ensembles. Our investigations show that the values of charge q and horizon radius \(r_{h}\) determine how thermal stability is influenced by the other parameters. We study thermal stability in canonical ensemble for two different types of solutions (Figs. 2, 4a). Let us focus on the first type exhibited in Fig. 2. As Fig. 2a shows, in highly nonlinear electrodynamics regime (small nonlinearity parameter b values), there is a minimum value for cubic coupling \(\lambda \) (\(\lambda _{\min }\)) that for values greater than it, the system is stable up to \(b_{\max }\) (Fig. 2b). For \(\lambda \) values a bit smaller than \(\lambda _{\min }\), the system is stable except for b values between \( b_{0}\) and \(b_{1}\) (Fig. 2a). As \(\lambda \) becomes smaller, \(b_{0}\) tends to zero and \(b_{1}\) becomes greater. So, the instability interval becomes larger. In this case, the system is stable for b values between \( b_{2}\) and \(b_{\max }\) (Fig. 2b). \(b_{2}\) and \(b_{\max }\) increase as \(\lambda \) decreases. Also, for b values greater than \(b_{\max }\) (including linear Maxwell case), the system is unstable. Note that we plot the curves in Fig. 2b just for two \(\lambda \) values in order not to have a messy plot. For other \(\lambda \) values, the qualitative behavior is the same. In addition, we plot the temperature in Figs. 2, 3 and 4 as well, to ensure it is positive in the regions under study. One could see that the behavior of temperature in terms of \(\lambda \) and b coincides with what we have discussed before i.e. temperature decreases as each of the latter parameters increases. In second type displayed in Fig. 4a, we again have a minimum value for \(\lambda \). However, in this case, there is no upper bound for b and the system is stable for all b values if \( \lambda \ge \lambda _{\min }\). For \(\lambda \) values a bit smaller than \( \lambda _{\min }\), there is a minimum value for b (\(b_{\min }\)) as well, that for \(b>b_{\min }\), we have a stable system. As \(\lambda \) becomes smaller, we have a critical value \(\lambda _{c}\) that for \(\lambda <\lambda _{c}\), the system is unstable for all b values.

We continue by discussing the thermal stability in grand canonical ensemble. For the first type represented in Fig. 3, in highly nonlinear electrodynamics regime (Fig. 3a), the system is unstable for some cubic couplings, however for some greater ones, the system becomes stable up to a specific value of b (\(b_{1}^{\prime }\)). For b values greater than \( b_{1}^{\prime }\), we have an unstable system between \(b_{1}^{\prime }\) and \( b_{2}^{\prime }\) (Fig. 3b). For \(b>b_{2}^{\prime }\) (including linear Maxwell case), the system becomes thermally stable. \(b_{1}^{\prime }\) (\(b_{2}^{\prime }\)) is greater (smaller) for greater values of \(\lambda \). In second type displayed in Fig. 4b, whereas the system is unstable for some values of \(\lambda \), it becomes totally stable as \(\lambda \) grows.

6 Summary and concluding remarks

In the present paper, we studied four-dimensional topological Born–Infeld (BI) charged black hole solutions in the context of Einsteinian cubic gravity (ECG) in the presence of a bare cosmological constant. ECG is the most general gravity up to cubic order in curvature which is independent of dimension and its linearized spectrum coincides with general relativity one. Also, both theoretical (open superstring theory) and experimental (photon-photon interaction experiments) evidences suggest nonlinear theories of electrodynamics such as BI model. On the other hand, topological black hole solutions are known to be dual to some thermal states in the holographic settings, with quite different thermodynamical features depending on the specific topology. To the best of our knowledge, this is the first study of topological black hole solutions in ECG with nonlinear electrodynamics.

We first introduced the action of the theory and obtained its corresponding field equations. Integrating these field equations, we obtained the electromagnetic potential as well as a second order differential equation for the metric function. Then, we presented the total mass M formula by employing Abbott–Deser–Tekin method and showed that the total mass depends on cubic coupling in general. We obtain the total charge Q via Gauss formula as well. Moreover, we discussed the conditions under which the model is unitary and perturbatively ghost-free and explored different possible cases. We found that if the bare cosmological constant vanishes, then the model is unitary only if the solution is asymptotically flat. Therefore, we focused on studying the thermodynamics of asymptotically flat solutions. These are necessarily of spherical topology. Expanding the metric function at the near horizon region up to second order in the field equations, we obtained the temperature T and the total mass as functions of horizon radius \(r_{h}\) and charge q. Next, we revealed the influences of cubic coupling \(\lambda \), nonlinearity parameter b and black hole charge on the behavior of the temperature. We found that the temperature decreases as the values of these parameters increase. It is remarkable to mention that, as b grows, the linear Maxwell electrodynamics is reproduced. In addition, the temperature of spherical asymptotically flat solutions have a maximum value as horizon radius changes. This maximum temperature is smaller for greater values of \(\lambda \), b and q. In order to check the satisfaction of thermodynamics first law, we turned to compute the entropy and the electric potential. We used the Wald formula in order to obtain the entropy of our solutions. We showed that the first law of thermodynamics is satisfied for spherical BI charged asymptotically flat solutions of ECG.

Finally, we analysed thermal stability of our solutions in both canonical and grand canonical ensembles. We showed that the values of q and \(r_{h}\) specify how the stability is affected by other parameters, in both ensembles. For one type of solutions, the system is unstable for b values larger than a maximum value \(b_{\max }\) (including linear Maxwell case) in canonical ensemble. In the small b regime (highly nonlinear electrodynamics), there is also a minimum \(\lambda \) value. The system is stable if \(\lambda \ge \lambda _{\min }\) for b in the range 0 to \(b_{\max }\). There exists another type of solution with \(\lambda _{\min }\) defined as before, but in which, black solutions are stable if \(\lambda \ge \lambda _{\min }\) with no upper bound for b. Furthermore, for \(\lambda \) values a bit smaller than \(\lambda _{\min }\), the system is stable for \(b>b_{\min }\). There is also a critical value for \( \lambda \) that if \(\lambda <\lambda _{c}\), the black hole solutions are thermally unstable for all b values. In grand canonical ensemble, the first type system is stable for large b values (including linear Maxwell case). Nevertheless, in small b regime, it is unstable for some \(\lambda \) values. In this regime, the system becomes stable for b values between 0 and an upper bound, as \(\lambda \) grows. The second type system, in grand canonical ensemble, is totally unstable for some \(\lambda \) values. However, it becomes totally stable, as \(\lambda \) increases.

In the present work, we focused on asymptotically flat solutions with vanishing bare cosmological constant according to discussions given in Sect. 3. In such a scenario, the higher curvature terms have no contribution in the total mass formula. However, if one considers possible (A)dS vacuum solutions in a unitary model, a correction due to cubic curvature terms appears in the total mass formula (see Eq. (33)). This fact may affect the thermdoynamical studies, drastically. It is fascinating to explore these kinds of charged solutions in future works. In addition, here, we considered Born–Infeld model as nonlinear electrodynamics. It is interesting to study the effects of other known nonlinear models of electrodynamics on ECG solutions. Moreover, it is worthwhile to explore the critical behavior of nonlinearly charged solutions of ECG in both extended [50] and non-extended [51] thermodynamics phase spaces. Some of these issues are under investigation and the results will appear elsewhere.