1 Introduction

The discovery of the recent accelerated expansion of the Universe [1, 2] has spurred much interest in modified theories of gravities [3,4,5], as an alternative to dark energy [6] and a possible cause of this late-time cosmic speed-up. In fact, a plethora of modified gravity theories have been explored in the literature that involve generalizations of the Hilbert–Einstein action, such as the f(R) gravity [7] and extensions [8,9,10,11,12], and the inclusion of second order curvature invariants such as \(R_{\mu \nu }R^{\mu \nu }\), \(R_{\alpha \beta \mu \nu }R^{\alpha \beta \mu \nu }\), \(C_{\alpha \beta \mu \nu }C^{{\alpha \beta \mu \nu }}\), \(\varepsilon ^{\alpha \beta \mu \nu }R_{\alpha \beta \gamma \delta }R_{\mu \nu }^{\gamma \delta }\), etc [13]. For further details on applications of certain modified gravities, we recommend the following references [14,15,16]. However, this class of Riemannian geometries has been recently extended to include new fundamental blocks to describe spacetime, such as the presence of the torsion scalar T [17] and nonmetricity, \({\mathbb {Q}}\) [18].

A fundamental result in differential geometry states that the general affine connection, which is a mathematical tool that plays a central role in defining the transport of tensors over the manifold and in defining the covariant derivative, may always be decomposed into the following three independent components [19, 20]:

$$\begin{aligned} {\bar{\Gamma }}{}_{\mu \nu }^{\beta }=\Gamma _{\mu \nu }^{\beta } +K_{\mu \nu }^{\beta }+L_{\mu \nu }^{\beta }, \end{aligned}$$
(1)

where the first term corresponds to the Levi-Civita connection

$$\begin{aligned} \Gamma _{\mu \nu }^{\beta } =\frac{1}{2}g^{\beta \alpha }\left( \partial _{\mu }g_{\nu \alpha } +\partial _{\nu }g_{\alpha \mu }-\partial _{\alpha }g_{\mu \nu }\right) , \end{aligned}$$
(2)

while the second term is known as the contorsion

$$\begin{aligned} K_{\mu \nu }^{\beta } =\frac{1}{2}T_{\mu \nu }^{\beta } +T_{(\mu \nu )}^{\beta }, \end{aligned}$$
(3)

and is defined in terms of the torsion tensor

$$\begin{aligned} T_{\mu \nu }^{\beta } =2\varGamma _{\left[ \mu \nu \right] }^{\beta } =-T_{\nu \mu }^{\beta }. \end{aligned}$$
(4)

Finally, we have the deformation tensor

$$\begin{aligned} L_{\mu \nu }^{\beta } =\frac{1}{2}Q_{\mu \nu }^{\beta } -Q_{(\mu \nu )}^{\beta } =L_{\nu \mu }^{\beta }, \end{aligned}$$
(5)

which is defined by the divergence of the metric tensor, i.e. the non-metricity tensor

$$\begin{aligned} Q_{\beta \mu \nu }\equiv \nabla _{\beta }g_{\mu \nu }\ne 0, \end{aligned}$$
(6)

where \(\nabla _\beta \) denotes the covariant derivative, which is given with respect to the affine connection (1).

In fact, by imposing assumptions on the affine connection, one specifies the metric-affine geometry. Thus, in the standard formulation of General Relativity (GR), one assumes a Levi-Civita connection, which implies a vanishing torsion and nonmetricity, while in the Teleparallel Equivalent of GR (TEGR), based on the torsion scalar T, one uses the Weitzenböck connection, implying zero curvature and nonmetricity [17, 21]. Another equivalent formulation of GR is the symmetric teleparallel equivalent of GR (STEGR), where it is the nonmetricity tensor \({\mathbb {Q}}\) that describes the gravitational interaction [18, 22]. These three frameworks have been denoted the “Geometric Trinity of Geometry” [23].

As interesting extension of STEGR, namely, \(f({\mathbb {Q}})\) gravity, is related with the inclusion of a boundary term \({\mathbb {B}}_Q\). It was found that the resulting theory, \(f({\mathbb {Q}},{\mathbb {B}}_Q)\) gravity, is dynamically equivalent to f(R) gravity, with \( f({\mathbb {Q}},{\mathbb {B}}_Q)=f({\mathbb {Q}}-{\mathbb {B}}_Q)\) [24]. In fact, it was also shown that in extended teleparallel gravity f(T), gravity, one also has \(f(T,{\tilde{B}})=f(-T+{\tilde{B}})\) gravity (where \({\mathbb {B}}_Q\ne {\tilde{B}}\)) [25]. Thus, it was suggested [24] that in this perspective, considering boundary terms in \( f({\mathbb {Q}})\) gravity represents the last ingredient towards an “Extended Geometric Trinity of Gravity”, where f(R), \(f(T,{\tilde{B}})\), and \(f({\mathbb {Q}},{\mathbb {B}}_Q)\) can be dealt with under the same standard.

Furthermore, there are several relevant studies on the theory \(f({\mathbb {Q}})\) with cosmological applications that have been investigated in [26,27,28,29,30]. In addition to the studies mentioned earlier on this theory, we also recommend consulting the following references to deepen your understanding [31,32,33,34,35,36,37,38]. A study on FLRW cosmology within the scope of the theory \(f({\mathbb {Q}},{\mathbb {B}}_Q)\), considering several families of connections, was developed in Ref. [39]. The behavior of cosmological models of dark energy described by the theory \(f({\mathbb {Q}},{\mathbb {B}}_Q)\) was investigated with perfect fluid in Ref. [40] and with quintessence in Ref. [41].

It this work, we explore static and spherically symmetric black hole and regular black hole solutions in \(f({\mathbb {Q}},{\mathbb {B}}_Q)\) gravity coupled to nonlinear electrodynamics (NLED). Since the discovery of non-singular black hole solutions [42, 43], a plethora of solutions have been analysed in this class of solutions (rather than provide an exhaustive list of the solutions, we refer the reader to [43,44,45,46,47,48,49,50] and the references therein).

This paper is organised in the following manner: in Sect. 2, we present the fundamentals of theories of gravity with non-metricity. In Sect. 3, we consider \(f({\mathbb {Q}},{\mathbb {B}}_Q)\) gravity coupled to NLED and present the equations of motion of \(f({\mathbb {Q}},{\mathbb {B}}_Q)\) that will be used throughout this work. In Sect. 4, we find several static and spherically symmetric black hole solutions. In Sect. 5, we explore regular black holes, and generalize the regular Bardeen and Culetu solutions. Finally, in Sect. 6, we summarize our results and conclude.

2 Theory of gravity with non-metricity

2.1 Fundamentals

In this section, we present the fundamentals of \(f({\mathbb {Q}},{\mathbb {B}}_Q)\) gravity, where the non-metricity tensor is given by \(Q_{\beta \mu \nu }\equiv \nabla _{\beta }g_{\mu \nu }\ne 0\), as mentioned above. The two characteristics of the non-metricity tensor are defined by their contractions:

$$\begin{aligned} Q_{\alpha }&=g^{\mu \nu }Q_{\alpha \mu \nu } =Q_{\alpha \nu }^{\nu }, \end{aligned}$$
(7a)
$$\begin{aligned} {\tilde{Q}}_{\alpha }&=g^{\mu \nu }Q_{\mu \alpha \nu } =Q_{\alpha \nu }^{\nu } . \end{aligned}$$
(7b)

It is useful to define a superpotential, given by

$$\begin{aligned} P_{\mu \nu }^{\beta }=-\frac{1}{2}L_{ \mu \nu }^{\beta }-\frac{1}{4}\left[ \left( {\tilde{Q}}^{\beta } -Q^{\beta }\right) g_{\mu \nu } +\delta _{(\mu }^{\beta }Q_{\nu )}\right] , \end{aligned}$$
(8)

so that the nonmetricity tensor (6) with the superpotential (8), defines the nonmetricity scalar from the contraction as follows

$$\begin{aligned} {\mathbb {Q}}=-Q_{\beta \mu \nu }P_{}^{\beta \mu \nu }. \end{aligned}$$
(9)

Applying the corresponding geometric configurations where the curvature is zero, i.e. the teleparallel condition, and where the geometry has no torsion, we obtain the following relation

$$\begin{aligned}&\overset{C}{R}{}_{\alpha \mu \nu }^{\beta } +\overset{C}{\nabla }_{\mu }L_{\nu \alpha }^{\beta } -\overset{C}{\nabla }_{\nu }L_{\mu \alpha }^{\beta } +L_{\mu \rho }^{\beta }L_{\nu \alpha }^{\rho } -L_{\nu \rho }^{\beta }L_{\mu \alpha }^{\rho } =0, \end{aligned}$$
(10)

where \(\overset{C}{\nabla }\) denotes the derivative that corresponds to the Christoffel symbol (2), as well as the curvature tensor \(\overset{C}{R}{}_{\alpha \mu \nu }^{\beta }\), which has its explicit form given by

$$\begin{aligned} \overset{C}{R}{}_{\alpha \mu \nu }^{\beta }= & {} \partial _{\alpha }\Gamma _{\nu \mu }^{\beta } -\partial _{\nu }\Gamma _{\alpha \mu }^{\beta } +\Gamma _{\alpha \rho }^{\beta }\Gamma _{ \nu \mu }^{\rho }-\Gamma _{\nu \rho }^{\beta }\Gamma _{ \alpha \mu }^{\rho }. \end{aligned}$$
(11)

The contraction of this tensor produces the Ricci tensor

$$\begin{aligned} \overset{C}{R}_{\mu \nu }=\overset{C}{R}{}_{ \mu \beta \nu }^{\beta }, \end{aligned}$$
(12)

and Eq. (10) reduces to the following expression after the appropriate contractions

$$\begin{aligned} \overset{C}{R}={\mathbb {Q}}-\overset{C}{\nabla }_{\beta } \left( Q_{}^{\beta } -{\tilde{Q}}{}_{}^{\beta }\right) , \end{aligned}$$
(13)

where

$$\begin{aligned} \overset{C}{R}= g^{\mu \nu }\overset{C}{R}_{\mu \nu }, \end{aligned}$$
(14)

is the contraction of Eq. (12), expressed in terms of (2).

Thus, taking into Eq. (13) we verify that the non-metricity scalar differs from the Ricci scalar by a term of the total derivation

$$\begin{aligned} {\mathbb {B}}_Q=\overset{C}{\nabla }_{\beta }\left( Q^{\beta } -{\tilde{Q}}^{\beta }\right) , \end{aligned}$$
(15)

which will be a fundamental quantity used throughout this work.

