1 Introduction

The use of differential three-form fields in the realm of cosmology has been gaining ever more attention over the last decade [1]. These form fields naturally emerge in fundamental theories, such as string theory [2,3,4] and, thus it is only reasonable to explore their existence under effective formulations of gravity. Its application in cosmology has already proven to be fruitful in explaining the early- and late-time acceleration periods of the cosmic history [5,6,7,8], reheating [9], screening solutions [10], generation of cosmological magnetic fields [11], among others. Primordial inflation driven by multiple three-form fields was studied in [12], considering several potential functions. One appealing consequence is that these models present distinct signatures when compared to the standard inflationary setting with a scalar degree of freedom, compatible with recent cosmological observations. Inflationary models in extra dimensional braneworld scenarios inhabited by a single three-form have also been explored in [13], through the use of dynamical systems analysis, and tested against the Planck data. The computation of non-Gaussianities produced by several three-form fields inflation has been examined in [14] through the analysis of curvature perturbations employing the \(\delta N\) formalism.

It is well known that in four spacetime dimensions a three-form field admits a dual scalar field representation [15, 16]. For example, assuming nonquadratic three-form potentials leads to an equivalent scalar representation exhibiting a noncanonical kinetic term. However, this mapping is nontrivial, and thus for several self-interaction choices, or for any nonminimal coupling, this dual representation breaks down [5]. Nonetheless, even in the cases where the dual scalar description exists, it is often quite complex to deal with and it becomes much more practical and intuitive to work in the form-representative framework. An interesting feature of three-forms, is that the standard Maxwell term, \(F=d A\), constructed from a massless three-form A, naturally induces a cosmological constant term, and therefore has been used to address the cosmological puzzle regarding the tiny value of \(\Lambda \) [17]. Furthermore, nonminimal interactions between dark energy, driven by a three-form field, and cold dark matter, were explored in [18, 19] along with the corresponding linear cosmological perturbations. It was shown that even small values for the coupling lead to substantial variations on the growth of matter fluctuations.

Screening solutions with three-form fields conformally coupled to matter were also studied in [10]. In [20] the authors investigate the existence of future abrupt cosmological events, particularly the little sibling of the Big Rip, and the possibility of avoiding these by considering interactions between the three-form field, portraying dark energy, and dark matter. With the aid of dynamical systems techniques, it was found that this can be achieved only by considering interactions not directly involving the dark matter species. The authors also shed some light on how to distinguish the quadratic and linear dark energy interactions, through the statefinder hierarchy diagnosis and computing the growth of matter perturbations, compatible with the observational SDSS III data. An alternate procedure to avoid these future cosmological abrupt events induced by the three-form is the quantization of the aforementioned system [21]. The Lagrangian formalism for a single three-form fluid with a noncanonical Maxwell term was also explored, in comparison with k-essence cosmology, in [22].

More recently, theories embracing three-forms have been extended to spherically symmetric and static spacetimes, in particular, wormhole geometries [23]. More specifically, solutions were found for the modified field equations, assuming a single static and radial-dependent three-form field, where the standard matter fields are allowed one to dwell within the entire wormhole domain without violating the null and weak energy conditions.

Indeed, the present work also deals with static three-forms supporting spherically symmetric spacetimes, however, in the context of black holes and naked singularities. In fact, black hole solutions are well known in many gravitational field models and, in particular, in standard scalar–tensor extensions of general relativity. For instance, black hole solutions were recently found in the scalar–tensor representation of the hybrid metric-Palatini gravitational theory [24], which is a combination of the metric and Palatini f(R) formalisms unifying local constraints at the Solar System level and the late-time cosmic acceleration [25,26,27,28]. Furthermore, many other exact analytical black hole solutions have been obtained and studied extensively for nonminimally coupled scalar fields [29,30,31,32,33,34] (for a review of the nonsingular general relativistic solutions with minimally coupled scalar fields see [35]).

Generally the latter solutions have been derived in the Einstein frame, without the assumption of the existence of any coupling between the scalar field and the Ricci scalar. Similarly to these scalar field models, the three-form field theory we have considered is also formulated in the Einstein frame. However, there are fundamental differences between the three-form field theory and scalar field models in the Einstein or Jordan frames. An interesting result in the Brans–Dicke-type scalar–tensor theories is that the solutions with zero scalar field potential have in general no horizons [36]. Our analytic and numerical investigations show that this is not the case in the three-form field theory. Another interesting result in scalar–tensor theories is that globally regular, asymptotically flat solutions are possible. These solutions correspond to at least partly negative potentials \(V(\phi )\), and they are solitons without horizons and with a regular center [37]. However, these results specific to scalar–tensor theories cannot be recovered in the three-form fields gravitational theory. On the other hand our analytical and numerical investigations do not indicate the possible existence of any globally regular solutions.

This work is outlined in the following manner: in Sect. 2, we present the general formalism of the three-form field extension of general relativity. In Sect. 3, considering a static and spherically symmetric background, we deduce the gravitational field equations. In Sect. 4, exact vacuum solutions with three-form fields are presented, for the specific cases of a zero potential and for a constant potential. In Sect. 5, numerical solutions of the field equations are explored, for the Higgs potential and the exponential potential. The thermodynamic properties of the black holes solutions obtained are studied in Sect. 6. Finally, in Sect. 7, we discuss and summarize our results.

2 Einstein gravity with a three-form field: general formalism

We start by considering the following action for Einstein gravity with a standard three-form field \(A_{\alpha \beta \gamma }\):

$$\begin{aligned} \mathcal {S}=\int d^4 x \sqrt{-g}\left( \frac{1}{2\kappa ^2}R+\mathcal {L}_A \right) , \end{aligned}$$
(1)

with R being the Ricci scalar, \(g=\) det \(g_{\mu \nu }\) the determinant of the metric tensor, \(\kappa ^2 = 8\pi G\) and \(\mathcal {L}_A\) stands for the Lagrangian density for our three-form, which reads [1, 5]

$$\begin{aligned} \mathcal {L}_A=-\frac{1}{48}F^2 - V(A^2), \end{aligned}$$
(2)

where we have used the notation

$$\begin{aligned} F^2 = F_{\alpha \beta \gamma \delta }F^{\alpha \beta \gamma \delta }\quad \quad \mathrm{and}\quad \quad A^2 = A_{\alpha \beta \gamma }A^{\alpha \beta \gamma }. \end{aligned}$$
(3)

Here \(V(A^2)\) is the three-form potential and \(\mathbf{F}=\mathbf{dA}\) is the field strength tensor [15, 20], a four-form, whose components can be written as

$$\begin{aligned} F_{\alpha \beta \gamma \delta } = \nabla _{\alpha }A_{\beta \gamma \delta }-\nabla _{\delta }A_{\alpha \beta \gamma }+\nabla _{\gamma }A_{\delta \alpha \beta }-\nabla _{\beta }A_{\gamma \delta \alpha }, \end{aligned}$$
(4)

with \(\nabla _{\mu }\) being the covariant derivative. The equations of motion for our three-form can be found by varying the action Eq. (1) with respect to \(A_{\alpha \beta \gamma }\). They read [6, 23]

$$\begin{aligned} \nabla _{\mu } F^{\mu }_{\,\,\,\,\,\,\alpha \beta \gamma } =12\frac{\partial V}{\partial (A^2)}A_{\alpha \beta \gamma }. \end{aligned}$$
(5)

Varying now Eq. (1) with respect to the metric \(g^{\mu \nu }\) one finds the following field equations:

$$\begin{aligned} G_{\mu \nu } =\kappa ^2 \, T_{\mu \nu }, \end{aligned}$$
(6)

where \(G_{\mu \nu }\) is the Einstein tensor and \(T_{\mu \nu }\) is the energy momentum tensor of our form field source:

$$\begin{aligned} T_{\mu \nu }= & {} -2 \frac{\delta \mathcal {L}_{A}}{\delta g^{\mu \nu }}+g_{\mu \nu }\mathcal {L}_{A} \nonumber \\= & {} \frac{1}{6}\left( F\circ F \right) _{\mu \nu } + 6 \frac{\partial V}{\partial (A^2)} \left( A\circ A \right) _{\mu \nu } +\mathcal {L}_A\,g_{\mu \nu }, \end{aligned}$$
(7)

with a circle denoting contraction of all indices but the first, i.e., \(\left( F\circ F \right) _{\mu \nu }=F_{\mu \alpha \beta \gamma }F_{\nu }^{\,\,\,\,\alpha \beta \gamma }\) .

We will construct our three-form with the aid of its dual vector (one-form) [23], via the Hodge star operator

$$\begin{aligned} B^{\mu } = (\star A)^{\delta } = \frac{1}{3!}\frac{1}{\sqrt{-g}}\,\epsilon ^{\mu \alpha \beta \gamma }A_{\alpha \beta \gamma }. \end{aligned}$$
(8)

Inverting the last identity, we express the three-form components in terms of its dual vector

$$\begin{aligned} A_{\alpha \beta \gamma } = \sqrt{-g}\,\epsilon _{\alpha \beta \gamma \delta }B^{\delta }. \end{aligned}$$
(9)

We now build the three-form components by parameterizing \(B^{\mu }\) in terms of a radial scalar function \(\zeta = \zeta (r)\), i.e.,

$$\begin{aligned} B^{\delta } = \left( 0,\zeta (r),0,0 \right) ^\mathrm{T}. \end{aligned}$$
(10)

Due to the antisymmetric nature of differential forms, once we attain a solution for \(\zeta (r)\), all the three-form components are automatically determined through Eq. (9).

Through Eqs. (4) and (9), we may rewrite the kinetic term in the action Eq. (1) as [7]

$$\begin{aligned} -\frac{1}{48}F^2=-\frac{1}{2} F_{0123}F^{0123} = \frac{1}{2}(\nabla _{\mu }B^{\mu })^2. \end{aligned}$$
(11)

We now consider applications of the general formalism obtained here to the specific case of static and spherically symmetric spacetimes.

3 Spherically symmetric and static background

3.1 Metric and field equations

Consider a static and spherically symmetric spacetime, given by the following line element:

$$\begin{aligned} ds^2 = -e^{\alpha (r)}dt^2 + e^{\beta (r)}dr^2 + r^2 \left( d\theta ^2 + \sin ^2\theta \,d\phi ^2 \right) .\nonumber \\ \end{aligned}$$
(12)

On this background geometry, the invariant \(A^2\) is given by

$$\begin{aligned} A^2 = -6\, e^{\beta (r)} \zeta (r)^2, \end{aligned}$$
(13)

and the term expressing the kinetic energy of the three-form is provided by