2.2 Coincident gauge

Let us now discuss an interesting symmetry property used in STEGR that simplifies this formulation somewhat, where the coordinate transformation considered here is known as the “coincident gauge” [22]. More specifically, a variety of flat geometry, where we have a flat connection with \(T^\beta _{\mu \nu }=0\), can be described by

$$\begin{aligned} {\bar{\Gamma }}{}_{\mu \nu }^{\beta } =\left( \Lambda ^{-1}\right) _{\alpha }^{\beta } \partial _{\mu }\Lambda _{\nu }^{\alpha }. \end{aligned}$$
(16)

where \(\Lambda _{\mu }^{\alpha }\) belongs to \(\textrm{GL}(4, \mathcal{{\mathbb {R}}})\).

The constraint in which the torsion is zero implies

$$\begin{aligned} \left( \Lambda ^{-1}\right) _{\beta }^{\alpha }\partial _{[\mu }\Lambda _{\nu ]}^{\beta }=0, \end{aligned}$$
(17)

so that it is possible to rewrite the connection (16) from a parametrization according to the transformation \(\Lambda _{\mu }^{\alpha }=\partial _{\mu }\xi ^{\alpha }\). Consequently, we obtain

$$\begin{aligned} {\bar{\Gamma }}{}_{\mu \nu }^{\beta } =\frac{\partial x^{\beta }}{\partial \xi ^{\rho }}\partial _{\mu }\partial _{\nu }\xi ^{\rho }, \end{aligned}$$
(18)

where \(\xi ^{\lambda }\) is a set of arbitrary functions of the coordinates \(x^{\lambda }\).

As mentioned above, this result makes it possible, by a suitable choice of coordinates, to make the connection flat, torsion-free, and therefore consequently vanishes. The choice of coordinates that cancels the connection is given by \(x^\lambda =\xi ^\lambda \), thus

$$\begin{aligned} {\bar{\Gamma }}{}_{\mu \nu }^{\beta }=0, \end{aligned}$$
(19)

due \(\partial _{\mu }\partial _{\nu }\xi ^{\rho }=0\).

The coordinate transformation that implies a connection described by the above result is called the “coincident gauge” [22, 51], so that we have

$$\begin{aligned} \Gamma _{\mu \nu }^{\beta }=-L_{\mu \nu }^{\beta }. \end{aligned}$$
(20)

The non-metricity tensor defined according to Eq. (6) now becomes

$$\begin{aligned} Q_{\beta \mu \nu }\equiv \partial _{\beta }g_{\mu \nu }. \end{aligned}$$
(21)

The coincident gauge is an extremely useful symmetry, which will be used to develop the solutions outlined throughout this work.

2.3 \(f({\mathbb {Q}})\) gravity

The STEGR describes the gravitational interaction by the non-metricity tensor, \(Q_{\alpha \mu \nu }\), and the action is defined as:

$$\begin{aligned} S_{\text {STEGR}}=\int \sqrt{-g}d^{4}x\left( {\mathbb {Q}} +2\kappa ^2\mathcal {L}_m\right) , \end{aligned}$$
(22)

where \(\kappa ^2=8\pi G/c^4\), G is the gravitational constant and \(\mathcal {L}_{m}\) is the Lagrangian of the matter field.

Note that from the relation (13), the action of STEGR differs from the Einstein–Hilbert action of GR by a boundary term \({\mathbb {B}}_Q\), given by Eq. (15), which means that STEGR is an equivalent formulation of GR.

A nonlinear generalization of STEGR is a proposal in which the action is described as follows

$$\begin{aligned} S_{ f}=\int \sqrt{-g}d^{4}x\left[ f\left( \mathbb {{Q}}\right) +2\kappa ^2\mathcal {L}_{m}\right] , \end{aligned}$$
(23)

where \(f({\mathbb {Q}})\) is an arbitrary function of the nonmetricity scalar \({\mathbb {Q}}\).

However, this theory raises some issues, for instance, if one chooses the coincident gauge with static and spherical symmetry, the equations of motion imply a theory with a linear \(f({\mathbb {Q}})\) function or constant scalar non-metricity, as verified in [52]. This is undoubtedly due to the use of a geometric object which is not a Lorentz invariant. To overcome this challenge, a theory must be constructed that accounts for this invariance.

2.4 \(f({\mathbb {Q}},{\mathbb {B}}_Q)\) gravity

In this work, we develop solutions based on a recent proposal by extending the theory \(f({\mathbb {Q}})\), namely, \(f({\mathbb {Q}},{\mathbb {B}}_Q)\) gravity [24, 53], where \({\mathbb {B}}_Q\) is the boundary term given by Eq. (15). In this approach, the action is defined as:

$$\begin{aligned} S_{ f({\mathbb {Q}},{\mathbb {B}}_Q)}= & {} \int \sqrt{-g}d^{4}x \Big [f\left( \mathbb {{Q}},{\mathbb {B}}_Q\right) +2\kappa ^2\mathcal {L}_{m}\Big ]. \end{aligned}$$
(24)

The variation of the action (24) with respect to the metric tensor, yields the following equation of motion:

$$\begin{aligned}{} & {} f_{Q}\left( {\mathbb {Q}}\right) \overset{C}{G}_{\mu \nu } -\frac{1}{2}g_{\mu \nu }\bigr [f({\mathbb {Q}}) -f_{Q}\left( {\mathbb {Q}}\right) {\mathbb {Q}}\nonumber \\{} & {} \quad -f_{B}\left( {\mathbb {B}}_Q\right) {\mathbb {B}}_Q\bigl ] +2P_{\mu \nu }^{\alpha }\partial _{\alpha } \bigr [f_{Q}\left( {\mathbb {Q}}\right) \nonumber \\{} & {} \quad +f_{B}\left( {\mathbb {B}}_Q\right) \bigl ] -g_{\mu \nu }\overset{C}{\square }\ f_{B}\left( {\mathbb {B}}_Q\right) \nonumber \\{} & {} \quad +\overset{C}{\nabla }_{\mu }\overset{C}{\nabla }_{\nu }f_{B} \left( {\mathbb {B}}_Q\right) =\kappa ^{2}\Theta _{\mu \nu }, \end{aligned}$$
(25)

where \(\overset{C}{G}_{\mu \nu }\) is the Einstein tensor of GR and we denote \(f_{Q}\left( {\mathbb {Q}},{\mathbb {B}}_Q\right) =\partial f\left( {\mathbb {Q}},{\mathbb {B}}_Q\right) /\partial {\mathbb {Q}}\), \(f_{B}\left( {\mathbb {Q}},{\mathbb {B}}_Q\right) = \partial f\left( {\mathbb {Q}},{\mathbb {B}}_Q\right) /\partial {\mathbb {B}}_Q\), the energy–momentum tensor is \(\Theta _{\mu \nu }\) and \(\square =\nabla ^{\mu }\nabla _{\mu }\). Note that the contribution of the boundary term leads to the existence of fourth order field equations in this theory [24].

Now, the variation of the action (24) with respect to the affine component provides the following field equation

$$\begin{aligned} \nabla _{\mu }\nabla _{\nu }\Bigl [\sqrt{-g} P_{\alpha }^{\mu \nu } \Bigl (f_{Q}\left( {\mathbb {Q}}\right) +f_{B}\left( {\mathbb {B}}_Q\right) \Bigr )\Bigr ]=0. \end{aligned}$$
(26)

If we substitute \({\mathbb {B}}_Q=0\) into the field equation (25), leading to \(f_{B}({\mathbb {Q}},{\mathbb {B}}_Q)\)=0, we arrive at the field equations of \(f({\mathbb {Q}})\) gravity. Note that the resulting equations of motion can be rewritten as [54,55,56]

$$\begin{aligned}{} & {} f_{{Q}}\left( {\mathbb {Q}}\right) \overset{C}{G}_{\mu \nu } -\frac{1}{2}g_{\mu \nu }\Bigl (f({\mathbb {Q}}) -f_{Q}\left( {\mathbb {Q}}\right) {\mathbb {Q}}\Bigr ) \nonumber \\{} & {} \quad +2P_{\mu \nu }^{\alpha }\partial _{\alpha } f_{Q}\left( {\mathbb {Q}}\right) =\kappa ^2\Theta _{\mu \nu }. \end{aligned}$$
(27)

An interesting property is that f(R) gravity is recovered by considering \(R={\mathbb {Q}}-{\mathbb {B}}_Q\), i.e., \(f(R)=f({\mathbb {Q}},{\mathbb {B}}_Q)=f({\mathbb {Q}}-{\mathbb {B}}_Q)\), which yields the following field equations

$$\begin{aligned}{} & {} \left( \overset{C}{G}_{\mu \nu }+\frac{1}{2}g_{\mu \nu }R +g_{\mu \nu }\overset{C}{\square } -\overset{C}{\nabla }_{\mu }\overset{C}{\nabla }_{\nu }\right) f_{R} \left( R\right) \nonumber \\{} & {} \quad -\frac{1}{2}g_{\mu \nu }f\left( R\right) =\kappa ^{2}\Theta _{\mu \nu }. \end{aligned}$$
(28)

3 \(f({\mathbb {Q}},{\mathbb {B}}_Q)\) gravity coupled to non-linear electrodynamics

3.1 Action and metric

In this section, we consider \(f({\mathbb {Q}},{\mathbb {B}}_Q)\) gravity coupled to NLED. The action is given by

$$\begin{aligned} S=\int d^{4}x\sqrt{-g}\Bigl [f({\mathbb {Q}},{\mathbb {B}}_{Q}) + 2\kappa ^2 {{\mathcal {L}}}_{\textrm{NLED}}({F})\Bigr ], \end{aligned}$$
(29)

where \({{\mathcal {L}}} _{\textrm{NLED}}(F)\) is the Lagrangian density describing the NLED that depends on the electromagnetic scalar F as

$$\begin{aligned} F=\frac{1}{4}F^{\mu \nu }F_{\mu \nu }, \end{aligned}$$
(30)

and the electromagnetic field \(F_{\mu \nu } \), defined as

$$\begin{aligned} F_{\mu \nu } = \partial _\mu A_\nu -\partial _\nu A_\mu , \end{aligned}$$
(31)

is the Maxwell–Faraday antisymmetric tensor, where \({A_\beta }\) is the magnetic vector potential.

Varying the action (29) with respect to \(A_\beta \), we have

$$\begin{aligned} \nabla _\mu [{{\mathcal {L}}}_F F^{\mu \nu }]=\partial _\mu [\sqrt{-g} {{\mathcal {L}}}_F F^{\mu \nu }]=0, \end{aligned}$$
(32)