$$\begin{aligned} F^2 = -6\left[ \zeta \left( \alpha ' + \beta ' +\frac{4}{r} \right) + 2\zeta ' \right] ^2, \end{aligned}$$
(14)

where a prime denotes the derivative with respect to the radial component.

The equations of motion (5) can now be written in terms of \(\zeta \), using the metric Eq. (12), as

$$\begin{aligned} 2\zeta '' + \left( \alpha '+\beta '+ \frac{4}{r}\right) \zeta '+\left( \alpha ''+\beta ''-\frac{4}{r^2} \right) \zeta +2V_{,\zeta }=0,\nonumber \\ \end{aligned}$$
(15)

where \(V_{,\zeta } = \partial V / \partial \zeta \).

Through Eq. (7), using Eq. (12), the components of the energy momentum tensor are then found to be

$$\begin{aligned}&T^t_{\,\,\,\,t} = -\rho = \frac{F^2}{48} - V + \zeta V_{,\zeta }, \end{aligned}$$
(16)
$$\begin{aligned}&T^r_{\,\,\,\,r} = p_r = \frac{F^2}{48} - V, \end{aligned}$$
(17)
$$\begin{aligned}&T^{\theta }_{\,\,\,\,\theta } = T^{\phi }_{\,\,\,\,\phi } = p = T^t_{\,\,\,\,t}, \end{aligned}$$
(18)

where the \(F^2\) term is given by Eq. (14) and we identify \(\rho \), \(p_r\) and p with the energy density, the radial pressure, and the tangential pressure, respectively, of the three-form. The components of the Einstein tensor can be written as

$$\begin{aligned} G^t_{\,\,\,\,t}= & {} \frac{1}{r^2}\frac{{\text {d}}}{{\text {d}}r}\left[ r\left( e^{-\beta }-1 \right) \right] = \frac{e^{-\beta }}{r^2}\left( 1-r\beta ' - e^{\beta } \right) , \end{aligned}$$
(19)
$$\begin{aligned} G^r_{\,\,\,\,r}= & {} \frac{e^{-\beta }}{r^2} \left( 1+r \alpha '-e^{\beta }\right) , \end{aligned}$$
(20)
$$\begin{aligned} G^{\theta }_{\,\,\,\,\theta }= & {} G^{\phi }_{\,\,\,\,\phi } = \frac{e^{-\beta }}{2}\left[ \alpha '' + \left( \frac{1}{r}+\frac{\alpha '}{2} \right) \left( \alpha '-\beta ' \right) \right] . \end{aligned}$$
(21)

At this point, we have four variables, namely, \(\alpha \), \(\beta \), \(\zeta \) and V, and four independent equations, namely, the equation of motion for \(\zeta \), i.e., Eq. (15) and the three field equations (setting \(\kappa = 1\) for simplicity):

$$\begin{aligned}&\frac{e^{-\beta }}{r^2}\left( 1-r\beta ' - e^{\beta } \right) = \frac{F^2}{48} - V + \zeta V_{,\zeta }, \end{aligned}$$
(22)
$$\begin{aligned}&\frac{e^{-\beta }}{r^2} \left( 1+r \alpha '-e^{\beta }\right) = \frac{F^2}{48} - V, \end{aligned}$$
(23)
$$\begin{aligned}&\frac{e^{-\beta }}{2}\left[ \alpha '' + \left( \frac{1}{r}+\frac{\alpha '}{2} \right) \left( \alpha '-\beta ' \right) \right] = \frac{F^2}{48} - V + \zeta V_{,\zeta }.\nonumber \\ \end{aligned}$$
(24)

Combining the first two field equations, Eqs. (22) and (23), one finds

$$\begin{aligned} \alpha ' + \beta ' = -re^{\beta }\zeta V_{,\zeta }. \end{aligned}$$
(25)

Using Eqs. (22) and  (24) yields

$$\begin{aligned} \frac{2}{r^2}\left( 1-r \beta '-e^{\beta }\right) = \alpha '' + \left( \frac{1}{r}+\frac{\alpha '}{2} \right) \left( \alpha '-\beta ' \right) . \end{aligned}$$
(26)

Moreover, from Eq. (22) we obtain

$$\begin{aligned} r\beta ^{\prime }=1-e^{\beta }\left[ 1+\left( \frac{F^2}{48}-V+\zeta V_{\zeta }\right) r^2\right] . \end{aligned}$$
(27)

These expressions will be useful below.

3.2 Dynamical system formulation

The field equation (22) can be rewritten as

$$\begin{aligned} \frac{{\text {d}}}{{\text {d}}r}\left( re^{-\beta }\right) =1-\left( V-\zeta V_{,\zeta }-\frac{F^2}{48}\right) r^2, \end{aligned}$$
(28)

and can immediately be integrated to give

$$\begin{aligned} e^{-\beta }=1-\frac{2M_{\mathrm{eff}}(r)}{r}, \end{aligned}$$
(29)

where the effective mass \(M_{\mathrm{eff}}(r)\) is defined as

$$\begin{aligned} M_{\mathrm{eff}}(r)=\frac{1}{2}\int _0^r{\left( V-\zeta V_{,\zeta }-\frac{F^2}{48}\right) r^2{\text {d}}r}, \end{aligned}$$
(30)

and satisfies the mass continuity-type equation

$$\begin{aligned} \frac{{\text {d}}M_{\mathrm{eff}}(r)}{{\text {d}}r}=\frac{1}{2}\left( V-\zeta V_{,\zeta }-\frac{F^2}{48}\right) r^2. \end{aligned}$$
(31)

From Eq. (23) we obtain the expression of \(\alpha '\):

$$\begin{aligned} \alpha '=\frac{r^3\left( F^2/48-V\right) +2M_{\mathrm{eff}}(r)}{r^2\left( 1-2M_{\mathrm{eff}}(r)/r\right) }. \end{aligned}$$
(32)

With the use of Eq. (25), Eq. (15) can be reformulated as

$$\begin{aligned} \zeta ''+\left[ \frac{2}{r}-\frac{r\zeta \left( V_{,\zeta }+\zeta V_{,\zeta \zeta }/2\right) }{1-2M_{\mathrm{eff}}(r)/r}\right] \zeta '-G\left( r,\zeta \right) \zeta +V_{,\zeta }=0, \end{aligned}$$
(33)

where we have denoted

$$\begin{aligned} G\left( r,\zeta \right)= & {} \frac{2}{r^{2}}+\frac{\zeta V_{,\zeta }}{ 2\left[ 1-2M_{\mathrm{eff}}(r)/r\right] } \nonumber \\&\times \Bigg [ 2-\frac{1+\left( F^{2}/48-V+\zeta V_{,\zeta }\right) r^{2}}{1-2M_{\mathrm{eff}}(r)/r}\Bigg ] . \end{aligned}$$
(34)

Finally, for the function \(F^2\) we obtain

$$\begin{aligned} F^2=-6\left[ \left( \frac{4}{r}-\frac{r\zeta V_{,\zeta }}{1-2M_{\mathrm{eff}}(r)/r}\right) \zeta +2\zeta '\right] ^2. \end{aligned}$$
(35)

Here we have presented the relevant equations, which set the stage for exploring solutions, both exactly and numerically.

4 Exact vacuum solutions with three-form fields

In the present section we will consider some solutions of the system of the vacuum field equations (22)–(24), respectively, which must be solved together with Eq. (15), once the functional expression of the three-form potential V is fixed.

4.1 First case: \(V=0\)

If the three-form field potential V identically vanishes, \(V\equiv 0\), then Eq. (25) can be immediately integrated, giving

$$\begin{aligned} \alpha =-\beta , \end{aligned}$$
(36)

where we have set, without loss in generality, the arbitrary integration constant equal to zero. Then Eq. (26) becomes

$$\begin{aligned} \alpha ''+\alpha '^2=\frac{2}{r^2}\left( 1-e^{-\alpha }\right) . \end{aligned}$$
(37)

By introducing a new variable \(\alpha =\ln f\), Eq. (37) becomes

$$\begin{aligned} f''-\frac{2f}{r^2}+\frac{2}{r^2}=0, \end{aligned}$$
(38)

with the general solution given by

$$\begin{aligned} f\left( r\right) =1+\frac{c_1}{r}+c_2 r^2, \end{aligned}$$
(39)

where \(c_1\) and \(c_2\) are arbitrary constants of integration. Since in the limit \(F^2\rightarrow 0\) the Schwarzschild solution of standard general relativity must be recovered, it follows that \(c_1=-2M\), where M is the mass of the gravitating body, while \(c_2=-\Lambda \) can be interpreted as the cosmological constant. Hence, we have recovered the Schwarzschild–de Sitter solution, with the cosmological constant naturally included,

$$\begin{aligned} e^{\alpha }=e^{-\beta }=1-\frac{2M}{r}-\Lambda r^2. \end{aligned}$$
(40)

Equation (15), giving the evolution of the function \(\zeta \), becomes

$$\begin{aligned} \zeta ''+\frac{2}{r}\zeta ^{\prime }-\frac{2}{r^2}\zeta =0, \end{aligned}$$
(41)

and it has the general solution

$$\begin{aligned} \zeta (r)= C_1r + \frac{C_2}{r^2}, \end{aligned}$$
(42)

where \(C_1\) and \(C_2\) are arbitrary constants of integration. With the use of Eq. (42), we obtain finally for \(F^2\) the expression

$$\begin{aligned} F^2=-216 C_1^2=\mathrm{constant}, \end{aligned}$$
(43)

as expected. On the other hand with the use of the field equation (23) we obtain \(F^2=144c_1=-288M\). By comparing the two expressions for \(F^2\) we obtain for \(C_1\) the representation \(C_1=\sqrt{4M/3}\).

4.2 The constant potential \(V=V_0=\mathrm{constant}\)

In the case of the constant potential \(V=V_0=\mathrm{constant}\), Eq. (33) take the same form as in the case of the vanishing scalar potential of the three-form field, Eq. (41), and its solution is given again by Eq. (42). As for \(F^2\), given by Eq. (35), we obtain again \(F^2=-216C_1^2\). By integrating the mass continuity equation, Eq. (31), we obtain

$$\begin{aligned} M_{\mathrm{eff}}(r)=\frac{1}{6}\left[ V_0+\frac{9}{2}C_1^2\right] r^3+c_1, \end{aligned}$$
(44)

where \(c_1\) is an arbitrary constant of integration. Hence for \(e^{-\beta }\) we immediately obtain

$$\begin{aligned} e^{-\beta }=1-\frac{c_1}{r}-\frac{1}{6}\left[ V_0+\frac{9}{2}C_1^2\right] r^2. \end{aligned}$$
(45)

Since also in the constant potential case the general relation \(\alpha +\beta =0\) holds, we obtain the metric tensor coefficient \(e^{\alpha }\) as

$$\begin{aligned} e^{\alpha }=1-\frac{c_1}{r}-\frac{1}{6}\left[ V_0+\frac{9}{2}C_1^2\right] r^2. \end{aligned}$$
(46)

The solution is again of the Schwarzschild–de Sitter type, with the potential \(V_0\) generating, together with \(F^2\), an effective cosmological constant. On the other hand, the arbitrary integration constant \(c_1\) is undetermined by the field equations, and must be chosen from physical considerations. If the integration constant \(C_1\) and the constant potential \(V_0\) vanish, the metric reduces to the standard Schwarzschild form. Moreover, there are no restrictions on the integration constant \(c_1\), whose sign and physical interpretation remains arbitrary.

5 Numerical solutions of the field equations

The system of three Eqs. (31), (32) and (33) for the three unknown functions \(M_{\mathrm{eff}}\), \(\alpha \) and \(\zeta \), representing a strongly nonlinear system of differential equations, determines the vacuum solutions of the three-form field model gravity. In order to integrate the equations we introduce a new independent variable \(\eta \), defined as

$$\begin{aligned} \eta =\frac{1}{r}. \end{aligned}$$
(47)

Then we can reformulate the gravitational field equations as the following first order dynamical system:

$$\begin{aligned}&\frac{{\text {d}}\zeta }{{\text {d}}\eta }=u, \end{aligned}$$
(48)
$$\begin{aligned}&\frac{{\text {d}}M_{\mathrm{eff}}}{{\text {d}}\eta }=\frac{1}{2}\left( \frac{F^2}{48}+\zeta V_{,\zeta }-V\right) \frac{1}{\eta ^4}, \end{aligned}$$
(49)
$$\begin{aligned}&\frac{{\text {d}}\alpha }{{\text {d}}\eta }=-\frac{F^2/48-V+2\eta ^3M_{\mathrm{eff}}}{\eta ^3\left( 1-2\eta M_{\mathrm{eff}}\right) }, \end{aligned}$$
(50)
$$\begin{aligned}&\frac{{\text {d}}u}{{\text {d}}\eta }=-\frac{V_{,\zeta }+\zeta V_{,\zeta \zeta }/2}{\eta ^{3}\left( 1-2\eta M_{\mathrm{eff}} \right) }u+\frac{1}{\eta ^{4}}G\left( \eta ,\zeta \right) \zeta -\frac{V_{,\zeta }}{\eta ^{4}}, \end{aligned}$$
(51)

where

$$\begin{aligned} G\left( \eta ,\zeta \right)= & {} 2\eta ^{2}+\frac{\zeta V_{,\zeta }}{2\left( 1-2\eta M_{\mathrm{eff}}\right) } \nonumber \\&\times \left[ 2-\frac{\eta ^{2}+\left( F^{2}/48-V+\zeta V_{,\zeta }\right) }{\eta ^{2}\left( 1-2\eta M_{\mathrm{eff}} \right) }\right] \end{aligned}$$
(52)

and

$$\begin{aligned} F^{2}=-6\left[ \left( 4\eta -\frac{\zeta V_{,\zeta }}{\eta \left( 1-2\eta M_{\mathrm{eff}}\right) }\right) \zeta -2\eta ^2u\right] ^{2}, \end{aligned}$$
(53)

respectively. In order to obtain the above equations we have used the mathematical relations \({\text {d}}\zeta /{\text {d}}r=-\eta ^2{\text {d}}\zeta /{\text {d}}\eta \), and

$$\begin{aligned} \frac{{\text {d}}^2 \zeta }{{\text {d}}r ^2}=\eta ^4\frac{{\text {d}}^2\zeta }{{\text {d}}\eta ^2}+2\eta ^3\frac{{\text {d}}\zeta }{{\text {d}}\eta }, \end{aligned}$$
(54)

respectively. The system of Eqs.  (48), (49), (50) and (51) must be integrated with the initial conditions at infinity, given by \(M_{\mathrm{eff}}(0)=M_{\mathrm{eff}}^{(0)}\), \(\alpha (0)=0\), \(\zeta (0)=\zeta _0\), and \(u(0)=u_0\), respectively.

5.1 The Higgs potential: \(V(\zeta )=\mu ^2\zeta ^2+\nu \zeta ^4\)

5.1.1 General considerations

The Higgs-type potential

$$\begin{aligned} V(\zeta )=\mu ^2\zeta ^2+\nu \zeta ^4 \end{aligned}$$
(55)

plays a fundamental role in elementary particle physics. From a physical point of view we may assume that \(-\mu ^2 \) represents the mass of the three-form field associated to the gravitational interaction. For the strong interaction case the Higgs self-coupling constant \(\nu \) takes the value \(\nu \approx 1/8\) [38], a value which follows from the analysis of accelerator experiments. But of course in the case of the gravitational models in the presence of a three-form field the values of both \(\mu ^2\) and \(\xi \) may be very different from those suggested by elementary particle physics.

In the case of the Higgs-type potential of the three-form field the vacuum gravitational field equations take the form

$$\begin{aligned}&\frac{{\text {d}}\zeta }{{\text {d}}\eta }=u, \quad \frac{{\text {d}}M_{\mathrm{eff}}}{{\text {d}}\eta }=\frac{1}{2}\left[ \frac{F^2}{48}+\left( \mu ^2+3\nu \zeta ^2\right) \zeta ^2\right] \frac{1}{\eta ^4}, \end{aligned}$$
(56)
$$\begin{aligned}&\frac{{\text {d}}\alpha }{{\text {d}}\eta }=-\frac{F^2/48-\left( \mu ^2+\nu \zeta ^2\right) \zeta ^2+2\eta ^3M_{\mathrm{eff}}}{\eta ^3\left( 1-2\eta M_{\mathrm{eff}}\right) }, \end{aligned}$$
(57)
$$\begin{aligned}&\frac{{\text {d}}u}{{\text {d}}\eta }=-\frac{\left( 3\mu ^2+10\nu \zeta ^2\right) \zeta }{\eta ^{3}\left( 1-2\eta M_{\mathrm{eff}} \right) }u+\frac{1}{\eta ^{4}}G\left( \eta ,\zeta \right) \zeta -\frac{2\left( \mu ^2+2\nu \zeta ^2\right) }{\eta ^{4}},\nonumber \\ \end{aligned}$$
(58)

where \(G(\eta , \zeta )\) and \(F^2\) are given by

$$\begin{aligned} G\left( \eta ,\zeta \right)= & {} 2\eta ^{2}+\frac{2\left( \mu ^2+2\nu \zeta ^2\right) \zeta ^2}{2\left( 1-2\eta M_{\mathrm{eff}}\right) } \nonumber \\&\times \left\{ 2-\frac{\eta ^{2}+\left[ F^{2}/48+\left( \mu ^2+3\nu \zeta ^2\right) \zeta ^2\right] }{\eta ^{2}\left( 1-2\eta M_{\mathrm{eff}} \right) }\right\} \nonumber \\ \end{aligned}$$
(59)

and

$$\begin{aligned} F^{2}=-6\left\{ \left[ 4\eta -\frac{2\left( \mu ^2+2\nu \zeta ^2 \right) \zeta ^2}{\eta \left( 1-2\eta M_{\mathrm{eff}}\right) }\right] \zeta -2\eta ^2u\right\} ^{2}, \end{aligned}$$
(60)

respectively. Equations (56)–(58) must be considered with the initial conditions at infinity \(M_{\mathrm{eff}}=1\), \(\zeta (0)=10^{-5}\), \(u(0)=40\), and \(\alpha (0)=0\), respectively. In Figs. 1, 2, 3 and 4 we present the variations with respect to \(\eta \) of the metric tensor coefficients \(e^{\alpha }\), \(e^{-\beta }\), of the effective mass \(M_{\mathrm{eff}}\) and of the radial scalar function \(\zeta \).

Fig. 1
figure 1

Specific case of the Higgs potential: Variation of the metric tensor coefficient \(e^{\alpha }\) as a function of the coordinate \(\eta \), for \(\nu =0.01\) and for different values of \(\mu ^2\): \(\mu ^2=0.0001\) (solid curve), \(\mu ^2=0.00012\) (dotted curve), \(\mu ^2=0.00013\) (short dashed curve), \(\mu ^2=0.000138\) (dashed curve), and \(\mu ^2=0.000141\) (long dashed curve), respectively

Fig. 2
figure 2

Specific case of the Higgs potential: Variation of the metric tensor coefficient \(e^{-\beta }\) as a function of the coordinate \(\eta \), for \(\nu =0.01\), and for different values of \(\mu ^2\): \(\mu ^2=0.0001\) (solid curve), \(\mu ^2=0.00012\) (dotted curve), \(\mu ^2=0.00013\) (short dashed curve), \(\mu ^2=0.000138\) (dashed curve), and \(\mu ^2=0.000141\) (long dashed curve), respectively

Fig. 3
figure 3

Specific case of the Higgs potential: Variation of the effective mass \(M_{\mathrm{eff}}\) as a function of the coordinate \(\eta \), for \(\nu =0.01\), and for different values of \(\mu ^2\): \(\mu ^2=0.0001\) (solid curve), \(\mu ^2=0.00012\) (dotted curve), \(\mu ^2=0.00013\) (short dashed curve), \(\mu ^2=0.000138\) (dashed curve), and \(\mu ^2=0.000141\) (long dashed curve), respectively

Fig. 4
figure 4

Specific case of the Higgs potential: Variation of the radial scalar function \(\zeta \) as a function of the coordinate \(\eta \), for \(\nu =0.01\), and for different values of \(\mu ^2\): \(\mu ^2=0.0001\) (solid curve), \(\mu ^2=0.00012\) (dotted curve), \(\mu ^2=0.00013\) (short dashed curve), \(\mu ^2=0.000138\) (dashed curve), and \(\mu ^2=0.000141\) (long dashed curve), respectively

The variation of the metric tensor coefficient \(e^{\alpha }\) is represented in Fig. 1. The metric function monotonically decreases from its constant, Minkowskian value at infinity, to zero, a value reached for finite values of \(\eta \), and which defines the singular surface of the black hole, or its event horizon. The position of the event horizon is strongly dependent on the numerical values of \(\mu ^2\). A similar behavior characterizes the metric tensor component \(e^{-\beta }\), whose variation with respect to \(\eta \) is represented in Fig. 2. The inverse of the metric tensor component decreases linearly from infinity to the event horizon of the black hole, where the metric tensor becomes singular. Similarly to \(e^{\alpha }\), the variation of \(e^{-\beta }\) is significantly influenced by the numerical values of \(\mu ^2\). The changes in the effective mass \(M_{\mathrm{eff}}\) are plotted in Fig. 3. The mass increases rapidly from its initial value at infinity to a maximum value, reached before the event horizon, an effect due to the presence of the three-form field, and its mass–energy contribution to the mass of the central object. After reaching its maximum value the effective mass decreases before reaching the event horizon. However, for some particular values of \(\mu ^2\), the mass becomes approximately constant beginning for some finite value of \(\eta \). This indicates that the model enters very quickly in an approximate Schwarzschild regime, with \(e^{-\beta }\approx 1-2M_{\mathrm{eff}}\eta \), with \(M_{\mathrm{eff}}\) a function of the parameters of the Higgs-type potential, and of the initial conditions for \(\zeta \). This dependence on the initial conditions and on the parameters of the potential also determines the modifications of the position of the event horizon of the black hole, with respect to its standard general relativistic value. The dependence of the scalar radial function \(\zeta \) is depicted in Fig. 4. The behavior of \(\zeta \) indicates a complex dynamics, with \(\zeta \) increasing initially, reaching a maximum value, and then becoming again zero at the event horizon. For some values of the Higgs potential the variation has a quasi-oscillatory behavior, characterized by an alternation of local maxima and minima.

Fig. 5
figure 5

Variation of the metric tensor coefficient \(e^{\alpha }\), for the Higgs potential, as a function of the coordinate \(\eta \), for \(\mu ^2=0.000095\), and for different values of \(\nu \): \(\nu =0.04\) (solid curve), \(\nu =0.08\) (dotted curve), \(\nu =0.12\) (short dashed curve), \(\nu =0.16\) (dashed curve), and \(\nu (0)=0.20\) (long dashed curve), respectively

Fig. 6
figure 6

Variation of the metric tensor coefficient \(e^{-\beta }\), for the Higgs potential, as a function of the coordinate \(\eta \), for \(\mu ^2=0.000095\), and for different values of \(\nu \): \(\nu =0.04\) (solid curve), \(\nu =0.08\) (dotted curve), \(\nu =0.12\) (short dashed curve), \(\nu =0.16\) (dashed curve), and \(\nu (0)=0.20\) (long dashed curve), respectively

Fig. 7
figure 7

Variation of the effective mass \(M_{\mathrm{eff}}\), for the Higgs potential, as a function of the coordinate \(\eta \), for \(\mu ^2=0.000095\), and for different values of \(\nu \): \(\nu =0.04\) (solid curve), \(\nu =0.08\) (dotted curve), \(\nu =0.12\) (short dashed curve), \(\nu =0.16\) (dashed curve), and \(\nu (0)=0.20\) (long dashed curve), respectively

Fig. 8
figure 8

Variation of the radial scalar function \(\zeta \), for the Higgs potential, as a function of the coordinate \(\eta \), for \(\mu ^2=0.000095\), and for different values of \(\nu \): \(\nu =0.04\) (solid curve), \(\nu =0.08\) (dotted curve), \(\nu =0.12\) (short dashed curve), \(\nu =0.16\) (dashed curve), and \(\nu (0)=0.20\) (long dashed curve), respectively

The behavior of the geometric and physical quantities in vacuum for the three-form supported black holes also depend sensitively on the numerical values of the self-coupling constant \(\nu \) of the Higgs potential. In Figs. 5, 6, 7 and 8 we present the behavior of the metric tensor components, of the effective mass and of the scalar function \(\zeta \) for different values of \(\nu \). To integrate the gravitational field equations we have fixed the values of \(\mu ^2=0.000095\), \(\zeta (0)=0\), \(u(0)=40\), \(M_{\mathrm{eff}}(0)=1\), and \(\alpha (0)=0\), respectively, and we have varied the values of \(\nu \).

The metric tensor components \(e^{\alpha }\) and \(e^{-\beta }\), represented in Figs. 5 and 6, decrease from their Minkowski values at infinity to zero, corresponding to a finite value of \(\eta \), indicating the presence of a singularity corresponding to the formation of a black hole. Their numerical values depend effectively on the numerical values of \(\nu \), which also strongly influence the position of the event horizon. The effective mass, shown in Fig. 7, increases from its value at infinity towards a maximum value reached far away from the event horizon. The value of the effective mass also strongly depends on the self-coupling constant \(\nu \). For some numerical values of \(\nu \) the mass becomes roughly a constant beginning from a finite \(\eta \), and the geometry near the compact object becomes quasi-Schwarzschild, with the effective mass of the black hole strongly dependent on the values of \(\nu \). The radial scalar function \(\zeta \), represented in Fig. 8, also shows an effective dependence on \(\nu \), reaching, similarly to the effective mass, a finite value at the event horizon of the black hole.

The numerical values of the event horizon \(\eta _S\) are presented, for a selected value of the initial conditions and of the model parameters, in Table 1. In the adopted system of units the position of the Schwarzschild singularity corresponds to \(\eta _S=1/2\). The position of the event horizon \(r_s\) of the black hole is found to be \(r_s=r_g/\eta _S\), where \(r_g\) is the Schwarzschild gravitational radius of the object, defined as \(r_g=GM/c^2\), where M is the total mass of the object. For example, the physical position of the event horizon of the black hole supported by a three-form field with Higgs potential having \(\eta _S=0.88\) is located at \(r_s=1.13r_g\), indicating a black hole more extended that its Schwarzschild counterpart. For a three-form field black hole with \(\eta _S=0.44\), the event horizon is located at \(r_s=2.27r_g\). The three-dimensional distribution of the event horizons of the black holes supported by three-form fields is represented in Fig. 9.

Table 1 The position of the event horizon \(\eta _S\) for selected values of the initial conditions and Higgs model parameters
Fig. 9
figure 9

Distribution of the event horizons of the three-form field supported black holes as a function of the parameters \(\mu ^2\) and \(\nu \) of the Higgs potential for \(\zeta _0=\{1, 1.105, 1.21, 1.315, 1.42\}\times 10^{-4}\), and \(u_0=\{1, 2, 3, 4, 5 \} \times 10^{-2}\)

5.1.2 Interpolating functions

In order to facilitate the further investigations of the properties of the three-form field supported black holes in the following we will present some explicit analytic expressions for the basic physical and geometrical quantities, obtained from the interpolation of the numerical results. For the effective mass function \(M_{\mathrm{eff}}(\eta )\) we assume a general expression of the form

$$\begin{aligned} M_{\mathrm{eff}}(\eta )=A_\mathrm{{MH}}+B_\mathrm{{MH}}\eta +C_\mathrm{{MH}}\eta ^2, \end{aligned}$$
(61)

where the coefficients \(A_\mathrm{{MH}}\), \(B_\mathrm{{MH}}\), and \(C_\mathrm{{MH}}\) are functions of \((\mu ^2,\nu ,\zeta _0,u_0)\). In the following analysis we will concentrate mostly on the dependence on the initial conditions, and hence we will fix the parameters of the Higgs potential as \(\mu ^2=10^{-4}\), and \(\nu =10^{-2}\), respectively. Then we obtain

$$\begin{aligned} Av= & {} -1.915 - 3031.75 \zeta _0 - 509.327 \frac{\sqrt{\zeta _0}}{u_0}\nonumber \\&+142182 \frac{\zeta _0}{u_0} +0.1444 u_0-0.0016 u_0^2, \end{aligned}$$
(62)

with the correlation coefficient \(R^2=0.999\),

$$\begin{aligned} B_\mathrm{{MH}}= & {} 70.627+81785.6 \zeta _0 - 5.312 \times 10^7 \zeta _0^2 \nonumber \\&- 2.98 \times 10^6 \frac{\zeta _0}{u_0} - 3.27 u_0 + 0.036 u_0^2, \end{aligned}$$
(63)

with \(R^2=0.989\), and

$$\begin{aligned} C_\mathrm{{MH}}= & {} -294.508 + 916417 \zeta _0 - 1.4889 \times 10^7 \frac{\zeta _0}{u_0}\nonumber \\&+12.86 u_0 - 12456.8 \zeta _0 u_0 - 0.129447 u_0^2, \end{aligned}$$
(64)

with \(R^2=0.99\). The general expression of the metric tensor component \(e^\alpha \) is

$$\begin{aligned} e^{\alpha (\eta )}=A_{\alpha H}+B_{\alpha H}\eta +C_{\alpha H}\eta ^2, \end{aligned}$$
(65)

with the coefficients \(A_{\alpha H}\), \(B_{\alpha H}\), \(C_{\alpha H}\) given, for fixed \(\mu ^2\) and \(\nu \), as functions of \(\left( \zeta _0,u_0\right) \) by

$$\begin{aligned} A_{\alpha H}= & {} 0.847 -125.916 \zeta _0 - 36.981 \frac{\sqrt{\zeta _0}}{u_0}+9656.83 \frac{\zeta _0}{u_0}\nonumber \\&+0.0069 u_0-0.000059 u_0^2, \end{aligned}$$
(66)

with \(R^2=0.999\),

$$\begin{aligned} B_{\alpha H}= & {} 0.62+2036.19 \zeta _0 - 4.162 \times 10^6 \zeta _0^2 \nonumber \\&- 102912 \times 10^6 \frac{\zeta _0}{u_0} - 0.0942 u_0 + 0.000678 u_0^2,\nonumber \\ \end{aligned}$$
(67)

with \(R^2=0.999\), and

$$\begin{aligned} C_{\alpha H}= & {} -2.754 + 90221.2 \zeta _0 - 1.96 \times 10^6 \frac{\zeta _0}{u_0}\nonumber \\&+ 0.00058 u_0 - 882.929 \zeta _0 u_0 - 0.00193 u_0^2, \end{aligned}$$
(68)

with \(R^2=0.998\).

For the metric tensor component \(e^\beta \) the general interpolating function can be taken as

$$\begin{aligned} e^{\beta (\eta )}=A_{\beta H}+B_{\beta H}\eta +C_{\beta H}\eta ^2, \end{aligned}$$
(69)

with the coefficients \(A_{\beta H}\), \(B_{\beta H}\), and \(C_{\beta H}\) given as functions of \(\left( \zeta _0,u_0\right) \) by

$$\begin{aligned} A_{\beta H}= & {} 1.653 + 730.135 \zeta _0 + 67.886 \frac{\sqrt{\zeta _0}}{u_0}- 35382 \frac{\zeta _0}{u_0} \nonumber \\&+0.0324 u_0-0.000353 u_0^2, \end{aligned}$$
(70)

with \(R^2=0.999\),

$$\begin{aligned} B_{\beta H}= & {} 12.688 - 16943.1 \zeta _0 + 9.355 \times 10^6 \zeta _0^2 + 737250 \frac{\zeta _0}{u_0} \nonumber \\&- 1.12 u_0 + 0.0296 u_0^2 - 0.000243 u_0^3, \end{aligned}$$
(71)

with \(R^2=0.976\), and

$$\begin{aligned} C_{\beta H}= & {} 59.249 + 212369 \zeta _0 + 3.0549 \times 10^6 \frac{\zeta _0}{u_0} - 2.506 u_0 \nonumber \\&+ 2762.06 \zeta _0 u_0 + 0.0206 u_0^2, \end{aligned}$$
(72)

with \(R^2=0.999\).

5.2 The exponential potential: \(V(\zeta )=V_0e^{\lambda \zeta }\)

As a second example of black hole solutions supported by a three-form field with non-zero potential, we consider the case of the exponential-type potential, \(V (\zeta ) = V_0e^{\lambda \zeta }\), where \(V_0\) and \(\lambda \) are constants. There are many physical processes in string theory and elementary particle physics described by this type of potential. For example, an exponential-type potential is obtained in string-type theories and in four-dimensional effective Kaluza–Klein theories from the compactification of the higher dimensions. Moduli fields and non-perturbative effects in quantum field theory such as gaugino condensation can also generate exponential-type potentials for scalar fields [39]. The role of the exponential potential has been intensively investigated especially in the framework of scalar field cosmological and gravitational models and for many field configurations, including the inhomogeneous and homogeneous scalar fields [40,41,42,43,44,45,46,47,48]. In the presence of an exponential potential the gravitational field equations (48)–(51) take the form

$$\begin{aligned}&\frac{{\text {d}}\zeta }{{\text {d}}\eta }=u, \quad \frac{{\text {d}}M_{\mathrm{eff}}}{{\text {d}}\eta }=\frac{1}{2}\left[ \frac{F^2}{48}+V_0\left( \lambda \zeta -1\right) e^{\lambda \zeta }\right] \frac{1}{\eta ^4}, \end{aligned}$$
(73)
$$\begin{aligned}&\frac{{\text {d}}\alpha }{{\text {d}}\eta }=-\frac{F^2/48-V_0e^{\lambda \zeta }+2\eta ^3M_{\mathrm{eff}}}{\eta ^3\left( 1-2\eta M_{\mathrm{eff}}\right) }, \end{aligned}$$
(74)
$$\begin{aligned}&\frac{{\text {d}}u}{{\text {d}}\eta }=-\frac{V_0\lambda \left( 1+\lambda \zeta /2\right) e^{\lambda \zeta } }{\eta ^{3}\left( 1-2\eta M_{\mathrm{eff}} \right) }u+\frac{1}{\eta ^{4}}G\left( \eta ,\zeta \right) \zeta -\frac{\lambda V_0e^{\lambda \zeta }}{\eta ^{4}}, \nonumber \\ \end{aligned}$$
(75)

where

$$\begin{aligned} G\left( \eta ,\zeta \right)= & {} 2\eta ^{2}+\frac{\lambda V_0\zeta e^{\lambda \zeta }}{2\left( 1-2\eta M_{\mathrm{eff}}\right) }\nonumber \\&\times \left\{ 2-\frac{\eta ^{2}+\left[ F^{2}/48+V_0\left( \lambda \zeta -1\right) e^{\lambda \zeta }\right] }{\eta ^{2}\left( 1-2\eta M_{\mathrm{eff}} \right) }\right\} \nonumber \\ \end{aligned}$$
(76)

and

$$\begin{aligned} F^{2}=-6\left\{ \left[ 4\eta -\frac{\lambda V_0\zeta e^{\lambda \zeta }}{\eta \left[ 1-2\eta M_{\mathrm{eff}}\right] }\right] \zeta -2\eta ^2u\right\} ^{2}, \end{aligned}$$
(77)

respectively.

5.2.1 Naked singularity solutions

The description of the state and structure of ordinary material systems, forming an initial regular distribution, after the gravitational collapse, is one of the most important theoretical and observational problems in general relativity. There are two questions one should consider when investigating the gravitational collapse. The first question is to find the initial conditions of the gravitational collapse that lead to the formation of a black hole. On the other hand a careful investigation of the gravitational collapse shows that it does not end always with the creation of a black hole. Depending on the initial conditions, another type of object, called a naked singularity, can also be born as the final state of the collapse [49,50,51,52,53,54]. For reviews of the naked singularity problem see [55, 56], respectively.

The second question one must also necessarily consider is the question if the physically realistic collapse solutions of the Einstein gravitational field equations that indicate the formation of naked singularities do really correspond to existing natural objects, which can be observed by astrophysical or astronomical methods. If detected observationally, the existence of the naked singularities would be counterexamples of the Cosmic Censorship Hypothesis, proposed by Roger Penrose [57]. The Cosmic Censorship Hypothesis conjectures that curvature singularities are always covered in asymptotically flat spacetimes by event horizons. In fact, one can formulate the Cosmic Censorship Hypothesis in a strong sense (in a geometry that is physically appropriate naked singularities cannot form), and in a weak sense (if naked singularities do really exist, they are securely covered by an event horizon, and therefore they cannot be detected by far-away observers). There have been many attempts to prove the Cosmic Censorship Hypothesis (see [58] for a review of the early investigations and results in this field). For the possibilities of observationally identifying naked singularities see [59], and the references therein.

For a certain range of parameters and of initial conditions, naked singularity solutions of the gravitational field equations in the presence of a three-form field with exponential potential can also be obtained. In Figs. 10, 11, 12 and 13 we present the behavior of the metric tensor coefficients \(e^{\alpha }\), \(e^{-\beta }\), of the effective mass \(M_{\mathrm{eff}}\), and of the radial scalar function \(\zeta \) for the initial conditions \(M(0)=1\), \(\alpha (0)=0\), \(\zeta (0)=10^{-5}\), and \(u(0)=10^{-4}\), respectively. For \(\lambda \) we have adopted the value \(\lambda =-10^{-3}\), and we have slightly varied the numerical values of \(V_0\).

Fig. 10
figure 10

Variation of the metric tensor coefficient \(e^{\alpha }\) as a function of the coordinate \(\eta \) for the case of the exponential three-form potential \(V=V_0e^{\lambda \zeta }\), for \(\lambda =-10^{-3}\), and for different values of \(V_0 \): \(V_0=0.001\) (solid curve), \(V_0=0.0012\) (dotted curve), \(V_0=0.0013\) (short dashed curve), \(V_0=0.0014\) (dashed curve), and \(V_0=0.0015\) (long dashed curve), respectively

Fig. 11
figure 11

Variation of the metric tensor coefficient \(e^{-\beta }\) as a function of the coordinate \(\eta \) for the case of the exponential three-form potential \(V=V_0e^{\lambda \zeta }\), for \(\lambda =-10^{-3}\), and for different values of \(V_0 \): \(V_0=0.001\) (solid curve), \(V_0=0.0012\) (dotted curve), \(V_0=0.0013\) (short dashed curve), \(V_0=0.0014\) (dashed curve), and \(V_0=0.0015\) (long dashed curve), respectively

Fig. 12
figure 12

Variation of the effective mass \(M_{\mathrm{eff}}\) as a function of the coordinate \(\eta \) for the case of the exponential three-form potential \(V=V_0e^{\lambda \zeta }\), for \(\lambda =-10^{-3}\), and for different values of \(V_0 \): \(V_0=0.001\) (solid curve), \(V_0=0.0012\) (dotted curve), \(V_0=0.0013\) (short dashed curve), \(V_0=0.0014\) (dashed curve), and \(V_0=0.0015\) (long dashed curve), respectively

Fig. 13
figure 13

Variation of the radial scalar function \(\zeta \) as a function of the coordinate \(\eta \) for the case of the exponential three-form potential \(V=V_0e^{\lambda \zeta }\), for \(\lambda =-10^{-3}\), and for different values of \(V_0 \): \(V_0=0.001\) (solid curve), \(V_0=0.0012\) (dotted curve), \(V_0=0.0013\) (short dashed curve), \(V_0=0.0014\) (dashed curve), and \(V_0=0.0015\) (long dashed curve), respectively

As one can see from Figs. 10 and 11, the metric tensor coefficients are monotonically increasing functions of \(\eta \), and they are singular only at the origin \(\eta \rightarrow \infty \), or, equivalently, \(r\rightarrow 0\). The correspondent massive object does not have an event horizon, and therefore it corresponds to a naked singularity, with the only singular point located at the center. The effective mass of the naked singularity, presented in Fig. 12, becomes negative at infinity, and takes a constant, negative value up to the singular center of the naked singularity. Hence the metric of this exotic object can be represented as

$$\begin{aligned} e^{\alpha (r)}=e^{-\beta (r)}=1+\frac{2M_{\mathrm{eff}}\left( V_0,\lambda , \zeta _0,\zeta '(0)\right) }{r}. \end{aligned}$$
(78)

The numerical values of the (negative) effective mass are determined by the initial conditions of the gravitational field equations at infinity, as well as by the parameters of the exponential potential. The solutions of the gravitational field equations depend sensitively on these parameters. The radial scalar function \(\zeta \), shown in Fig. 13, diverges at the center of the naked singularity. Its behavior is also dependent on the initial conditions used to solve the gravitational field equations, and on the parameters of the exponential potential.

5.2.2 Black hole solutions

The static spherically symmetric vacuum gravitational field equations in the presence of a three-form field also admit black hole-type solutions. In Figs. 14, 15, 16 and 17 we present the results of the numerical integration of the gravitational field equations for \(M(0)=1\), \(\alpha (0)=0\), \(\zeta (0)=10^{-2}\), \(\zeta '(0)=10^{-1}\), \(V_0=9.9\times 10^{-10}\), and different values of \(\lambda \).

Fig. 14
figure 14

Variation of the metric tensor coefficient \(e^{\alpha }\) as a function of the coordinate \(\eta \) for the case of the exponential three-form potential \(V=V_0e^{\lambda \zeta }\), for \(V_0 =9.9\times 10^{-10}\), and for different values of \(\lambda \): \(\lambda =-40\) (solid curve), \(\lambda =-120\) (dotted curve), \(\lambda =-200\) (short dashed curve), \(\lambda =-280\) (dashed curve), and \(\lambda =-360\) (long dashed curve), respectively

Fig. 15
figure 15

Variation of the metric tensor coefficient \(e^{-\beta }\) as a function of the coordinate \(\eta \) for the case of the exponential three-form potential \(V=V_0e^{\lambda \zeta }\), for \(V_0 =9.9\times 10^{-10}\), and for different values of \(\lambda \): \(\lambda =-40\) (solid curve), \(\lambda =-120\) (dotted curve), \(\lambda =-200\) (short dashed curve), \(\lambda =-280\) (dashed curve), and \(\lambda =-360\) (long dashed curve), respectively

Fig. 16
figure 16

Variation of the effective mass \(M_{\mathrm{eff}}\) as a function of the coordinate \(\eta \) for the case of the exponential three-form potential \(V=V_0e^{\lambda \zeta }\), for \(V_0 =9.9\times 10^{-10}\), and for different values of \(\lambda \): \(\lambda =-40\) (solid curve), \(\lambda =-120\) (dotted curve), \(\lambda =-200\) (short dashed curve), \(\lambda =-280\) (dashed curve), and \(\lambda =-360\) (long dashed curve), respectively

Fig. 17
figure 17

Variation of the radial scalar function \(\zeta \) as a function of the coordinate \(\eta \) for the case of the exponential three-form potential \(V=V_0e^{\lambda \zeta }\), for \(V_0 =9.9\times 10^{-10}\), and for different values of \(\lambda \): \(\lambda =-40\) (solid curve), \(\lambda =-120\) (dotted curve), \(\lambda =-200\) (short dashed curve), \(\lambda =-280\) (dashed curve), and \(\lambda =-360\) (long dashed curve), respectively

As one can see from Figs. 14 and 15, the metric tensor coefficients \(e^{\alpha }\) and \(e^{-\beta }\) decrease monotonically from their Minkowskian values at infinity to zero, a value reached for a finite value of \(\eta =\eta _S\). Hence the compact object possesses an event horizon, and is thus a black hole. The effective mass \(M_{\mathrm{eff}}\), depicted in Fig. 16, decreases very quickly from its initial value at infinity, and becomes a constant, having the same numerical value from infinity to the event horizon of the black hole. Hence the metric is of the Schwarzschild type, with

$$\begin{aligned} e^{\alpha (r)}=e^{-\beta (r)}=1-\frac{2M_{\mathrm{eff}}\left( V_0,\lambda , \zeta (0),\zeta '(0)\right) }{r}, \end{aligned}$$
(79)

with the effective mass, and the position of the event horizon depending on the initial conditions at infinity, and on the parameters of the exponential potential. The radial scalar function increases rapidly from infinity when approaching the event horizon, and takes a finite value for \(\eta =\eta _S\).

We can obtain an interpolating expression for the effective mass function \(M_{\mathrm{eff}}\) in the case of the exponential potential,

$$\begin{aligned} M_{\mathrm{eff}}(\eta )\approx A_\mathrm{{MH}}^{(\mathrm exp)}, \end{aligned}$$
(80)

with the coefficient \(A_\mathrm{{MH}}^{(\mathrm exp)}\) given, for fixed \(V_0\) and \(\lambda \), by

$$\begin{aligned} A_\mathrm{{MH}}^{(\mathrm exp)}\approx & {} 0.74 + 9.78 \zeta _0 + 15.812 \frac{\sqrt{\zeta _0}}{u_0} - 146.854 \frac{\zeta _0}{u_0} \nonumber \\&+ 0.00134 u_0 - 0.0000528 u_0^2, \end{aligned}$$
(81)

with \(R^2=0.997\).

Fig. 18
figure 18

Variation of the position of the event horizons of the three-form field supported black holes as a function of the parameters \(V_0\) and \(\lambda \) of the exponential potential for \(\zeta _0=\{0.0004, 0.0064, 0.0124, 0.0184, 0.0244\}\) and \(u_0=\{4,8,16,32,64\}\)

The distribution of the position of the black holes event horizon supported by a three-form field with exponential potential are represented, for fixed \(V_0\) and \(\lambda \), in Fig. 18.

6 Thermodynamic properties of black holes

In the present section we consider the thermodynamic properties of the black hole solutions supported by a three-form field. In particular we will concentrate on the surface gravity of the black holes, their Hawking temperature, their specific heat, their entropy and their Hawking luminosity. In our investigation of the vacuum field equations in the three-form fields model we have adopted the simplifying assumption that the effective mass function and the lapse function \(e^{\alpha }\) are functions of the radial coordinate r only. Hence the spacetime is static and a timelike Killing vector \(t^{\mu }\) exists [60, 61]. In the present section, for the sake of clarity, we will restore the physical units in all mathematical expressions.

6.1 Brief summary of black hole thermodynamics

For a static black hole that possesses a Killing horizon the definition of the surface gravity \(\tilde{\kappa }\) is given by [60, 61]

$$\begin{aligned} t^{\mu }\nabla _{\mu }t^{\nu }=t^{\nu }\tilde{\kappa }, \end{aligned}$$
(82)

where \(t^{\mu }\) is a Killing vector, and \(\nabla _{\mu }\) denotes the covariant derivative with respect to the metric. For a static, spherically symmetric geometry, with the line element given by

$$\begin{aligned} ds^2=-\tilde{\sigma } ^2 (r)f(r)c^2dt^2+\frac{dr^2}{f(r)}+r^2d\Omega ^2, \end{aligned}$$
(83)

we can adopt a suitable normalized Killing vector defined as \(t^{\mu }=\left( 1/\tilde{\sigma }_{\infty },0,0,0\right) \). Then the surface gravity of the black hole is found to be [61]

$$\begin{aligned} \tilde{\kappa }=\left( \frac{\tilde{\sigma } _\mathrm{hor}}{\tilde{\sigma } _{\infty }}\right) \frac{c^4}{4GM_\mathrm{hor}}\left. \left[ 1-\frac{2GM'(r)}{c^2}\right] \right| _\mathrm{hor}. \end{aligned}$$
(84)

The subscript hor requires that all physical quantities must be evaluated on the outer apparent horizon. If the function \(\tilde{\sigma }\equiv 1\), and \(M=\mathrm{constant}\), from the above definition we reobtain the standard result of the surface gravity of a Schwarzschild black hole, which is given by [60]

$$\begin{aligned} \tilde{\kappa }=\frac{c^4}{4GM_\mathrm{hor}}. \end{aligned}$$
(85)

The temperature \(T_\mathrm{{BH}}\) of the black hole is found to be

$$\begin{aligned} T_\mathrm{{BH}}=\frac{\hbar }{2\pi ck_\mathrm{{B}}} \tilde{\kappa }, \end{aligned}$$
(86)

where \(k_\mathrm{{B}}\) is Boltzmann’s constant. Equivalently, in the variable \(r=r_g/\eta \) we obtain for the Hawking temperature of the black hole the expression

$$\begin{aligned} T_\mathrm{{BH}}= & {} \frac{T_H}{M_{\mathrm{eff}}\left( \eta _S\right) }\left. \left( 1+\eta ^2\frac{{\text {d}}M_{\mathrm{eff}}\left( \eta \right) }{{\text {d}}\eta }\right) \right| _{\eta =\eta _S} \nonumber \\= & {} T_H\left. \theta (\eta )\right| _{\eta =\eta _S}, \end{aligned}$$
(87)

where

$$\begin{aligned} T_H=\frac{\hbar c^3}{8\pi Gk_\mathrm{{B}}M_{0}}, \end{aligned}$$
(88)

\(M_0\) is the standard general relativistic mass of the black hole, and we have denoted

$$\begin{aligned} \theta (\eta )=\frac{1}{M_{\mathrm{eff}}\left( \eta \right) }\left( 1+\eta ^2\frac{{\text {d}}M_{\mathrm{eff}}\left( \eta \right) }{{\text {d}}\eta }\right) . \end{aligned}$$
(89)

The specific heat \(C_\mathrm{{BH}}\) of the black hole is defined as

$$\begin{aligned} C_\mathrm{{BH}}= & {} \frac{{\text {d}}M}{{\text {d}}T_\mathrm{{BH}}}=\left. \frac{{\text {d}}M}{{\text {d}}r}\frac{{\text {d}}r}{{\text {d}}T_\mathrm{{BH}}}\right| _{r=r_\mathrm{hor}}, \end{aligned}$$
(90)

and it takes the dimensionless form

$$\begin{aligned} C_\mathrm{{BH}}=\frac{M_0}{T_H} \left. \frac{{\text {d}}M_{\mathrm{eff}}\left( \eta \right) }{{\text {d}}\eta }\frac{{\text {d}}\eta }{{\text {d}}\theta }\right| _{\eta =\eta _S}. \end{aligned}$$
(91)

The Hawking entropy \(S_\mathrm{{BH}}\) of the black hole is found to be

$$\begin{aligned} S_\mathrm{{BH}}= & {} \int _{r_{in}}^{r_\mathrm{hor}}{\frac{{\text {d}}M}{T_\mathrm{{BH}}}}=\int _{r_{in}}^{r_\mathrm{hor}}{\frac{1}{T_\mathrm{{BH}}}\frac{{\text {d}}M}{{\text {d}}r}{\text {d}}r}, \end{aligned}$$
(92)

or, in a dimensionless form,

$$\begin{aligned} S_\mathrm{{BH}}\left( \eta _S\right) =C_H\int _0^{\eta _S}{\frac{1}{\theta \left( \eta \right) }\frac{{\text {d}}M_{\mathrm{eff}}\left( \eta \right) }{{\text {d}}\eta }{\text {d}}\eta }. \end{aligned}$$
(93)

The black hole luminosity due to Hawking evaporation can be found to be

$$\begin{aligned} L_\mathrm{{BH}}=-\frac{{\text {d}}M}{{\text {d}}t}=-\sigma A_\mathrm{{BH}}T_\mathrm{{BH}}^4, \end{aligned}$$
(94)

where \(\sigma \) is a parameter depending on the adopted physical model, while

$$\begin{aligned} A_\mathrm{{BH}}=4\pi r_\mathrm{hor}^2 \end{aligned}$$
(95)

is the surface area of the event horizon. Then for the black hole evaporation time \(\tau \) we find

$$\begin{aligned} \tau= & {} \int _{t_{in}}^{t_{fin}}{{\text {d}}t}=-\frac{1}{4\pi \sigma }\int _{t_{in}}^{t_{fin}}{\frac{{\text {d}}M}{r_\mathrm{hor}^2T_\mathrm{{BH}}^4}}. \end{aligned}$$
(96)

Hence for the black hole evaporation time \(\tau \) we obtain the expression

$$\begin{aligned} \tau= & {} \int _{t_{in}}^{t_{fin}}{{\text {d}}t}=-\frac{1}{4\pi \sigma }\int _{t_{in}}^{t_{fin}}{\frac{{\text {d}}M}{r_\mathrm{hor}^2T_\mathrm{{BH}}^4}}, \end{aligned}$$
(97)

or, in an equivalent dimensionless form,

$$\begin{aligned} \tau \left( \eta _S\right) =-\tau _H\int _0^{\eta _S}{\frac{1}{\eta ^2\theta ^4\left( \eta \right) }\frac{{\text {d}}M_{\mathrm{eff}}\left( \eta \right) }{{\text {d}}\eta }{\text {d}}\eta }, \end{aligned}$$
(98)

where we have denoted

$$\begin{aligned} \tau _H=\frac{c^4}{8\pi G^2\sigma M_{0}T_\mathrm{{BH}}^4}. \end{aligned}$$
(99)

6.2 Thermodynamics of the Higgs-type black holes

With the help of the interpolating function for the effective mass the Hawking temperature of a black hole supported by a three-form field is obtained explicitly:

$$\begin{aligned} T_\mathrm{{BH}}\left( \eta _S\right)= & {} T_H\left. \frac{1+\eta ^2\left( B_\mathrm{{MH}}+2C_\mathrm{{MH}}\eta \right) }{A_\mathrm{{MH}}+B_\mathrm{{MH}}\eta +C_\mathrm{{MH}}\eta ^2}\right| _{\eta =\eta _S}. \end{aligned}$$
(100)

The variations of the Hawking temperature of the three-form field black holes are represented in Fig. 19.

Fig. 19
figure 19

Variation of the dimensionless Hawking temperature of the three-form field supported black holes as a function of the event horizon radius \(\eta _S\) for the case of the Higgs-type three-form field potential, for \(\mu ^2=10^{-4}\), \(\nu =10^{-2}\), \(u_0 \in [40,60]\), and for different values of \(\zeta _0 \): \(\zeta _0=10^{-5}\) (solid curve), \(\zeta _0=2\times 10^{-5}\) (dotted curve), \(\zeta _0 =4\times 10^{-5}\) (short dashed curve), \(\zeta _0=8\times 10^{-5}\) (dashed curve), and \(\zeta _0 =16\times 10^{-5}\) (long dashed curve), respectively

The Hawking temperature depends on the position of the event horizon of the three-form black hole, and it monotonically increases in time with \(\eta _S\). For \(\eta _S=0.80\), \(T_\mathrm{{BH}}\approx 1.6T_H\), while for \(\eta _S=0.60\), \(T_\mathrm{{BH}}\approx 1.2T_H\). The standard Hawking temperature is of the order of \(T_H=6.169\times 10^{-8}\times \left( M_{\odot }/M_0\right) \), and for astrophysical-type three-form field black holes, having large masses, the shifts in the position of the event horizon produce negligible effects.

The specific heat of the Higgs-type three-form field black hole can be found to be

$$\begin{aligned} C_\mathrm{{BH}}\left( \eta _S\right)= & {} \Big \{C_H\left( B_\mathrm{{MH}}+2 C_\mathrm{{MH}} \eta \right) \times \nonumber \\&\left[ A_\mathrm{{MH}} +\eta \left( B_\mathrm{{MH}}+C_\mathrm{{MH}} \eta \right) \right] ^2 \Big \}\Big / \nonumber \\&\qquad \Big \{B_\mathrm{{MH}} \left( 2 A_\mathrm{{MH}} \eta +4 C_\mathrm{{MH}} \eta ^3 -1\right) \nonumber \\&\qquad +\eta \left[ C_\mathrm{{MH}} \left( 6 A_\mathrm{{MH}} \eta -2\right) +2 C_\mathrm{{MH}}^2 \eta ^3\right] \nonumber \\&\qquad +B_\mathrm{{MH}}^2 \eta ^2 \Big \}, \end{aligned}$$
(101)

and its variation with respect to the event horizon is represented in Fig. 20.

Fig. 20
figure 20

Variation of the dimensionless specific heat of the three-form field supported black holes as a function of the event horizon radius \(\eta _S\) for the case of the Higgs-type three-form field potential, for \(\mu ^2=10^{-4}\), \(\nu =10^{-2}\), \(u_0 \in [40,60]\), and for different values of \(\zeta _0 \): \(\zeta _0=10^{-5}\) (solid curve), \(\zeta _0=2\times 10^{-5}\) (dotted curve), \(\zeta _0 =4\times 10^{-5}\) (short dashed curve), \(\zeta _0=8\times 10^{-5}\) (dashed curve), and \(\zeta _0 =16\times 10^{-5}\) (long dashed curve), respectively

The specific heat of the three-form black holes shows a complicated behavior. After initially increasing as a function of \(\eta _S\), \(C_\mathrm{{BH}}\) reaches a maximum, and then it monotonically decreases towards a minimum value reached at \(\eta )S\approx 0.90\). In the range \(\eta _S\in (0.50,0.65)\) there is an increase in the numerical values of \(C_\mathrm{{BH}}\) as compared to the standard general relativistic case, so that, for \(\eta _S=0.65\), \(C_\mathrm{{BH}}\approx 1.45C_H\).

The variation of the black hole entropy as a function of the event horizon, as given by Eq. (93), is represented in Fig. 21. The behavior of the Hawking entropy of the three-form field black holes has a similar behavior like their specific heat. The entropies are monotonically increasing functions for small \(\eta _S\), they reach a maximum, and they decrease for larger values of \(\eta _S\). The maximum values of the entropy are of the order \(S_\mathrm{{BH}}\approx 1.5S_H\).

Fig. 21
figure 21

Variation of the dimensionless Hawking entropy of the three-form field supported black holes as a function of the event horizon radius \(\eta _S\) for the case of the Higgs-type three-form field potential, for \(\mu ^2=10^{-4}\), \(\nu =10^{-2}\), \(u_0 \in [40,60]\), and for different values of \(\zeta _0 \): \(\zeta _0=10^{-5}\) (solid curve), \(\zeta _0=2\times 10^{-5}\) (dotted curve), \(\zeta _0 =4\times 10^{-5}\) (short dashed curve), \(\zeta _0=8\times 10^{-5}\) (dashed curve), and \(\zeta _0 =16\times 10^{-5}\) (long dashed curve), respectively

Fig. 22
figure 22

Variation of the dimensionless evaporation time of the three-form field supported black holes as a function of the event horizon radius \(\eta _S\) for the case of the Higgs-type three-form field potential, for \(\mu ^2=10^{-4}\), \(\nu =10^{-2}\), \(u_0 \exists [40,60]\), and for different values of \(\zeta _0 \): \(\zeta _0=10^{-5}\) (solid curve), \(\zeta _0=2\times 10^{-5}\) (dotted curve), \(\zeta _0 =4\times 10^{-5}\) (short dashed curve), \(\zeta _0=8\times 10^{-5}\) (dashed curve), and \(\zeta _0 =16\times 10^{-5}\) (long dashed curve), respectively

With the use of the general Eq. (98), the ratio of the black hole evaporation time and of the Hawking evaporation time is represented in Fig. 22. As a function of \(\eta _S\) the evaporation time monotonically decreases from a maximum value \(\tau _\mathrm{{BH}}\approx 120 \tau _H\), reached for \(\eta _S\approx 0.55\), to a value of \(\tau _\mathrm{{BH}}\approx 20 \tau _H\) for \(\eta _S\approx 0.85\). Since the standard Hawking evaporation time of a black hole is of the order \(\tau _H\approx 4.8\times 10^{-27}\times \left( M_0/\mathrm{g}\right) ^3\), it turns out that the evaporation time for three-form field supported black holes can be one or two orders of magnitude higher. Even so, the evaporation time for astrophysical size objects remains very high, with a three-form field black hole having one solar mass completely evaporating via Hawking radiation in around \(10^{62}-10^{63}\) years.

6.3 Exponential potential-type black holes

In the presence of an exponential potential of the three-form field, the metric of the black holes are quasi-Schwarzschild, with the effective mass a constant for most of the range of the variation of \(\eta \). Hence the thermodynamical properties of the black holes can be obtained by the thermodynamic properties of the Schwarzschild black holes, with the mass substituted by the effective mass \(M_{\mathrm{eff}}\), as given by Eq. (80). For the Hawking temperature of the exponential-type black hole we obtain

$$\begin{aligned} T_\mathrm{{BH}}=\frac{\hbar c^3}{8\pi Gk_\mathrm{{B}}M_{\mathrm{eff}}\left( \eta _S\right) }\approx \frac{\hbar c^3}{8\pi Gk_\mathrm{{B}}A_\mathrm{{MH}}^{(\mathrm exp)}\left( V_0,\lambda , \zeta _0,u_0\right) }. \end{aligned}$$
(102)

Hence the Hawking temperature of the black hole is dependent on the parameters of the potential, and of the initial conditions at infinity of the field \(\zeta \). For the specific heat of the black hole we obtain

$$\begin{aligned} C_\mathrm{{BH}}=-\frac{8\pi k_\mathrm{{B}}G}{\hbar G}M_{\mathrm{eff}}\approx -\frac{8\pi k_\mathrm{{B}}G}{\hbar G}A_\mathrm{{MH}}^{(\mathrm exp)}\left( V_0,\lambda , \zeta _0,u_0\right) . \end{aligned}$$
(103)

There is a dependence of the specific heat on the potential parameters, and on the initial values for \(\zeta \). The negative sign indicates that, as a black hole loses mass, and hence energy, its temperature increases. For the entropy of the black hole we can write down the standard expression,

$$\begin{aligned} S_\mathrm{{BH}}= & {} \frac{k_\mathrm{{B}} c^3}{4\hbar G}A_\mathrm{{BH}}=\frac{\pi k_\mathrm{{B}}c^3}{\hbar G}r_\mathrm{hor}^2 \nonumber \\= & {} \frac{\pi k_\mathrm{{B}} G M_0^2}{\hbar c}\frac{1}{\eta _S^2\left( V_0,\lambda , \zeta _0,u_0\right) }. \end{aligned}$$
(104)

The numerical value of the black hole entropy is determined by the mass of the central compact object, as well as of the initial condition at infinity of the radial scalar function \(\zeta \). Finally, for the rate of the mass loss we have the relation

$$\begin{aligned} \frac{{\text {d}}M}{{\text {d}}t}=\frac{\hbar c^4}{15360 \pi G^2}\frac{1}{M^2}, \end{aligned}$$
(105)

which yields for the lifetime of the black hole the expression

$$\begin{aligned} t= & {} \frac{5120\pi G^2M_0^3}{\hbar c^4}M_{\mathrm{eff}}^3=\frac{5120\pi G^2 M_0^3}{\hbar c^4}M_{\mathrm{eff}}^3 \nonumber \\= & {} \frac{5120\pi G^2M_0^3}{\hbar c^4}\left[ A_\mathrm{{MH}}^{(\mathrm exp)}\left( V_0,\lambda , \zeta _0,u_0\right) \right] ^3. \end{aligned}$$
(106)

The corrections to the black hole lifetime are given by the third power of the function \(A_\mathrm{{MH}}^{(\mathrm exp)}\left( V_0,\lambda , \zeta _0,u_0\right) \), determined by the parameters of the potential and the initial conditions at infinity. However, the evaporation time of the black holes is not significantly influenced by the potential parameters, and the initial conditions.

7 Discussion and final remarks

In the present paper we have investigated the possible existence of massive compact astrophysical objects, described by black hole and naked singularity-type geometries, in the framework of the three-form field gravitational theory, in which the standard Hilbert–Einstein action of general relativity is extended by the addition of the Lagrangian of a three-form field. In order to investigate the gravitational properties of the model we have considered the simplest case, corresponding to a vacuum static and spherically symmetric geometry. In this case the system of gravitational field equations depend on the scalar radial function \(\zeta \), the radial component of the dual vector \(B^{\delta }\) of \(A_{\alpha \beta \gamma }\), and on the arbitrary potential \(V(\zeta )\) of the three-form field.

Even within the simple vacuum spherically symmetric static model the field equations of the theory become extremely complicated. However, the exact solution of the field equations can be obtained in the case of a constant potential. In this case the metric functions can be obtained in a form similar to the Schwarzschild–de Sitter geometry, but with the solution containing two arbitrary integration constants \(c_1\) and \(c_2\). Similarly to the metric, the radial scalar function also depends on two arbitrary constants. Depending on the choice of these constants, several types of compact geometries can be obtained. We can first reproduce the standard Schwarzschild–de Sitter geometry, in which the three-form field generates a cosmological constant. However, different choices of the constants are also possible, with the solutions corresponding to the Schwarzschild–anti de Sitter geometry, or, more interestingly, to naked singularities having no event horizon, and with the singularity located at the center \(r=0\) of the massive object.

For arbitrary non-constant potentials \(V(\zeta )\) in order to obtain solutions of the vacuum field equations one must use numerical methods. In order to investigate numerically the static spherically symmetric Einstein field equations in the presence of a three-form field we have reformulated them in a dimensionless form, and, moreover, we have introduced as the independent variable \(\eta \) the inverse of the radial coordinate \(\eta =r_g/r\). In this representation the numerical integration procedure of the field equations is significantly simplified. On the other hand to proceed with the numerical integration we need to fix the numerical values of the effective mass, of the radial scalar field \(\zeta \), and of its derivative \(\zeta '\) at infinity. In our present numerical approach we have assumed that at infinity the geometry is asymptotically flat, and therefore the metric tensor components take their Minkowskian values when \(r\rightarrow \infty \). But at infinity the values of \(\zeta \), and of \(\zeta '\) may be arbitrary (but small), and this choice is consistent with the interpretation of the radial scalar field as a cosmological constant.

When integrating numerically the gravitational field equations two possible behaviors are detected. The presence of a singular behavior at finite \(\eta =\eta _S\) in the field equations, or, more precisely, in the variation with the distance of the metric tensor coefficients, is explained as indicating the existence of an event horizon. We detect the presence of the horizon from the conditions \(e^{\alpha \left( \eta _S\right) }=0\) and \(e^{-\beta \left( \eta _S\right) }=0\), respectively. Consequently, a singularity at finite \(\eta \) corresponds to a black hole-type massive astrophysical object. The physical mass of the black hole corresponds to the effective mass of the three-form field model, which is found to be the sum of the standard mass of the black hole plus the energy contribution from the scalar component of the three-form field. The second situation is related to the location of the singularity at the center of the massive object \(r=0\) only, with the corresponding object representing a naked singularity, an object not covered by an event horizon.

We have considered two classes of numerical solutions of the gravitational field equations in the three-form field theory, corresponding to two choices of the three-form field potential \(V (\zeta )\). The potentials we have chosen for our investigations are the Higgs-type potential, and the exponential-type potential, respectively. In the case of the Higgs potential we have restricted our investigations to initial conditions that lead to the formation of an event horizon, and consequently of black holes. The locations of the event horizons are strongly dependent on the values at infinity of the radial scalar field and of its derivative (the initial conditions), indicating the existence of an important and subtle relation between the initial values of the three-form field at infinity, and the black hole properties. The location of the event horizon is also strongly dependent on the parameters \(\mu ^2\) and \(\nu \) of the potential Higgs-type potential, indicating a complex multi-parametric dependence of the position of the event horizon, and of the black hole properties. The event horizon for three-form field models with Higgs potential can be located at distances of the order of \(2r_g\) from the massive object center, at distance much higher than those of the standard Schwarzschild black holes, indicating the formation of more massive black holes than predicted by standard general relativity. In the case of the Higgs potential the numerical results can be fitted well by some simple analytical that depend on the initial conditions at infinity.

In the case of the exponential potential we have identified numerically two distinct classes of solutions, depending on the initial values of \(\zeta \) and \(\zeta '\) at infinity, and on the potential parameters. These two classes correspond to naked singularity-type solutions, with the singularity located at the center, and no event horizon present, and to standard black holes, characterized by the presence of an event horizon. In both cases the numerical solutions for the metric functions can be approximated well by a Schwarzschild-type form with \(e^{\alpha }=e^{-\beta }=1\pm 2M_{\mathrm{eff}}/r\), where \(M_{\mathrm{eff}}\) depends on the initial conditions at infinity, and on the parameters of the exponential potential. The positive sign indicates the presence of a naked singularity, while the negative sign corresponds to an object with an event horizon. The analytical fittings of the numerical results are extremely useful in the study of the thermodynamic properties of the three-form field black holes. They also greatly simplify the study of the dynamics and motion of matter particles around black holes and naked singularities.

One important question in the physics of the naked singularities is if the energy conditions are satisfied for this type of objects. For example, in [53] it was shown that in the case of a mixture of a null charged strange quark fluid and radiation collapsing in a Vaidya spacetime a naked singularity can be formed, with the matter satisfying the weak, strong and dominant energy conditions. The formation of a naked singularity or of a black hole depends on the initial distribution of the density and velocity, and on the constitutive nature of the collapsing matter. In the case of the naked singularities generated by the three form fields one could also formulate a number of energy conditions, similar to the case of ordinary matter, as \(\rho >0\), \(p_r>0\), \(p>0\), \(\rho>p_r>0\), and \(\rho>p>0\). However, from a physical point of view, these conditions refer to some effective quantities constructed from the three form field, and generally they are not satisfied for arbitrary three form field configurations.

The problem of the existence in nature of the naked singularities is closely related to their stability properties. The Reissner–Nordström naked singularity with \(|Q|>M\) is unstable under small linear perturbations, and a similar instability occurs for the rotating (Kerr) solution with angular momentum \(a>M\) [62]. These results suggest that these spacetimes cannot be the endpoint of physical gravitational collapse. In [63] it was shown that as a result of the gravitational collapse of a spherically symmetric scalar field the formation of the naked singularities occurs, but this phenomenon is unstable. However, as compared to the above-mentioned types of naked singularities, those formed in the presence of a three form field are of a modified Schwarzschild type, with the sign of the term M/r changed from negative to positive. The problem of their stability/instability is a very interesting one, due to the different nature of the metric. It is well known that the Schwarzschild black hole is linearly stable under gravitational and electromagnetic perturbations [64], so that one may conjecture that a similar property holds for a three form field supported by Schwarzschild-type naked singularities. Hence the investigation of the stability properties of the three form field naked singularities may lead to different results as compared to the Reissner–Nordström, Kerr or scalar field collapse cases.

We have also investigated in detail the thermodynamic properties of the three-form field black hole solutions obtained by numerical methods. The Hawking temperature is an interesting and essential physical property of black holes. The horizon temperature of the three-form black holes indicates a strong dependence on the initial conditions at infinity of the radial scalar field, and of the properties of the radial scalar field potential. This is very different from the properties of the standard general relativistic Hawking temperature, which depends only on the mass of the black hole, and is independent of the asymptotic conditions at infinity. Similar properties characterize the behavior of the specific heat, entropy and evaporation time of the three-form field black holes. The numerical results show that these quantities are also strongly dependent on the initial conditions of the radial scalar field \(\zeta \) at infinity.

Moreover, in the three-form field gravitational theory the black hole evaporation times may be very different as compared to the similar results in standard general relativity. But we should note that our results on the thermodynamics of black holes, obtained for only two radial scalar field potentials and for a limited range of initial conditions at infinity may be considered of qualitative nature only. But even at this qualitative level they show the complexity of the compact objects supported by the three-form fields, and of the interesting physical and astrophysical processes related to them. In particular, the analytic representations of the numerical results may be applied for the study of the electromagnetic properties of the thin accretion disks existing around black holes. These properties may help in discriminating the three-form field black holes from other similar theoretical objects as well as from their general relativistic counterparts, and for allowing to obtain observational constraints on the model parameters.

The no-hair theorem [65,66,67,68] is an important result in black hole physics. It states that asymptotically flat black holes do not allow for the presence of external nontrivial scalar fields, with non-negative field potential \(V (\phi )\). The results obtained in our present investigations indicate that in its standard formulation the no-hair theorem cannot be extended to the static spherically symmetric solutions of the three-form field gravitational theory. All the numerical black hole solutions we have obtained are asymptotically flat, and three-form fields with positive radial scalar field potentials exist around them. Similar results have been obtained for the black hole solutions in the hybrid metric-Palatini gravity theory [24], suggesting that the no-hair theorems may not be valid in some modified gravity theories. But the answer to the question if such properties are a result of the particular choice of the three-form field radial scalar potentials, of the initial conditions and of the scalar potential parameters, or that they are some generic properties of the theory, requires further and detailed investigations, at both theoretical and computational levels.

Different types of physical fields may play an important role in astrophysics and cosmology, especially as potential constituents of the dark energy and dark matter components of the Universe. In particular the role of the scalar fields has been intensively investigated. From a qualitative point of view, since there is a dual representation of the three-form theory to a scalar field [15, 16] one may explore particular astrophysical settings describing objects such as oscillatons [69, 70] or maybe, their complex relatives, boson stars [71, 72], in the scalar representative frame. In the latter case, since boson stars are usually constructed from a complex scalar field (decomposed as two real scalar fields), the existence of the dual description from a single three-form field could be arguable. However, it is important to note that this dual nature, three-form \(\leftrightarrow \) scalar field, breaks down in some cases, even for fairly simple self-interactions [5], depending on the choice for the potential \(V(A^2)\). This fact, however, is not problematic, and, on the contrary, simply suggests that three-forms can provide us with new physics upon richer cosmological and astrophysical settings, undoubtedly worth exploring.

In particular, one may mention Bose–Einstein condensates consisting of ultralight bosons, and which can form localized and coherently oscillating stellar-type configurations [73]. For bosons having a higher mass, the bounded configurations may have typical masses and sizes of the order of magnitude similar to those of the neutron stars. Such objects formed from a primordial scalar field are known as boson stars (for reviews on the structure and properties of boson stars see [71] and [72]). Usually boson stars are constructed from a complex scalar field coupled to gravity, with the scalar field \(\phi \left( t,\vec {r}\right) \) decomposed into two real scalar fields \(\phi _R\left( t,\vec {r}\right) \) and \(\phi _I\left( t,\vec {r}\right) \) so that \(\phi \left( t,\vec {r}\right) =\phi _R\left( t,\vec {r}\right) +i\phi _I\left( t,\vec {r}\right) \) [72]. The evolution and properties of a boson star can be obtained from the action [72]

$$\begin{aligned} S_\mathrm{{BS}}=\int d^4 x \sqrt{-g}\left( \frac{1}{2\kappa ^2}R+\mathcal {L}_\mathrm{{SF}} \right) , \end{aligned}$$
(107)

where

$$\begin{aligned} \mathcal {L}_\mathrm{{SF}} = - \frac{1}{2} \left[ g^{\mu \nu } \nabla _{\mu } \bar{\phi }\, \nabla _{\nu } \phi + V\left( \left| \phi \right| ^2\right) \right] , \end{aligned}$$
(108)

where \(\bar{\phi }\) is the complex conjugate of the field, while \(V(|\phi |^2)\) is the scalar field potential depending only on the magnitude of the scalar field. There is a formal analogy between scalar field models, and the three form field model investigated in the present paper. In both cases the action is constructed from the standard Hilbert–Einstein term plus the field Lagrangian, which in our approach is given by Eq. (2). Hence if one could map in the complex plane the three-form Lagrangian \(\mathcal {L}_A\mapsto \mathcal {L}_{SF}\), then the present three form field extension of Einstein gravity would become equivalent with a scalar field model, and the corresponding solutions, assuming that the mapping three forms field \(\rightarrow \) scalar field does exist, would describe boson stars of different types.

In our analysis we have considered as a particular case of the three-form field potential the Higgs-type potential, given by Eq. (55). In the case of boson stars the gravitationally bounded configuration with the quartic self-interaction potential \(V\left( \left| \phi \right| ^2 \right) = m^2 \left| \phi \right| ^2 \, + \frac{\lambda }{2} \left| \phi \right| ^4\) was investigated in [74]. The physical properties of such a potential can be parameterized by the quantity \(\Lambda =\lambda M_\mathrm{Planck}^2/4\pi G m^2\), where \(M_\mathrm{Planck}\) is the Planck mass. The maximum mass of a boson star with quartic self-interaction potential is given by \( M_{\max } \approx 0.22 \Lambda ^{1/2} M_\mathrm{Planck}/m \) [72, 74]. It is also interesting to point out that in the Thomas–Fermi limit the scalar field becomes equivalent with a fluid, which in the low or moderate density limit has the equation of state \(P\propto \rho ^2\) [74].

By assuming that the three form field-scalar field correspondence is valid, by parameterizing the Higgs-type potential with the help of the parameter \(\Lambda =\nu /4\pi G\mu ^2\), the maximum mass of a stable three form field star is given by \(M_{\max }\propto \sqrt{\nu }/\sqrt{4\pi G}\mu ^3\), which could lead to masses of the same order of magnitude as the Chandrasekhar limit for neutron stars, \(M_{\max }\approx 3M_{\odot }\). However, the similarity (or equivalence) between three-form field stars and boson stars strongly relies on the fulfilling of the condition \(p_r\propto \rho ^2\), which, with the use of Eqs. (16) and (17) becomes

$$\begin{aligned} F^2-V\propto \left( F^2/48-V+\zeta V_{,\zeta }\right) ^2. \end{aligned}$$
(109)

If this differential equation has a solution for \(\zeta \), then the corresponding three form field stellar model has similar properties to a boson star. On the other hand, if no such solution for the radial function \(\zeta \) exists, then in the absence of ordinary matter the modified Einstein gravity in the presence of the three form fields admits only black hole-type solutions. We would like to point out that the above conclusions about the relations between boson and three form stars are mostly of qualitative nature, and in order to fully clarify these issues a detailed investigation of the properties of the three form field stars is required. From a technical point of view one can identify the formation of a star-like object from the conditions \(e^{\alpha \left( \eta _S\right) }\ne 0\) and \(e^{-\beta \left( \eta _S\right) }\ne 0\), respectively. Thus, as mentioned above, since the physical properties of the solutions of the Einstein field equations in the presence of a three form field are strongly dependent on the parameters \(\mu ^2\) and \(\nu \) of the Higgs potential the formation of stellar-type objects with no event horizon is possible.

In summary, the three-form field gravitational theory possess a rich mathematical structure, with theoretical properties that generate an intricate external dynamics. Consequently, the properties of the black holes in the three-form field theory are more complex as compared with the standard general relativistic black holes. These properties are related to the intrinsic properties of the three-form fields, which, from a mathematical point of view, lead to very complicated, strongly nonlinear, gravitational field equations. The new physical and geometrical effects generated by the presence of the three-form fields can also lead to some specific astrophysical and cosmological effects, whose observational detection may open new perspectives in the testing of gravitational theories. The astrophysical and observational implications of the black holes supported by a three-form field will be investigated in a future publication.