where we denote \({{\mathcal {L}}}_F = \partial {{\mathcal {L}}} _{\textrm{NLED}}(F)/\partial F\).

Finally, varying the action (29) with respect to the metric, we obtain the equations of motion (25) and the matter contributions expressed by the energy–momentum tensor, with the following contribution

$$\begin{aligned} \Theta _{\nu }^\mu = \overset{F}{\Theta } {_{\nu }^\mu } , \end{aligned}$$
(33)

where the nonlinear electromagnetic energy–momentum tensor \(\overset{F}{\Theta }{}_{\nu }^{\mu }\) is defined as

$$\begin{aligned} \overset{F}{\Theta }{}_{\nu }^{\mu }= \delta ^\mu _\nu {{\mathcal {L}}}_{\textrm{NLED}} (F) +{{\mathcal {L}}}_F F^{\mu \alpha } F_{\nu \alpha }. \end{aligned}$$
(34)

The solutions that we discuss later in the context of black holes and regular black holes will be obtained using the following static and spherically symmetric metric

$$\begin{aligned} ds^2=A(r)dt^2-\frac{1}{A(r)}dr^2-\Sigma ^2(r)\left( d\theta ^{2} +\sin ^{2}\theta \, d\phi ^{2}\right) , \end{aligned}$$
(35)

where A(r) and \(\Sigma (r)\) are functions of the radial coordinate r.

Here, we only consider the magnetic charge, so that the only non-zero component of the tensor \(F_{\mu \nu }\) is

$$\begin{aligned} F_{23}=q \sin \theta , \end{aligned}$$
(36)

and the electromagnetic scalar is given by

$$\begin{aligned} F=\frac{q^2}{2 \Sigma ^4(r)}. \end{aligned}$$
(37)

In this way, using the metric (35), we can describe boundary term (15) as follows

$$\begin{aligned} {\mathbb {B}}_Q(r)= & {} A''(r)+ \frac{6 A'(r) \Sigma '(r)+4 A(r) \Sigma ''(r)}{\Sigma (r)} \nonumber \\{} & {} +\frac{4A(r)\Sigma '^{2} (r)-2}{\Sigma (r)}, \end{aligned}$$
(38)

and the nonmetricity scalar (9) now becomes

$$\begin{aligned} {\mathbb {Q}} (r)=-\frac{2 \Sigma '(r) \Bigl (\Sigma (r) A'(r)+A(r) \Sigma '(r)\Bigr )}{\Sigma ^2 (r)}. \end{aligned}$$
(39)

3.2 Equations of motion

Thus, taking into account the field equation (25), the line element (35), and the energy–momentum tensor (33), the equations of motion for \(f({\mathbb {Q}},{\mathbb {B}}_{Q})\) gravity coupled with NLED are given by

$$\begin{aligned}{} & {} \frac{1}{2}f_{Q}(r)\left[ {\mathbb {Q}}(r)-\frac{2 \Bigl (\Sigma (r)A'(r)\Sigma '(r) +A(r)\Sigma '(r)^{2}-1\Bigr )}{\Sigma (r)^{2}}\right] + \frac{1}{2}{\mathbb {B}}_Q(r)f_{B}(r) -\frac{f(r)}{2} -\frac{f_{B}^{\prime }(r)\Big (\Sigma (r)A'(r) +4A(r)\Sigma '(r)\Big )}{2\Sigma (r)} \nonumber \\{} & {} \quad = \kappa ^{2}{{\mathcal {L}}}_{\textrm{NLED}}(r), \end{aligned}$$
(40)
$$\begin{aligned}{} & {} \frac{1}{2} A(r) \cot \theta \left( f_B'(r)+f_Q'(r)\right) =0, \end{aligned}$$
(41)
$$\begin{aligned}{} & {} \frac{\cot \theta \left( f_B'(r)+f_Q'(r)\right) }{2 \Sigma ^2 (r)}=0, \end{aligned}$$
(42)
$$\begin{aligned}{} & {} \frac{1}{2}f_{Q}(r)\left\{ {\mathbb {Q}}(r) -\frac{2\left[ \Sigma (r)\Bigl (A^{\prime }(r)\Sigma ^{\prime }(r)+2A(r) \Sigma ^{\prime \prime }(r)\Bigr ) +A(r)\Sigma ^{\prime 2}(r)-1\right] }{\Sigma (r)^{2}}\right\} -\frac{2A(r)\Sigma ^{\prime }(r)\left( f_{B}^{\prime }(r) +f_{Q}^{\prime }(r)\right) }{\Sigma (r)} \nonumber \\{} & {} \quad -\frac{f(r)}{2} +\frac{1}{2}\left[ -\frac{f_{B}^{\prime }(r)\Bigl (\Sigma (r)A'(r) +4A(r)\Sigma ^{\prime }(r)\Bigr )}{\Sigma (r)} -2A(r)f_{B}^{\prime \prime }(r)+{\mathbb {B}}_Q(r)f_{B}(r)\right] \nonumber \\{} & {} \quad = \kappa ^{2}{{\mathcal {L}}}_{\textrm{NLED}}(r), -\frac{f_{B}^{\prime }(r)\left( \Sigma (r)A^{\prime }(r) +A(r)\Sigma ^{\prime }(r)\right) }{\Sigma (r)}-A(r)f_{B}^{\prime \prime }(r)\end{aligned}$$
(43)
$$\begin{aligned}{} & {} \quad +\frac{1}{2}f_{Q}(r)\left( {\mathbb {Q}}(r) -\frac{\Sigma (r)A^{\prime \prime }(r) +2A^{\prime }(r)\Sigma ^{\prime }(r)+2A(r)\Sigma ^{\prime \prime } (r)}{\Sigma (r)}\right) +\frac{1}{2}{\mathbb {B}}_Q(r)f_{B}(r) \nonumber \\{} & {} \quad -\frac{1}{2}\left[ f(r)+\frac{\Big (\Sigma (r)A'(r) +2A(r)\Sigma '(r)\Big )\Big (f_{B}^{\prime }(r) +f_{Q}^{\prime }(r)\Big )}{\Sigma (r)}\right] \nonumber \\{} & {} \quad =\kappa ^{2}\left( {{\mathcal {L}}}_{\textrm{NLED}}(r)-\frac{q^{2}L_{F}(r)}{\Sigma ^{4}(r)}\right) , \end{aligned}$$
(44)

where the prime denotes a derivative with respect to the radial coordinate r.

4 Black hole solutions in \(f(\mathbb {Q,B_Q})\) gravity

In this section, in order to deduce black hole solutions, consider the following metric

$$\begin{aligned} ds^2= & {} A(r)dt^2-\frac{1}{A(r)}dr^2\nonumber \\{} & {} -r^2 \left( d\theta ^{2} +\sin ^{2}\theta \,d\phi ^{2}\right) , \end{aligned}$$
(45)

where we have used \( \Sigma (r) = r\) for simplicity.

Thus the boundary term (38) is now given as

$$\begin{aligned} {\mathbb {B}} _Q(r)=-\frac{r^2 A''(r)+6 r A'(r)+4 A(r)-2}{r^2}, \end{aligned}$$
(46)

and the non-metricity scalar (39) takes the form

$$\begin{aligned} {\mathbb {Q}}(r)=-\frac{2 \Bigl (r A'(r)+A(r) \Bigr )}{r^2}, \end{aligned}$$
(47)

where Eqs. (46) and (47) satisfy condition (13).

To obtain black hole solutions, and in order to ensure that our solutions preserve Lorentz symmetry, we assume the following relation, throughout this work:

$$\begin{aligned} f_Q = -f_B. \end{aligned}$$
(48)

Equation (48) contains the particular case of the symmetric teleparallel theory \(f({\mathbb {Q}},{\mathbb {B}}_Q)={\mathbb {Q}}-{\mathbb {B}}_Q\), and even if \(f_Q=1\) and \(f_B=-1\) this equation is identically satisfied. In the coincident gauge, condition (48) must always be enforced, where the equations of motion (41) and (42) are identically satisfied. If Eq. (48) is not imposed for the entire radial coordinate range, two of the equations of motion, i.e., (41) and (42), will not be satisfied for some value or interval of r. We conclude that condition (48) is a necessary imposition for the coincident gauge in the theory \(f({\mathbb {Q}},{\mathbb {B}}_Q)\) with static spherically symmetry. It is likely that this condition can be relaxed for metrics that also depend on time, such as the Lemaître–Tolman–Bond metric, as is the case with the f(T) theory [57].

Thus, solving the equations of motion (43)–(44), we obtain the following quantities:

$$\begin{aligned}{} & {} {{\mathcal {L}}} _{\textrm{NLED}} (r)= \frac{1}{2 \kappa ^2}\left[ \frac{r^2 f'(r) \Bigl (r A''(r)+2 A'(r)\Bigr )}{r^3 A^{3}(r)+6 r^2 A''(r)-2 r A'(r)-8 A(r)+4} -f(r) \right] , \end{aligned}$$
(49)
$$\begin{aligned}{} & {} {{\mathcal {L}}} _F (r)= \frac{r^5 \Bigl (r^2 A''(r)-2 A(r)+2\Bigr ) f'(r)}{2 \kappa ^2 q^2 \Bigl (r^3 A^{3}(r)+6 r^2 A''(r)-2 r A'(r)-8 A(r)+4\Bigr )}. \end{aligned}$$
(50)

4.1 First model of black hole

It is an extremely difficult task to solve the equations of motion, and therefore we propose several simplifications. As a starting point, we model the boundary term described by Eq. (46) and assume that it is a regular quantity, so that we can solve the differential equation to find an expression for the metric function A(r).

Thus, consider the following boundary term

$$\begin{aligned} {\mathbb {B}}_Q(r)=\frac{b_0}{r^2}, \end{aligned}$$
(51)

so that the black hole solution can therefore be determined from Eq. (46) as follows:

$$\begin{aligned} -\frac{r^2 A''(r)+6 r A'(r)+4 A(r)-2}{r^2}=\frac{b_0}{r^2}. \end{aligned}$$
(52)

which provides the following metric function

$$\begin{aligned} A(r)=\frac{2-b_0}{4}+\frac{b_1}{r^4}-\frac{2 M}{r}, \end{aligned}$$
(53)

where \(b_0\) an \(b_1\) are constants. Note that we obtain the Schwarzschild metric function by imposing \(b_0=-2\) and \(b_1=0\) into Eq. (53). In the limit \(r \rightarrow 0\), the metric function diverges, while in the limit \(r \rightarrow \infty \), it is asymptotically flat. If we set \(b_0=-2\), the metric function (53) is asymptotically Minkowskian.

The non-metricity scalar for this solution is given by

$$\begin{aligned} {\mathbb {Q}}(r)=\frac{b_0-2}{2 r^2}+\frac{6 b_1}{r^6}. \end{aligned}$$
(54)

Note that the non-metricity scalar and the boundary term satisfy the condition (13). Although \({\mathbb {Q}}(r)\) and \({\mathbb {B}}_Q(r)\) are divergent in the limit \(r\rightarrow 0\), the combination \({\mathbb {Q}}(r)-{\mathbb {B}}_Q(r)=R\) is what characterizes whether the solution is regular or not. Thus, as \({\mathbb {Q}}(r)-{\mathbb {B}}_Q(r)\) is singular at \(r=0\), the geometry of the spacetime is that of a singular solution (we refer the reader to the comment in Ref. [52] under Eq. (46)).

To determine the existence of event horizons in a given solution, it is necessary to solve the following equation:

$$\begin{aligned} A(r_H)=0, \end{aligned}$$
(55)

where \(r_{H}\) is the radius of the horizon.

In addition, we can use the following condition to obtain the values of the critical parameters

$$\begin{aligned} \frac{d A{(r)}}{dr}\bigg |_{r=r_H}=0. \end{aligned}$$
(56)

We can therefore solve Eqs.  (55) and (56) simultaneously to obtain algebraic solutions for the radius of the event horizon and the value of a certain critical parameter in the model. However, depending on the model, it may be difficult to find analytical solutions for these quantities. Therefore, we will approach this problem numerically by assigning values to the constants to investigate the presence of event horizons through graphical representations.

Therefore, to find numerical solutions, we set the following values for the parameters: \(b_0=0.5\) and \(M=1.0\). If the two conditions (55) and (56) are satisfied simultaneously, we obtain the critical parameter value \(b_{1c}=32\). Figure 1 shows the behavior of the metric function (53), where \(b_{1c}\) takes the values \(b_1>b_{1c}\), \(b_1=b_{1c}\) and \(b_1<b_{1c}\). For \(b_1>b_{1c}\) there are no event horizons; for \(b_1=b_{1c}\) there is only one horizon; and for \(b_1<b_{1c}\) two event horizons are present. We can see in Fig. 1 that all the curves for the different values of \(b_1\) with \(r\gg 1\) match.

The quantities \({{\mathcal {L}}}_{\textrm{NLED}}\) and \({{\mathcal {L}}}_F\) take the form

$$\begin{aligned}{} & {} {{\mathcal {L}}}_{\textrm{NLED}}(r)=-\frac{1}{2 \kappa ^2} \Biggl ( f(r) -\frac{6 b_1 f'(r)}{b_0 r^3}\Biggr ), \end{aligned}$$
(57)
$$\begin{aligned}{} & {} {{\mathcal {L}}}_F(r)=\frac{r \Bigl [\bigl (b_0+2\bigr ) r^4 +36 b_1\Bigr ] f'(r)}{8 b_0 \kappa ^2 q^2}, \end{aligned}$$
(58)

respectively, and indeed, these obey the consistency relation

$$\begin{aligned} {{\mathcal {L}}}_F-\frac{\partial {{\mathcal {L}}} _{\textrm{NLED}}}{\partial r} \bigg (\frac{\partial F}{\partial r}\bigg )^{-1}=0. \end{aligned}$$
(59)

If we solve the consistency condition (59), we get the following solution

$$\begin{aligned} f(r)=-\frac{f_0 r \root 4 \of {3}\; \Gamma \!\left( \frac{1}{4},-\frac{(b_0-2) r^4}{48 b_1}\right) }{2 \root 4 \of {-\frac{(b_0-2) r^4}{b_1}}}+f_1, \end{aligned}$$
(60)

where \(f_0\) and \(f_1\) are the constants of the integration.

Consequently, we can express the function as follows

$$\begin{aligned}{} & {} f({\mathbb {Q}},{\mathbb {B}}_Q) =-\frac{\sqrt{\frac{\xi }{{\mathbb {Q}} -{\mathbb {B}}_{Q}}}}{4\sqrt{6}}\text {E}_{\frac{3}{4}}\left( -\frac{(b_{0}-2)\xi ^{2}}{1728\,b_{1} \left( {\mathbb {Q}}-{\mathbb {B}}_{Q}\right) ^{2}}\right) .\nonumber \\ \end{aligned}$$
(61)

where \({{\text {E}}}_k\) is the exponential integral function

$$\begin{aligned} {\text {E}}_{k}\left( w\right) =\int _{1}^{\infty }\frac{e^{-wt}}{t^{k}}dt, \end{aligned}$$
(62)

and \(\xi \) is given by

$$\begin{aligned} \xi= & {} \root 3 \of {36\sqrt{b_{0}({\mathbb {Q}} -{\mathbb {B}}_Q)^{2} \left( 324b_{1}({\mathbb {Q}} -{\mathbb {B}}_Q)^{2}-(b_{0}+2)^{3}\right) } -(b_{0}+2)^{3}+648b_{1}({\mathbb {Q}}-{\mathbb {B}}_Q)^{2}} \nonumber \\{} & {} +\frac{(b_{0}+2)^{2}}{\root 3 \of {36\sqrt{b_1({\mathbb {Q}} -{\mathbb {B}}_Q)^{2}\left( 324 b_1({\mathbb {Q}}-{\mathbb {B}}_Q)^{2}-(b_{0}+2)^{3}\right) } -(b_{0}+2)^{3}+648b_{1}({\mathbb {Q}}-{\mathbb {B}}_Q)^{2}}}-b_{0}-2. \end{aligned}$$
(63)

By choosing \(b_0=-2\) and \(b_1=1\) in Eq. (61), we get the form of

$$\begin{aligned} f({\mathbb {Q}},{\mathbb {B}}_Q) = -\frac{\root 6 \of {3}\; {\text {E}}_{\frac{3}{4}} \left( \frac{1}{2 \root 3 \of {6} ({\mathbb {Q}}-{\mathbb {B}}_Q)^{2/3}}\right) }{2\ 2^{5/6} \root 6 \of {{\mathbb {Q}}-{\mathbb {B}}_Q}} \end{aligned}$$
(64)

that permits us see more clearly that the Eq. (64) satisfies \(f_Q=-f_B\).

We have investigated the asymptotic behavior of function (61) for small and very large values of \({\mathbb {Q}}-{\mathbb {B}}_Q\). We found that there is no linear behavior in these two cases. To illustrate this, in Fig. 2 we compared function (61) for small values of \({\mathbb {Q}}-{\mathbb {B}}_Q\) using the blue curve with the linear case represented by a dashed red line. However, if we analyze very small intervals of \({\mathbb {Q}}-{\mathbb {B}}_Q\), the behavior we observe is practically linear.

Fig. 1
figure 1

The plot depicts the metric function (53), i.e., A(r), for the specific values of \(b_0=0.5\) and \(M =1.0\), where Eq. (56) determines the critical value of the constant \(b_1\), given by \(b_{1c} = 32.0\). For \(b_1 > b_{1c}\) there is no horizon, for \(b_{1} = b_{1c}\), there is only one horizon and for \(b_{1} < b_{1c}\), there are two horizons

Full size image

From the function f(r) at our disposal, the quantities are now expressed as:

$$\begin{aligned}{} & {} f_B(r)=-\frac{r^3 e^{\frac{(b_0-2) r^4}{48 b_1}}}{2 b_0}, \end{aligned}$$
(65)
$$\begin{aligned}{} & {} {{\mathcal {L}}}_{\textrm{NLED}}(r)=\frac{r \, {\text {E}}_{\frac{3}{4}}\left( -\frac{(b_0-2) r^4}{48 b_1}\right) }{8 \kappa ^2}+\frac{3 b_1 e^{\frac{(b_0-2) r^4}{48 b_1}}}{b_0 \kappa ^2 r^3}, \end{aligned}$$
(66)
$$\begin{aligned}{} & {} {{\mathcal {L}}}_F(r)=\frac{r e^{\frac{(b_0-2) r^4}{48 b_1}} \left( (b_0+2) r^4+36 b_1\right) }{8 b_0 \kappa ^2 q^2}, \end{aligned}$$
(67)

where Eqs. (66) and (67) satisfy the consistency condition (59).

Finally, using Eq. (30), we derive the expression for r(F), and consequently, we are able to express the electromagnetic Lagrangian in terms of F, which is given by

$$\begin{aligned}{} & {} {{\mathcal {L}}} _{\textrm{NLED}} (F) \nonumber \\{} & {} \quad =\frac{b_0 q^2 \, {\text {E}}_{\frac{3}{4}}\left( -\frac{(b_0-2) q^2}{96 b_1 F}\right) +48 b_1 F e^{\frac{(b_0-2) q^2}{96 b_1 F}}}{8 \root 4 \of {2} b_0 \root 4 \of {F} \kappa ^2 q^{3/2}}. \end{aligned}$$
(68)

We illustrate the behavior of the Lagrangian (68) in Fig. 3 for \(q=0.7\), \(b_0=0.5\) and three different values of \(b_1\). A note that setting \(b = 2\) is an inconsistent value, as it lies outside the domain of our solution.

Let us now analyze the asymptotic limits of the Lagrangian described by Eq. (68). For instance, for \(F\gg 1\), we observe that

$$\begin{aligned} {{\mathcal {L}}} _{\textrm{NLED}} (F)\sim \frac{3\ 2^{3/4} b_1 F^{3/4}}{b_0\kappa ^2 q^{3/2}}. \end{aligned}$$
(69)

Whereas for \(F\ll 1\),

$$\begin{aligned}{} & {} {{\mathcal {L}}} _{\textrm{NLED}} (F)\sim \exp \left[ \frac{\left( b_{0}-2\right) q^{2}}{96b_{1}F}\right] \nonumber \\{} & {} \quad \left[ -\frac{3\left( 2^{3/4}\left( b_{0} +2\right) b_{1}\right) F^{3/4}}{\left( b_{0}-2\right) b_{0}q^{3/2}\kappa ^{2}}\right] . \end{aligned}$$
(70)

Therefore, based on results (69) and (70), we find that we have not recovered the linear case.

Fig. 2
figure 2

Graphic representation of function \(f({\mathbb {Q}},{\mathbb {B}}_Q)\), described by expression (61). Where we consider \(M=10\), \(b_0=-2\) and \(b_1=1\). The blue curve represents the behavior of our model for \(({\mathbb {Q}}-{\mathbb {B}}_Q)\ll 1\), while the dashed red line represents the linear case

Full size image
Fig. 3
figure 3

Graphic representation of \({{\mathcal {L}}}_{\textrm{NLED}}(F)\), described by expression (68), with respect to F. Where we consider \(q=0.7\) and \(b_0=0.5\)

Full size image

If we analyze the graph in Fig. 3, we find that the Lagrangian \({{\mathcal {L}}}_{\textrm{NLED}} (F)\) described by Eq. (68), is nonlinear. Another important feature is that it is not “function multivalued”. We also note that this Lagrangian is a bijective function, since each element in the domain maps to only one element in the codomain. This is a peculiarity of magnetically charged solutions where we have a unique representation, which means that a single Lagrangian is sufficient to describe a single solution. This aspect does not apply to electrically charged solutions, where the Lagrangians are no longer bijective, as for one value of F there can be one or more values of \({{\mathcal {L}}}_{\textrm{NLED}}(F)\) (see Refs. [45, 58]). Furthermore, we see in the graph that the Lagrangian \({{\mathcal {L}}}_{\textrm{NLED}}(F)\) has no peaks, which is another characteristic feature of magnetically charged solutions. We can also observe that by increasing the value of \(b_1\) while maintaining q and \(b_0\) fixed, the Lagrangian \({{\mathcal {L}}}_{\textrm{NLED}}(F)\) assumes increasing values.

4.2 Second model of black hole

In this model we adopt the symmetry presented in Eq. (48), as before, since it gives us analytical solutions. We model the non-metricity scalar from Eq. (47) so that it becomes a regular quantity, which allows us to obtain a differential equation to find an expression for A(r). For instance, consider the following non-metricity scalar

$$\begin{aligned} {\mathbb {Q}}(r)=\frac{q_0}{r^2+r_0^2}, \end{aligned}$$
(71)

so that the solution is therefore characterized by Eq. (47), and is given by

$$\begin{aligned} -\frac{2 \Bigl (r A'(r)+A(r)\Bigr )}{r^2}=\frac{q_0}{r^2+r_0^2}. \end{aligned}$$
(72)

which yields the following metric function

$$\begin{aligned} A(r)=-\frac{2 M}{r}-\frac{q_0}{2} \left[ 1-\frac{r_0}{r} \tan ^{-1}\left( \frac{r}{r_0}\right) \right] . \end{aligned}$$
(73)

where \(q_0\) an \(r_0\) are constants. The metric function (73) diverges when we take the limit \(r\rightarrow 0\), however, for the limit \(r\rightarrow \infty \), this A(r) function is asymptotically flat. If we set \(q_0=-2\) and \(r_0=0\) in Eq. (73), we recover the Schwarzschild solution.

We may consider once again Eqs. (55) and (56) to determine the radius of the event horizon and the critical parameter by numerical solutions. However, for the specific model described by the metric function (73), we did not find an extreme value for \(q_{0}\). In Fig. 4, we depict the behavior of A(r). When \(q_0=0\), we observe a curve that is consistently negative and tends towards zero. Conversely, for \(q_0> 0\) there is no horizon, while for \(q_0 < 0\) there is an horizon. Equation (73) shows that the function A(r) tends to \(-q_0/2\) in the limit when r tends to infinity. This behavior can be illustrated in Fig. 4, where for \(q_0 > 0\) the metric function A(r) stabilizes at a fixed negative value (as shown by the green curve), while for \(q_0 < 0\) A(r) tends to a fixed positive value (as shown by the red curve).

Fig. 4
figure 4

The plot depicts the metric function (73), i.e., A(r), for the specific values of \(b_0=0.1\) and \(M =1.0\). For \(q_0 > 0\) there is no horizon, for \(q_{0} = 0\), the curve is negative and tends to zero and for \(q_{0} < 0\), there is only one horizon

Full size image

In this model the boundary term is given by

$$\begin{aligned} {\mathbb {B}}_Q(r)=\frac{2 (q_0+1) r^4+(3 q_0+4) r^2 r_0^2 +2 r_0^4}{r^2 \left( r^2+r_0^2\right) ^2}. \end{aligned}$$
(74)

Equations (71) and (74) satisfy relation (13). Although \({\mathbb {Q}}(r)\) is regular for for all values of r and \(q_0\), the combination \({\mathbb {Q}}(r)-{\mathbb {B}}_Q(r)=R\) in this model exhibits a singularity in the limit of \(r\rightarrow 0\).

We also find that the quantities \({{\mathcal {L}}}_{\textrm{NLED}}\) and \({{\mathcal {L}}}_F\) take the following forms

$$\begin{aligned}{} & {} {{\mathcal {L}}}_{\textrm{NLED}}(r)=-\frac{f(r)}{2\kappa ^{2}} \nonumber \\{} & {} \quad -\frac{1}{8\kappa ^{2}}\left[ \frac{q_{0}r^{3}r_{0}^{2} \left( r^{2}+r_{0}^{2}\right) f'(r)}{(q_{0}+1)r^{6} +(2q_{0}+3)r^{4}r_{0}^{2}+3r^{2}r_{0}^{4}+r_{0}^{6}}\right] , \nonumber \\ \end{aligned}$$
(75)

and

$$\begin{aligned} {{\mathcal {L}}}_F(r)=\frac{r^{5}\left( r^{2}+r_{0}^{2} \right) f'(r)\Bigl [(q_{0}+2)r^{4} +4r^{2}r_{0}^{2}+2r_{0}^{4}\Bigr ]}{8\kappa ^{2}q^{2}\left[ (q_{0}+1)r^{6}+(2q_{0}+3)r^{4}r_{0}^{2} +3r^{2}r_{0}^{4}+r_{0}^{6}\right] }, \nonumber \\ \end{aligned}$$
(76)

respectively.

Solving the consistency condition (59), we deduce

$$\begin{aligned}{} & {} f(r)=\frac{1}{4}f_{0}\left\{ -\frac{eq_{0}}{r_{0}^{2}}\, \text {Ei}\left( -\frac{r^{2}}{r_{0}^{2}}-1\right) \right. \nonumber \\{} & {} \quad -\frac{2}{r_{0}^{2}} \,\text {Ei}\left( -\frac{r^{2}}{r_{0}^{2}}\right) -\frac{e^{-\frac{r^{2}}{r_{0}^{2}}} \big [(q_{0}+2)r^{4}+2(q_{0}+2)r^{2}r_{0}^{2}+2r_{0}^{4}\big ]}{r^{2}\left( r^{2}+r_{0}^{2}\right) ^{2}}\Biggr \} \nonumber \\{} & {} +f_{1}, \end{aligned}$$
(77)

where \(f_0\) and \(f_1\) are the integration constants, and \(\textrm{Ei}(\omega \)) is the exponential integral function

$$\begin{aligned} \text {Ei}\left( w\right) =\int _{-w}^{\infty }\frac{e^{-t}}{t}dt. \end{aligned}$$
(78)

With the function f(r) at our disposal, the quantities are now expressed in the following form:

$$\begin{aligned} f_Q(r)= & {} \frac{1}{4} e^{-\frac{r^2}{r_0^2}}, \end{aligned}$$
(79)
$$\begin{aligned} {{\mathcal {L}}}_{\textrm{NLED}}(r)= & {} \frac{1}{8\kappa ^{2}}\left\{ \frac{eq_{0}}{r_{0}^{2}}\, \text {Ei}\left( -\frac{r^{2}}{r_{0}^{2}}-1\right) \right. \nonumber \\{} & {} +\frac{2}{r_{0}^{2}}\,\text {Ei}\left( -\frac{r^{2}}{r_{0}^{2}}\right) \nonumber \\{} & {} +\frac{e^{-\frac{r^{2}}{r_{0}^{2}}}\big [(q_{0}+2)r^{2}+2r_{0}^{2}\big ]}{r^{2}\left( r^{2}+r_{0}^{2}\right) }\Biggr \}, \end{aligned}$$
(80)
$$\begin{aligned} {{\mathcal {L}}}_F(r)= & {} \frac{e^{-\frac{r^{2}}{r_{0}^{2}}}\left[ (q_{0}+2)r^{6}+4r^{4}r_{0}^{2}+2r^{2}r_{0}^{4}\right] }{8\kappa ^{2}q^{2}\left( r^{2}+r_{0}^{2}\right) ^{2}}, \end{aligned}$$
(81)

where Eqs. (80) and (81) satisfy the condition (59).

We can also represent the Lagrangian density as follows

$$\begin{aligned} {{\mathcal {L}}} _{\textrm{NLED}}(F)= & {} \frac{1}{8\kappa ^{2}}\bigg [ \frac{e^{-\frac{q}{\sqrt{2F}r_{0}^{2}}} \left( 2\sqrt{F}q(q_{0}+2)+4\sqrt{2}Fr_{0}^{2}\right) }{q\left( 2\sqrt{F}r_{0}^{2}+\sqrt{2}q\right) } \nonumber \\{} & {} +\frac{eq_{0}}{r_{0}^{2}} \text {Ei}\left( -\frac{q}{\sqrt{2F}r_{0}^{2}}-1\right) \nonumber \\{} & {} + \frac{2}{r_{0}^{2}}\text {Ei}\left( -\frac{q}{\sqrt{2F}r_{0}^{2}}\right) \bigg ]. \end{aligned}$$
(82)

The asymptotic limits of the Lagrangian (82) for \(F\gg 1\) are

$$\begin{aligned} {{\mathcal {L}}} _{\textrm{NLED}} (F)&\sim -0.219\frac{eq_{0}}{8\kappa ^{2}r_{0}^{2}}+\frac{\sqrt{F}}{2\sqrt{2}\kappa ^{2}q} \nonumber \\&\quad +\frac{-\ln (F)+q_{0}+2\gamma -2}{8\kappa ^{2}r_{0}^{2}} -\frac{\sqrt{\frac{1}{F}}q\left( 2q_{0}+1\right) }{8\left( \sqrt{2}\kappa ^{2}r_{0}^{4}\right) } \nonumber \\&\quad +\frac{q^{2}(15q_{0}+1)}{96F\kappa ^{2}r_{0}^{6}}, \end{aligned}$$
(83)

where \(\gamma \) is the Euler constant. The numerical value 0.219 came from an approximation of the integral exponential function. We have constructed a plot to demonstrate the asymptotic limit of Lagrangian (82) for \(F\ll 1\). To illustrate this behavior, we compared the curve of our model with a Lagrangian proportional to F, i.e. the linear case, represented by the blue curve and the dashed red line, respectively, as depicted in Fig. 5.

Fig. 5
figure 5

Graphic representation of \({{\mathcal {L}}}_{\textrm{NLED}}(F)\), described by expression (82), with respect to F. Where we consider \(q=0.7\) and \(r_0=0.5\). The blue curve represents the behavior of our model for \(F\ll 1\), while the dashed red line represents the linear case

Full size image

We see that we cannot recover the linear case if we take the asymptotic limits.

The behavior of Lagrangian (82) for three different scenarios of \(q_0\) (\(q_0 > 0\), \(q_0 = 0\), \(q_0 < 0\)) is described in Fig. 6.

Fig. 6
figure 6

Graphic representation of \({{\mathcal {L}}}_{\textrm{NLED}}(F)\), described by expression (82), with respect to F. Where we consider \(q=0.7\) and \(r_0=0.5\)

Full size image

In Fig. 6, we see that we get different curves of \({{\mathcal {L}}}_{\textrm{NLED}}(F)\) (Eq. (82)) for different values of \(q_0\) if we keep the values of q and \(r_0\) constant. Similarly, we could hold the values of \(q_0\) and q constant and vary only the value of \(r_0\) or alternatively hold \(r_0\) and \(q_0\) constant and vary q. From the graph we can also see that the Lagrangian \({{\mathcal {L}}}_{\textrm{NLED}}(F)\) (Eq. (82)) is a monotonically increasing function. This could be a characteristic of a solution in nonlinear magnetically charged electrodynamics.

4.3 Third model of black hole

In constructing the third model, we assume again the symmetry provided by (48), and consider the following regular non-metricity scalar

$$\begin{aligned} {\mathbb {Q}}(r)=\frac{q_0}{(r+r_0)^2}. \end{aligned}$$
(84)

Thus, from Eq. (47), we have

$$\begin{aligned} -\frac{2 \Bigl (r A'(r)+A(r)\Bigr )}{r^2}=\frac{q_0}{(r+r_0)^2}, \end{aligned}$$
(85)

which yields the following solution

$$\begin{aligned} A(r)=-\frac{2M}{r}-\frac{q_{0}}{2r}\left[ r-\frac{r_{0}^{2}}{r+r_{0}} -2r_{0}\ln (r+r_{0})\right] . \end{aligned}$$
(86)

The metric function (86) diverges in the limit of \(r\rightarrow 0\), but becomes asymptotically flat at spatial infinity, \(r\rightarrow \infty \).

Here, too, we analyze the presence of horizons and the critical parameter \(q_{0c}\) using Eqs. (55) and (56), based on numerical solutions. Using these equations and the values for the constants \(r_0=0.1\) and \(M=1\), we find the critical parameter value \(q_{0c}=-55.475\). Under these assumptions, we plot the metric function (86) in Fig. 7 for three scenarios: \(q_0 > q_{0c}\), \(q_0 = q_{0c}\), and \(q_0 < q_{0c}\). When \(q_0\) exceeds the critical value, \(q_0 > q_{0c}\), we observe the presence of an horizon. For \(q_0 = q_{0c}\), we observe a light-type singularity, and if \(q_0 < q_{0c}\), no horizon is formed.

Fig. 7
figure 7

The plot depicts the metric function (53), i.e., A(r), for the specific values of \(r_0=0.1\) and \(M =5.0\), the critical value of the constant \(q_0\), given by \(q_{0c} = -55.475\). For \(q_0 > q_{0c}\) there is only one horizon, for \(q_0 = q_{0c}\), there is light-type singularity and for \(q_0 < q_{0c}\), there is no horizon

Full size image

In this scenario, the boundary term is given by

$$\begin{aligned} {\mathbb {B}}_Q(r)=&\frac{2 (q_0+1) r^3+3 (q_0+2) r^2 r_0+6 r r_0^2+2 r_0^3}{r^2 (r+r_0)^3}, \end{aligned}$$
(87)

so that the scalar (84) and boundary term (87) satisfy the condition (13). In this model, we also find that the combination \({\mathbb {Q}}(r)-{\mathbb {B}}_Q(r)\) is singular in the limit of \(r\rightarrow 0\).

Now, Eqs. (49) and (50) are expressed as

$$\begin{aligned}&{{\mathcal {L}}}_{\textrm{NLED}}(r)\nonumber \\&=-\frac{1}{2\kappa ^{2}} \left[ \frac{q_{0}r^{3}r_{0}(r+r_{0})f'(r)}{4(q_{0}+1) r^{4} +(7q_{0}+16)r^{3}r_{0}+24r^{2}r_{0}^{2} +16rr_{0}^{3}+4r_{0}^{4}}\right. \nonumber \\ {}&\left. \qquad +f(r)\right] , \end{aligned}$$
(88)
$$\begin{aligned}&{{\mathcal {L}}}_F(r)\nonumber \\&=\frac{r^{5}(r+r_{0})f'(r) \Bigl [(q_{0}+2)r^{3}+6r^{2}r_{0} +6rr_{0}^{2} +2r_{0}^{3}\Bigr ]}{2\kappa ^{2}q^{2} \Bigl [4(q_{0}+1)r^{4} +(7q_{0}+16)r^{3}r_{0} +24r^{2}r_{0}^{2}+16rr_{0}^{3}+4r_{0}^{4}\Bigr ]}. \end{aligned}$$
(89)

Using Eqs. (88) and (89) in the relation (59), allows us to determine

$$\begin{aligned}{} & {} f(r)=\frac{f_{0}}{r_{0}^{2}}\left[ 4e^{2}q_{0} \text {Ei}\left( -\frac{2(r+r_{0})}{r_{0}}\right) \right. +8\text {Ei}\left( -\frac{2r}{r_{0}}\right) \nonumber \\{} & {} \left. +\frac{r_{0}e^{-\frac{2r}{r_{0}}} \left( 2(q_{0}+2)r^{4} +(3q_{0}+10)r^{3}r_{0}+6r^{2}r_{0}^{2}-2rr_{0}^{3}-2r_{0}^{4}\right) }{r^{2}(r+r_{0})^{3}}\right] \nonumber \\ {}{} & {} \qquad +f_1, \end{aligned}$$
(90)

where \(f_0\) and \(f_1\) are integration constants, so that Eqs. (88) and (89) are given by

$$\begin{aligned} {{\mathcal {L}}}_{\textrm{NLED}} (r)= & {} \frac{1}{2\kappa ^{2}r_{0}^{2}}\left\{ -4e^{2}q_{0}\text {Ei}\left( -\frac{2(r+r_{0})}{r_{0}}\right) \right. -8\text {Ei}\left( -\frac{2r}{r_{0}}\right) \nonumber \\{} & {} \left. +\frac{r_{0}e^{-\frac{2r}{r_{0}}}\left[ -2(q_{0}+2)r^{3} -(q_{0}+6)r^{2}r_{0}+2r_{0}^{3}\right] }{r^{2}(r+r_{0})^{2}}\right\} ,\nonumber \\ \end{aligned}$$
(91)
$$\begin{aligned} {{\mathcal {L}}}_F (r) =\frac{r^{2}e^{-\frac{2r}{r_{0}}}\Bigl [(q_{0}+2)r^{3} +6r^{2}r_{0}+6rr_{0}^{2}+2r_{0}^{3}\Bigr ]}{2\kappa ^{2}q^{2}(r+r_{0})^{3}}, \end{aligned}$$
(92)

respectively.

Finally, using Eq. (30), we derive the expression for r(F), so that

$$\begin{aligned}&{{\mathcal {L}}}_{\textrm{NLED}} (F) \nonumber \\&= \frac{e^{-\frac{2^{3/4} \sqrt{q}}{\root 4 \of {F} r_0}} \left( -2\ 2^{3/4} \root 4 \of {F} q^{3/2} (q_0+2)-2 \sqrt{F} q (q_0+6) r_0+4 \sqrt{2} F r_0^3\right) }{\kappa ^2 q r_0 \left( 2 \root 4 \of {F} r_0+2^{3/4} \sqrt{q}\right) ^2} \nonumber \\&\quad -\frac{2 }{\kappa ^2 r_0^2} \left[ e^2 q_0 \text {Ei}\left( -\frac{2^{3/4} \sqrt{q}}{\root 4 \of {F} r_0}-2\right) \right. \nonumber \\&\quad \left. +2 \text {Ei}\left( -\frac{2^{3/4} \sqrt{q}}{\root 4 \of {F} r_0}\right) \right] , \end{aligned}$$
(93)

which is depicted in in Fig. 8 for different values of \(q_0\). We note that \({{\mathcal {L}}}_{\textrm{NLED}} (F)\) is non-linear and has the properties of an increasing bijective function. If we fix the values of \(r_0\) and q, we find different models for varying values of \(q_0\).

Fig. 8
figure 8

The plot depicts \({{\mathcal {L}}}_{\textrm{NLED}}(F)\), described by Eq. (93), for the values \(q_0=-0.1\), \(q_0=-1.0\) and \(q_0=-1.9\)

Full size image

The asymptotic limits of the Lagrangian (8) for \(F\gg 1\) are

$$\begin{aligned} {{\mathcal {L}}} _{\textrm{NLED}} (F)&\sim 0.01\frac{e^{2}q_{0}}{\kappa ^{2}r_{0}^{2}} +\frac{\sqrt{2}\sqrt{F}}{\kappa ^{2}q} -\frac{4\root 4 \of {2}\root 4 \of {F}}{\kappa ^{2} \sqrt{q}r_{0}}\nonumber \\&\quad +\frac{\root 4 \of {\frac{1}{F}}\sqrt{q}(3q_{0}+8)}{3\root 4 \of {2}\kappa ^{2}r_{0}^{3}}\nonumber \\&\quad +\frac{2\ln (F) -q_{0}-8\gamma +12}{2\kappa ^{2}r_{0}^{2}}\nonumber \\&\quad -\frac{2\left( \frac{1}{F}\right) ^{3/4}\left( 15\root 4 \of {2}q^{3/2}q_{0} -2\root 4 \of {2}q^{3/2}\right) }{45\left( \kappa ^{2}r_{0}^{5}\right) } \nonumber \\&\quad +\frac{375q^{2}q_{0}-4q^{2}}{180F\kappa ^{2}r_{0}^{6}} \nonumber \\&\quad -\frac{\sqrt{\frac{1}{F}}(q(3q_{0}+4))}{6\left( \sqrt{2}\kappa ^{2}r_{0}^{4} \right) }, \end{aligned}$$
(94)

where \(\gamma \) is the Euler constant. The numerical value 0.01 came from an approximation of the integral exponential function. We illustrate again the behavior described by the Lagrangian of our model, Eq. (93), with a Lagrangian that is proportional to the linear case, represented by the blue curve and the dashed red line, respectively, as described in Fig. 9.

Fig. 9
figure 9

Graphic representation of \({{\mathcal {L}}}_{\textrm{NLED}}(F)\), described by expression (93), with respect to F. Where we consider \(q=0.7\), \(q_0=-1.0\) and \(r_0=0.5\). The blue curve represents the behavior of our model for \(F\ll 1\), while the dashed red line represents the linear case

Full size image
Fig. 10
figure 10

The left plot depicts the graphical representation of the expression \({\text {Ln}}[f_Q(r)]\), described by Eq. (102), with respect to the coordinate r. The right plot represents the graphic representation of \({\text {Ln}}[f_Q({\mathbb {Q}})]\) with respect to the boundary term \({\mathbb {Q}}\). We have considered \(M=10\) and \(q=1\)

Full size image

5 Regular black hole solutions in \(f(\mathbb {Q,B_Q})\) gravity

In this section, we explore solutions for regular black holes, taking into account the metric functions of the Bardeen [42] and Culetu solutions [59, 60]. To obtain these solutions, we use again the symmetry imposed in Eq. (48), and the line element (45).

5.1 Bardeen-type solutions

For the first solution, consider the metric function described by Bardeen’s model [42], given by

$$\begin{aligned} A(r)=1-\frac{2Mr^{2}}{(q^{2}+r^{2})^{3/2}}, \end{aligned}$$
(95)

so that Eqs. (38) and (39) take the following forms:

$$\begin{aligned} {\mathbb {B}}_Q(r)= & {} \frac{6 M q^2 r^2}{\left( q^2+r^2\right) ^{7/2}}\nonumber \\{} & {} +\frac{36 M q^4}{\left( q^2+r^2\right) ^{7/2}}-\frac{2}{r^2}, \end{aligned}$$
(96)
$$\begin{aligned} {\mathbb {Q}}(r)= & {} \frac{12 M q^2}{\left( q^2+r^2\right) ^{5/2}}-\frac{2}{r^2}, \end{aligned}$$
(97)

respectively. We verify that in the limit of \(r\rightarrow 0\) the combination of \({\mathbb {Q}}(r)-{\mathbb {B}}_Q(r)\) from Eqs. (96) and (97) is regular, so that the spacetime geometry is regular.

The following electromagnetic relations are given by

$$\begin{aligned}{} & {} {{\mathcal {L}}}_{\textrm{NLED}} (r) =\frac{1}{2 \kappa ^2} \Biggl [\frac{3 M q^2 r^3 \left( 3 r^2-2 q^2\right) f'(r)}{\left( q^2+r^2\right) ^{7/2}-15 M q^2 r^4}-2 f(r)\Biggr ], \nonumber \\ \end{aligned}$$
(98)
$$\begin{aligned}{} & {} {{\mathcal {L}}}_F (r) =\frac{15 M r^9 f'(r)}{4 \kappa ^2 \left[ \left( q^2+r^2\right) ^{7/2}-15 M q^2 r^4\right] }, \end{aligned}$$
(99)

so that from the consistency relation (59), we obtain the following expression:

$$\begin{aligned} f(r)=f_{0}\int \left( \int \exp \left( \zeta \right) dr\right) dr+f_{1}, \end{aligned}$$
(100)

where \(f_0\) and \(f_1\) are constants, to simplify this expression and subsequent ones, we define

$$\begin{aligned} \zeta= & {} -\frac{1}{3Mq^{2}r^{3}\left( 2q^{2} -3r^{2}\right) \left( \left( q^{2}+r^{2}\right) ^{7/2}-15Mq^{2}r^{4}\right) } \nonumber \\{} & {} \times \Biggl [2q^{14}+14q^{12}r^{2}+2r^{14}+ \nonumber \\{} & {} +q^{8}\left( 70r^{6}-123Mr^{4}\sqrt{q^{2}+r^{2}}\right) \nonumber \\{} & {} +3q^{4}r^{8}\left( 345M^{2}-101M\sqrt{q^{2}+r^{2}}+14r^{2}\right) \nonumber \\{} & {} +2q^{2}r^{10}\left( 7r^{2}-36M\sqrt{q^{2}+r^{2}}\right) \nonumber \\{} & {} +6q^{10}\left( 3Mr^{2}\sqrt{q^{2}+r^{2}} +7r^{4}\right) \nonumber \\{} & {} +2q^{6}r^{6}\left( 45M^{2}-186M\sqrt{q^{2}+r^{2}}+35r^{2} \right) \Biggl ]. \end{aligned}$$
(101)

From Eq. (100), we obtain the derivative of \(f_Q(r)\), given by

$$\begin{aligned} f_{Q}(r)= & {} \frac{r^3 \left( q^2+r^2\right) ^{7/2} }{4 \left( \left( q^2+r^2\right) ^{7/2}-15 M q^2 r^4\right) }\nonumber \\{} & {} \times \int \exp \left( \zeta \right) \, dr. \end{aligned}$$
(102)

We depict the linearity of the function \(f_Q(r)\), described by Eq. (102), with respect to the non-metricity scalar in two ways in Fig. 10. On the left, the graph illustrates the variation of the natural logarithm of \(f_Q(r)\) with respect to r, highlighting a non-linear dependence on \({\mathbb {Q}}\). On the right, it is evident that the function \(f_Q({\mathbb {Q}})\) is non-linear with respect to \({\mathbb {Q}}\). Furthermore, we observe that the function \(f({\mathbb {Q}},{\mathbb {B}}_Q)\) is also non-linear with respect to \({\mathbb {Q}}\) (Fig. 11).

With these expressions we find that the Lagrangian density is now given by

$$\begin{aligned} {{\mathcal {L}}} _{\textrm{NLED}}(r)= & {} \frac{1}{4\kappa ^{2}} \left[ \frac{M\left( 9q^{2}r^{5}-6q^{4}r^{3}\right) }{\left( q^{2}+r^{2}\right) ^{7/2}-15Mq^{2}r^{4}}\right. \nonumber \\{} & {} \times \left. \int \exp \left( \zeta \right) \,dr-2f(r)\right] , \end{aligned}$$
(103)

where f(r) in the above equation is given by Eq. (100), and \({{\mathcal {L}}}_F(r)\) is

$$\begin{aligned} {{\mathcal {L}}}_F(r)= & {} \frac{15Mr^{9}}{4\kappa ^{2}\left[ \left( q^{2} +r^{2}\right) ^{7/2}-15Mq^{2}r^{4}\right] }\nonumber \\{} & {} \times \int \exp \left( \zeta \right) dr. \end{aligned}$$
(104)

Equation (104) is depicted in Fig. 12, where we verify that \({{\mathcal {L}}}_F(F)\) is a non-linear and bijective function with respect to F.

In Fig. 11 we show a comparison of the behavior of the Lagrangian obtained in our model (103) with a Lagrangian from the linear case.

Fig. 11
figure 11

Graphic representation of \({{\mathcal {L}}}_{\textrm{NLED}}(F)\), described by expression (103), with respect to F. We consider here \(M=10\), \(q=1\), \(a_1=2.45\times 10^7\) and \(a_2=-72.45\). The blue curve represents the behavior of our model for \(F\ll 1\), while the dashed red line represents the linear case, both described by the logarithmic function

Full size image
Fig. 12
figure 12

Graphic representation of \({{\mathcal {L}}}_F(F)\), described by Eq. (104), with respect to F, where we have considered \(M=10\) and \(q=1\)

Full size image

5.2 Culetu-type solutions

In this second approach, we use the model proposed by Culetu [59, 60]:

$$\begin{aligned} A(r)=1-\frac{2 M }{r} \exp \left( -\frac{q^2}{2 M r}\right) , \end{aligned}$$
(105)

so that Eqs. (38) and (39) take the following forms:

$$\begin{aligned}{} & {} {\mathbb {B}}_Q(r)= \frac{e^{-\frac{q^2}{2 M r}} \left( -4 M r^3 e^{\frac{q^2}{2 M r}}+4 M q^2 r+q^4\right) }{2 M r^5}, \end{aligned}$$
(106)
$$\begin{aligned}{} & {} {\mathbb {Q}}(r)=\frac{2 q^2 e^{-\frac{q^2}{2 M r}}-2 r^2}{r^4}, \end{aligned}$$
(107)

respectively. We verify that, in the limit \(r \rightarrow 0\), the combination of \({\mathbb {Q}}(r) - {\mathbb {B}}_Q(r)\) from Eqs. (106) and (107) is regular, rendering a regular spacetime geometry.

We also deduce the following relations

$$\begin{aligned}{} & {} {{\mathcal {L}}}_{\textrm{NLED}}(r)=-\frac{1}{2 \kappa ^2} \Biggl (f(r)\nonumber \\{} & {} \quad -\frac{2 M q^2 r^2 f'(r) \left( q^2-4 M r\right) }{16 M^2 r^4 e^{\frac{q^2}{2 M r}}-32 M^2 q^2 r^2-6 M q^4 r+q^6}\Biggr ), \end{aligned}$$
(108)
$$\begin{aligned}{} & {} {{\mathcal {L}}}_F (r) =\frac{M r^6 f'(r) \left( q^2-8 M r\right) }{\kappa ^2 \left( 16 M^2 r^4 e^{\frac{q^2}{2 M r}}-32 M^2 q^2 r^2-6 M q^4 r+q^6\right) },\nonumber \\ \end{aligned}$$
(109)

which from the consistency condition (59), provides the following

$$\begin{aligned}{} & {} f(r)=\int \exp \nonumber \\{} & {} \left\{ -\frac{e^{2}q^{4}\text {Ei}\left( \frac{q^{2}}{2Mr}-2\right) +q^{4}\left[ \frac{4M}{r}-5\text {Ei} \left( \frac{q^{2}}{2Mr}\right) \right] +8Mre^{\frac{q^{2}}{2Mr}}\left( Mr+q^{2}\right) }{8M^{2}q^{2}}\right\} \nonumber \\{} & {} \times \frac{f_{0}\left( q^{2}-4Mr\right) ^{2} \left( -32M^{2}q^{2}r^{2}+16M^{2}r^{4}e^{\frac{q^{2}}{2Mr}} -6Mq^{4}r+q^{6}\right) }{r^{5}}dr +f_1,\nonumber \\ \end{aligned}$$
(110)

where \(f_0\) and \(f_1\) are constants.

Unfortunately, we were unable to obtain an analytic solution for the function f(r) directly from the consistency relation. However, we succeeded in determining the function \(f_B(r)\), which is described as

$$\begin{aligned}{} & {} f_B(r) =4 M^2 r^2 \left( q^2-4 M r\right) ^2 \nonumber \\{} & {} \times \exp \left\{ -\frac{q^4 \left[ e^2 \text {Ei}\left( \frac{q^2}{2 M r}-2\right) -5 \text {Ei}\left( \frac{q^2}{2 M r}\right) \right] +8 M r e^{\frac{q^2}{2 M r}} \left( M r+q^2\right) }{8 M^2 q^2}\right\} .\nonumber \\ \end{aligned}$$
(111)

The qualitative behavior of Eq. (111) with respect to the radial coordinate r is depicted in the left plot of Fig. 13, where the function \(f_B(r)\) is non-linear in relation to r, indicating non-linearity with respect to the boundary term \({\mathbb {B}}\). This can be directly verified in the right plot of Fig. 13, where the function \(f_B\) is clearly non-linear. Therefore, we conclude that the function \(f({\mathbb {Q}},{\mathbb {B}})\) must be non-linear in the boundary term \({\mathbb {B}}_Q\).

Fig. 13
figure 13

The left plot depicts \({\text {Ln}}[f_B(r)]\), described by expression (111), with respect to the radial coordinate r. The right plot represents \({\text {Ln}}[f_B(r)]\) with respect to the boundary term \({\mathbb {B}}_Q\). We consider \(M=10\) and \(q=1\)

Full size image

Taking into account Eqs. (110) and (111), we find that the Lagrangian density is now given by

$$\begin{aligned}{} & {} {{\mathcal {L}}} _{\textrm{NLED}}(r)\nonumber \\{} & {} \quad =-\frac{1}{2\kappa ^{2}}\left[ -\exp \left\{ -\frac{q^{2}\left[ -5r\text {Ei}\left( \frac{q^{2}}{2Mr}\right) +e^{2}r\text {Ei}\left( \frac{q^{2}}{2Mr}-2\right) +4M\right] }{8M^{2}r}\right. \right. \nonumber \\{} & {} \quad \left. -\frac{re^{\frac{q^{2}}{2Mr}}\left( Mr+q^{2}\right) }{Mq^{2}}\right\} \nonumber \\{} & {} \quad \left. \times \frac{2Mq^{2}\left( q^{2}-4Mr\right) ^{3}}{r^{3}}+f(r)\right] \end{aligned}$$
(112)

and its derivative is

$$\begin{aligned} {{\mathcal {L}}}_F (r)= & {} \frac{M r \left( q^2-8 M r\right) \left( q^2-4 M r\right) ^2 }{\kappa ^2}\nonumber \\{} & {} \times \exp \left[ -\frac{q^2 \left( -5 r \text {Ei}\left( \frac{q^2}{2 M r}\right) +e^2 r \text {Ei}\left( \frac{q^2}{2 M r}-2\right) +4 M\right) }{8 M^2 r}\right. \nonumber \\{} & {} \quad \left. -\frac{r e^{\frac{q^2}{2 M r}} \left( M r+q^2\right) }{M q^2}\right] . \end{aligned}$$
(113)

From Eq. (30), we obain r(F), which allows us to write \({{\mathcal {L}}}_F\) in terms of F, which is given by

$$\begin{aligned} {{\mathcal {L}}}_F(F)= & {} \exp \left\{ -\frac{q^{2}}{8M^{2}} \left[ e^{2}\text {Ei} \left( \frac{\root 4 \of {F}q^{3/2}}{2^{3/4}M}-2\right) \right. \right. \nonumber \\{} & {} \left. -5\text {Ei}\left( \frac{\root 4 \of {F}q^{3/2}}{2^{3/4}M}\right) \right] \nonumber \\{} & {} +\frac{\root 4 \of {2F}q^{3/2}}{2M} \nonumber \\{} & {} \left. +\frac{\root 4 \of {2F}q^{3/2}+M}{\sqrt{2F}Mq}\; e^{\frac{\root 4 \of {F}q^{3/2}}{2^{3/4}M}}\right\} \nonumber \\{} & {} \times \frac{M}{2F\kappa ^{2}} \left( \root 4 \of {2F}q^{3/2}-8M\right) \left( \root 4 \of {2F}q^{5/2}-4Mq\right) ^{2}.\nonumber \\ \end{aligned}$$
(114)

We present the behavior of Eq. (114) in Fig. 14, which depicts \({{\mathcal {L}}}_F(F)\) as a non-linear bijective function with respect to F.

For a more comprehensive understanding, we depict the asymptotic limits of the Lagrangian (112) in two distinct graphs. In Fig. 15, we illustrate the behavior of \({{\mathcal {L}}} (F)\) for \(F\ll 1\). In this context, we compare our Culetu model with the linear case, depicted by the blue curve and the dashed red line, respectively. Meanwhile, we explore the behavior of \({{\mathcal {L}}} (F)\) for \(F\gg 1\), once again comparing it with the linear case.

6 Summary and conclusion

In this work we investigated solutions for black holes and regular black holes in the newly proposed theory of \(f({\mathbb {Q}},{\mathbb {B}}_Q)\) gravity. This theory is a generalization of the \(f({\mathbb {Q}})\) theory, which describes gravitational effects through the non-metricity tensor \(Q_{\beta \mu \nu }\). More specifically, we couple \(f({\mathbb {Q}},{\mathbb {B}}_Q)\) gravity with NLED and find generalizations of solutions for black holes and regular black holes. In particular, we thus further extend the class of regular solutions for models of the Bardeen and Culetu type black holes found in the literature. It is interesting to note that in the \(f({\mathbb {Q}})\) theory, in a coincident gauge, a constraint was found that renders the theory linear [52], while in \(f({\mathbb {Q}},{\mathbb {B}}_Q)\) gravity, we gain a certain freedom which is not necessarily linear.

In the first model, the solution (53) can be asymptotically Minkowskian if \(b_0=-2\). However, it is also possible to choose values for \(b_0\) that are smaller than \(-2\), but in this case, the metric function (53) remains asymptotically flat and the metric signature becomes (\(-,+,-,-\)). This solution also generalizes GR with a term proportional to \(1/r^4\). When \(b_0=-2\) and \(b_1=0\), we arrive at the Schwarzschild solution. However, if \(b_1\ne 0\), we get a different geometry that is still asymptotically flat. Therefore, it is interesting to study the physical properties of this solution, such as shadows and gravitational lensing, for example. Although the non-metricity scalar in Eq. (54) is regular, the combination \({\mathbb {Q}}-{\mathbb {B}}_Q\) is singular, rendering a singular spacetime geometry.

Fig. 14
figure 14

Graphical representation of the derivative of the electromagnetic Lagrangian \({{\mathcal {L}}} _F(F)\), described by expression (114), as a function of F. We consider here \(M=10\) and \(q=1\)

Full size image
Fig. 15
figure 15

Graphic representation of \({{\mathcal {L}}}(F)\), described by expression (112). We consider here \(M=10\) and \(q=1\). The blue curve represents the behavior of our model for \(F\ll 1\), while the dashed red line represents the linear case

Full size image

In the second model, described by the metric function (73), the latter diverges in the limit \(r\rightarrow 0\), while, as \(r\rightarrow \infty \), A(r) becomes asymptotically flat. If we choose \(q_0=-2\), the metric function A(r) becomes asymptotically Minkowskian. Although the non-metricity scalar is regular Eq. (71), this does not guarantee that the spacetime solution is regular, as the combination of \({\mathbb {Q}}\) and \({\mathbb {B}}_Q\) for this model ensures that our solution is singular.

The equations of motion of the \(f({\mathbb {Q}},{\mathbb {B}}_Q )\) theory are quite complicated, coupled non-linear differential equations. Taking the vacuum as a possibility, there are two ways to proceed, one is to consider a functional form for the function \(f({\mathbb {Q}},{\mathbb {B}}_Q )\), and then integrate the equations to obtain the metric functions, first fixing a gauge. Another is to leave \(f({\mathbb {Q}},{\mathbb {B}}_Q )\) generic, fix a gauge, and integrate the equations to determine the metric functions and the functional form of \(f({\mathbb {Q}},{\mathbb {B}}_Q )\) together. These two paths seem impractical due to the complexity of the equations. Our strategy was to couple an NLED Lagrangian, leave the functions \(\mathcal {L}\) and \(\mathcal {L}_F\) free, solve the equations to determine these functions with respect to \(f({\mathbb {Q}},{\mathbb {B}}_Q )\) and the metric functions, and then set up the consistency equation such that \(\mathcal {L}_F\) is the derivative of \(\mathcal {L}\) with respect to F. Then we define a model for the metric functions, integrate the consistency equation to determine \(f({\mathbb {Q}},{\mathbb {B}}_Q )\), and finally the functional form of \(\mathcal {L}(F)\), if possible. This procedure has already been successfully applied to several other modified gravitational theories [46, 47, 61,62,63,64]. The solutions obtained are generally of an NLED that does not become linear for the weak field limit, i.e. \(F\ll 1\), but this also happens for some solutions in general relativity, such as the first one suggested by Bardeen [43].

In the third black hole model, we obtain the metric function described by Eq. (86). In the limit of \(r\rightarrow 0\), this metric function diverges; however, if in the limit of \(r\rightarrow \infty \), we find that the metric function is asymptotically flat. However, if we choose \(q_0=-2\), this metric function becomes asymptotically Minkowskian. Again, we find that the combination of the quantities \({\mathbb {Q}}\) and \({\mathbb {B}}_Q\) is singular, rendering a singular geometry. Analogously to regular black hole solutions that exhibit a regular spacetime, the derivative of the Lagrangian, \({{\mathcal {L}}}_F\), for both models of regular black holes is also distinct from its counterparts in GR.

Considering that the theory addressed in this work is of recent development, we see in this proposal a significant opportunity to explore new approaches. Our intention is to continue the study of \(f({\mathbb {Q}},{\mathbb {B}}_Q)\), exploring alternatives similar to those analysed in other theories, such as studies of black hole thermodynamics in f(R) gravity [65,66,67], f(T) gravity [68,69,70,71,72] and \(f(T,{\tilde{B}})\) gravitational theories [73], perturbation theory in f(R) [74,75,76], f(R, T) [77, 78], f(T) [79, 80] and \(f({\mathbb {Q}})\) gravities [51, 81], as well as the study of black hole shadows and the analysis of gravitational lensing in f(R) [82,83,84,85], f(R, T) [86] and f(T) gravity [87]. Furthermore, we intend to explore new solutions of black holes and regular black holes in the context of nonlinear electrodynamics and with scalar field, through other classes of connections, i.e., beyond the coincident gauge, as presented and applied in [39, 53, 88] These are some of the areas we intend to address in the future to deepen our understanding and contribute to the development of this new modified theory of gravity.