1 Introduction

The experimental detection of gravitational waves by the LIGO and VIRGO scientific collaborations [1, 2] has opened a new window on the intricate physical processes that control gravitational phenomena, leading to a better understanding of the properties of compact objects. The important GW170817 event [3] initiated the Multimessenger Era, with the signal, originating from the shell elliptical galaxy NGC 4993, and produced by the merging of two neutron stars, detected by more than 60 instruments worldwide. GW170817 suggests a constraint on the mass of the nonrotating neutrons star given by \(M\le 2.3M_{\odot }\). For a review of the merger of binary neutron star systems we refer the reader to [4]. The merger of neutron stars that takes place in the conditions of extreme gravity leads to the emission of intense fluxes of gravitational waves, associated with complicated microphysical and electromagnetic processes. The resulting astrophysical signatures can be observable even at the highest redshifts.

The detection of gravitational waves has also lead to some observations that are likely to change some basic paradigms in astrophysics. One such observation is related to the GW190814 event [5], which has shown a very intriguing structure of the mass components, with one of the individual masses in the range 2.5–2.67\(M_{\odot }\) (90% confidence). No optical counterpart to the gravitational wave was observed. If one assumes that this object is a neutron star, its high mass would strongly contradict the paradigm of the existence of a standard \( 1.4M_{\odot }\) mass scale for compact stars. In fact, a Bayesian statistical inference, performed in [6], evaluating the likelihood of the proposed Gaussian peaks by using 54 measured points obtained in a variety of systems, has already concluded to the existence of a bimodal distribution of the masses, with the first peak around 1.37 \(M_{\odot }\), and a much wider second peak at 1.73 \(M_{\odot }\).

However, the mass observed in the GW190814 event is much higher than these peaks. On the other hand, another accurate measurement of the mass of a compact object, using Shapiro delay, yielded the value \(2.14 ^{+ 0.10}_{ -0.09}M_{\odot }\) for the mass of the millisecond pulsar MSP J0740+6620 [7]. A few higher measured mass values also exist. These observations raise the question of the maximum mass \( M_\mathrm{max}\) of stable compact general relativistic objects, since in order to accommodate the observed values one must significantly enlarge the allowed range of \(M_\mathrm{max}\). In turn, this would require important modifications for the physical properties of the dense matter at densities higher than the saturation density, and a corresponding modification of the equation of state. Even more interesting, and with important theoretical consequences, is the recent discovery of a companion, having a mass of around \(3M_{\odot }\), of V723 Mon, a nearby evolved red giant in a high-mass function, \(f(M)=1.72\pm 0.01M_{\odot }\), nearly circular binary [8]. If confirmed, this discovery will raise new questions about present-day knowledge of the mass distribution of massive compact objects, and on the transition of neutron stars to black holes.

There are a number of effects that could lead to the increase of the masses of compact stars. One such effect is related to the assumption of a phase transition in the dense matter, leading to the formation of a quark star. In [9] the possibility that stellar mass black holes, with masses in the range of 3.8\(M_{\odot }\) and 6\(M_{\odot }\), could be in fact quark stars in the Color–Flavor Locked (CFL) phase was investigated in detail. It was shown that, depending on the value of the gap parameter, rapidly rotating CFL quark stars may have much higher masses than ordinary neutron stars. On the other hand quark stars have a very low luminosity, and an almost completely absorbing surface, due to the fact that infalling matter on the surface of the quark star is fully transformed into quark matter.

A possibility of distinguishing quark stars in standard or CFL phase from low-mass black holes or neutron stars could be through the study of thin accretion disks around rapidly rotating stars (neutron or quark), and Kerr-type black holes, respectively. It was already suggested that the GW190814 event resulted from the merging of a black hole–strange quark star system [10,11,12]. Some compact astrophysical objects may contain a significant part of their matter in the form of a Bose–Einstein condensate. The basic astrophysical parameters (mass and radius) of the neutron stars sensitively depend on the mass of the condensed particle, and on the scattering length. Hence the recently observed neutron stars with masses in the range of 2–2.5 \(M_{\odot }\) could be Bose–Einstein condensate stars, containing a large amount of superfluid matter [13].

An alternative explanation of the higher masses of some classes of neutron stars is related to the possible modification of the very nature of the gravitational force at very high densities, implying that the structure of neutron stars is described by some modified theories of gravity. In the presence of a cosmological constant \(\Lambda \) the mass–radius M/R ratio of compact objects satisfy the constraint \(2M/R\le \left( 1-8\pi \Lambda R^2/3 \right) \left[ 1-\left( 1-2\Lambda /{\bar{\rho }}\right) ^2/9\left( 1-8\pi \Lambda R^2/3\right) \right] \) [14]. Upper and lower bounds on the mass–radius ratio of stable compact objects in extended gravity theories, in which modifications of the gravitational dynamics with respect to standard general relativity are described by an effective contribution to the matter energy-momentum tensor \(T_{\mu }^{\nu }\), given by a tensor \(\theta _{\mu }^{\nu }\), were derived in [15]. By introducing the effective density inside the star defined as \(\rho _{\mathrm {eff}}c^{2} = \rho c^{2}/G + \theta _{0}^{0}\), and the effective mass \(m_{\mathrm {eff}}=4\pi \int _0^r{r^2\rho _{\mathrm {eff}}dr}\), one obtains for the generalized Buchdahl limit for extended gravitational theories the expression

$$\begin{aligned} \frac{2m_{\mathrm {eff}}(r)}{r}\le 1-\left[ 1+\frac{2\left( 1+f(r)\right) }{1+4\pi w_{\mathrm {eff}}(r)}\right] ^{-2}, \end{aligned}$$
(1)

where

$$\begin{aligned} f(r) = 4\pi \frac{\Delta (r)}{\langle \rho _{\mathrm {eff}}\rangle (r)}\left\{ \frac{\arcsin \left[ \sqrt{2m_{\mathrm {eff}}(r)/r}\right] }{\sqrt{2m_{\mathrm {eff}}(r)/r}} -1\right\} , \end{aligned}$$
(2)

\(w_{\mathrm {eff}}(r)=p_{\mathrm {eff}}/\langle \rho _{\mathrm {eff}}\rangle (r)\), and \(\Delta =\left( G/c^4\right) \left( \theta _1^1-\theta _2^2\right) \), respectively. Hence, the extra contributions to the matter energy-momentum tensor due to the modifications of the gravitational force could lead to a significant increase in the mass of the neutron star.

The structure and physical properties of specific classes of neutron, quark and “exotic” stars in Eddington-inspired Born–Infeld gravity were considered in [16]. The latter reduces to standard general relativity in vacuum, but presents a different behavior of the gravitational field in the presence of matter. The equilibrium equations for a spherically symmetric configuration (mass continuity and Tolman–Oppenheimer–Volkoff) were derived, and their solutions were obtained numerically for different equations of state of neutron and quark matter. The internal structure and the physical properties of specific classes of neutron, quark and Bose–Einstein condensate stars in the hybrid metric-Palatini gravity theory [17,18,19,20,21], which is a combination of the metric and Palatini f(R) formalisms, was considered in [22]. As a general result it was found that, for all the considered equations of state, hybrid metric-Palatini gravity stars are more massive than their general relativistic counterparts. The properties of neutron stars in f(RT) gravity [23] for the case \(R+2\lambda T\), where R is the Ricci scalar and T is the trace of the energy-momentum tensor were investigated in [24]. The hydrostatic equilibrium equations have been solved by considering realistic equations of state. It was also found that using several relativistic and non-relativistic models the variation on the mass and radius of the neutron star is almost the same for all considered equations of state, indicating that the results are independent of the high density part of the equation of state. Hence the stellar masses and radii depend only on the crust, where the equation of state is essentially the same for all the models.

The above-mentioned modified theories of gravity can be represented in a respective scalar-tensor theory. In fact, scalar fields in cosmology and gravitation have unquestionably played a crucial role over the last decades [25, 26]. Nonetheless, so far we have only detected one scalar particle in nature, namely, the Higgs boson [27]. On a smooth manifold \(\mathcal {M}\), scalar fields are part of a more general class of fields, denoted as n-forms, with scalars being the \(n=0\) case. These differential forms inhabit smooth sections of the nth exterior power of the cotangent bundle \({T^* \mathcal {M}\overset{\pi }{\longrightarrow } \mathcal {M}}\), and thus naturally exist in most geometrical settings. When performing calculus on manifolds, these fields have an enormous importance, since their algebraic structure allows us to construct diffeomorphism invariant objects and carry out integration in a coordinate-free fashion. Pioneered by Élie Cartan in the beginning of the 20th century, they are also the central objects of De Rham cohomology, linking the topological and differentiable structures of manifolds. Since the geometrical anatomy of a manifold is intimately connected with the dynamical behavior of theories defined within it, it is reasonable to assume that differential forms play a fundamental role in physics, particularly while formulating (gauge) field theories in a form-representative framework. Therefore it becomes relevant to explore the effects of considering higher spin fields in gravitation and cosmology.

Here, we will focus on three-form fields [28]. Three-forms naturally exist in fundamental theories, such as string theory and supergravity [29,30,31,32,33], so it is not unreasonable to expect their emergence in low energy effective actions. Their first link to cosmology may be traced back to Hawking and Turok [34], when trying to explain the tiny value of the cosmological constant and thus realizing that an action encompassing a four-form, constructed from a three-form gauge field, behaves exactly as a cosmological constant. Ten years later, models of three-forms with the addition of self-interacting potentials were studied, where it was shown that these may give rise to self-accelerating attractor solutions [35], useful to explain primordial inflation [36,37,38,39] and dark energy [40,41,42]. These models also present distinct observational signatures, so, in principle, it would be possible to distinguish between three-form and standard scalar driven models [43].

Three-forms were further applied to other scenarios, such as a mechanism for magnetogenesis [44] considering U(1) couplings to the electromagnetic field, and three-form screening mechanisms around dense objects [45]. Three-forms cosmological models were analyzed in detail in [46], where it was shown that in this kind of models there are fixed points at infinity. This result was obtained by introducing appropriate compactifications, as well as by defining a new time variable that eliminates any potential divergence of the system. Normally hyperbolic non-isolated fixed points were also identified.

The employment of three-form fields in gravitation was analyzed in [47], where wormhole solutions in a static and spherically symmetric spacetime continuum were explored. It was found that it is possible for the ordinary matter fields threading the wormhole to obey the classical energy conditions throughout the spacetime, while the three-form field holds the wormhole open, violating the null and weak energy conditions. It is known [37] that a three-form admits a dual scalar field representation. However, it is important to notice that this mapping easily breaks down when considering even fairly simple self interactions, nonminimal couplings, or noncanonical kinetic terms for the three-form [35]. When this dual representation is well defined, the scalar representation is often complicated and untractable. Hence the three-form formalism may provide for new and rich frameworks to explore physical phenomena. For a model of three-form fields in a homogeneous and isotropic universe a canonical quantization procedure of the Wheeler–DeWitt type was discussed in [48], by obtaining first the Hamiltonian description of the model. This formalism was applied to a Little Sibling of the Big Rip (LSBR), an event that is known to appear in several minimally coupled three-form fields for a variety of potentials. A set of analytical solutions of the Wheeler–DeWitt equation was obtained, their physical implications were discussed. It was also shown that there are quantum states where the wave function of the universe vanishes, i.e. the DeWitt condition is fulfilled for them.

The static and spherically symmetric vacuum solutions in the three-form field gravity theory were also investigated in [49]. For the case of the vanishing three-form field potential the gravitational field equations can be solved exactly. However, for arbitrary potentials, due to their mathematical complexity, numerical approaches were adopted for studying the behavior of the metric functions and the three-form field. The formation of a black hole was detected from the presence of a Killing horizon for the time-like Killing vector in the metric tensor components. Several models, corresponding to different functional forms of the three-field potential, namely, the Higgs and exponential type, were considered. In particular, naked singularity solutions were also obtained for the exponential potential case. The thermodynamic properties of the black hole solutions, such as the horizon temperature, specific heat, entropy and evaporation time due to the Hawking luminosity, were investigated in detail. In [50] it was shown that a minimally coupled 3-form endowed with a proper potential can support a regular black hole interior, and that by choosing an appropriate form for the potential of the 3-form field, one can construct an interior black hole spacetime that is everywhere regular. Moreover, the singularity is replaced with a Nariai-type spacetime, whose topology is \(dS_2\times S_2\). Despite the negative potential of the 3-form field, the 2-dimensional de Sitter geometry always appears. Due to the violation of the null energy condition this particular geometry is singularity-free.

As an extension of the above-mentioned work, it is the goal of the present paper to consider the properties of compact high density stars in the three-form field gravity theory. To simplify the mathematical formalism we adopt a spherically symmetric static geometry, and a matter source. After writing down the gravitational field equations in their full generality, as a first step in our study we obtain the generalized mass continuity equation, the Tolman–Oppenheimer–Volkoff equation, describing the hydrostatic equilibrium of the star, as well as the field equation describing the variation inside the compact object of the non-zero component of the three-form field. These three equations fully describe the macroscopic properties of the three-form field stars. The system of the structure equations of the three-form field gravity theory is then solved numerically for several equations of state of the high density matter.

In our study we consider three-form field stars whose ordinary baryonic matter content is described by the stiff-fluid (Zeldovich) equation of state, satisfying the causality condition requiring that the speed of sound in the dense matter does not exceed the speed of light, by the photon gas equation of state, describing a radiation fluid with the property, that the trace of the energy-momentum tensor identically vanishes; the strange quark matter equation of state, derived from the MIT bag model, and, the equation of state of the Bose–Einstein condensates, given by a polytropic equation of state with polytropic index \(n = 1\). For all these equations of state of the dense matter the astrophysical parameters of the stars (density, pressure and mass distribution, radius and total mass), as well as the behavior of the three-form field and of its potential are obtained numerically, and compared with the results of similar standard general relativistic models. Hence this strategy allows us to perform an in-depth comparison of the two gravitational theories, and of the impact of modified gravity on the description of the properties of the stellar structures. We can formulate a general conclusion of our study by pointing out that three-form field gravity theory predicts the existence of more massive high density stars, as compared to standard general relativity.

The present paper is organized as follows. The action and the field equations of the three-form field gravity theory is briefly introduced in Sect. 2. The system of gravitational field equations, describing the interior of static spherically symmetric stars, are written down in Sect. 3, where the structure equations of compact objects (mass continuity, Tolman–Oppenheimer–Volkoff, and three-form field evolution equation) are also derived, and rewritten in a dimensionless form. The global astrophysical properties of the Zeldovich (stiff-fluid), photon, strange matter, and Bose–Einstein condensate stars are obtained, by numerically integrating the structure equations of the static spherically symmetric three-form field gravity theory, in Sect. 4. We discuss and conclude our results in Sect. 5.

2 Action and field equations

Let us start by considering a single canonical three-form field \(\mathbf{A}\), with action given by [28, 35]

$$\begin{aligned} \mathcal {S}_A = - \int \omega _g\left[ \frac{1}{48}F^{\alpha \beta \gamma \delta }F_{\alpha \beta \gamma \delta } + V(A^{\alpha \beta \gamma } A_{\alpha \beta \gamma }) \right] , \end{aligned}$$
(3)

where \(\omega _g = \sqrt{-g}\,\mathrm{d}^4x\) is the metric volume form, with g being the determinant of the metric, and we have assumed a self-interacting potential V, which is a function of the invariant \(A^{\alpha \beta \gamma } A_{\alpha \beta \gamma }=:A^2\,\). The first kinetic term in Eq. (3),

$$\begin{aligned} -\frac{1}{48} \int \omega _g F^{\alpha \beta \gamma \delta }F_{\alpha \beta \gamma \delta } = -\frac{1}{2} \int \mathbf {F}\wedge \star \mathbf {F}, \end{aligned}$$
(4)

comprises the 4-form strength tensor \(\mathbf{{F}=d\mathbf {A}\,}\) [41, 43], which generalizes Maxwell’s 2-form from classical electromagnetism, with components \({F_{\alpha \beta \gamma \delta } = 4\nabla _{[\alpha }A_{\beta \gamma \delta ]}\,}\), and it is a closed form, \({\mathbf{{dF}}=0\,}\). Note that, since the action Eq. (3) bears a self-interaction term, \(V(A^2)\), gauge invariance may be broken. However, it is possible to restore this symmetry through the introduction of a Stückelberg form [35, 37].

The total action of our theory minimally coupled to Einstein’s gravity, can thus be written as

$$\begin{aligned} \mathcal {S} = \int \omega _g \left[ \frac{1}{16\pi G}R + \mathcal {L}_m (g_{\mu \nu },\psi ,\nabla \psi ) \right] +\mathcal {S}_A, \end{aligned}$$
(5)

where G is Newton’s constant, R the curvature scalar and \(\mathcal {L}_m\) stands for an anisotropic distribution of matter threading our spactime. The field action \(\mathcal {S}_A\) is given by \({\mathcal {S}_A = \int \omega _g\,\mathcal {L}_A}\), where \(\mathcal {L}_A\) is the field’s Lagrangian density.

To find the dynamics governing this model, we start by varying the action Eq. (5) with respect to the three-form, and we find the following equations of motion:

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

The energy-momentum tensor relative to the three-form reads

$$\begin{aligned} T^{(A)}_{\mu \nu }= & {} -2 \frac{\delta \mathcal {L}_{A}}{\delta g^{\mu \nu }}+g_{\mu \nu }\mathcal {L}_{A} \nonumber \\= & {} \frac{1}{6}F_{\mu \alpha \beta \gamma }F_{\nu }{}^{\alpha \beta \gamma } + 6 \frac{\partial V}{\partial (A^2)} A_{\mu \alpha \beta }A_{\nu }{}^{\alpha \beta } +\mathcal {L}_A\,g_{\mu \nu }.\qquad \nonumber \\ \end{aligned}$$
(7)

Note, however, that the equations of motion (6) could equivalently be deduced from the contracted Bianchi identities [42],

$$\begin{aligned} \nabla _{\mu }T^{(A)\,\mu }{}_{\nu }=\frac{1}{6}F_{\nu \alpha \beta \gamma }E^{\alpha \beta \gamma }=0. \end{aligned}$$
(8)

The modified field equations can be computed by varying Eq. (5) with respect to \(g^{\mu \nu }\), leading to

$$\begin{aligned} G_{\mu \nu } = 8\pi G\left( T^{(m)}_{\mu \nu } + T^{(A)}_{\mu \nu } \right) , \end{aligned}$$
(9)

with \(G_{\mu \nu }\) being the components of the Einstein tensor.

Taking the divergence of Eq. (9) and using the field equation (6) one can find that the energy-momentum tensor of matter is conserved,

$$\begin{aligned} \nabla ^\mu \, T^{(m)}_{\mu \nu }=0. \end{aligned}$$
(10)

Regarding the matter source, we assume an anisotropic distribution of matter, where the components of the energy-momentum tensor can be written as

$$\begin{aligned} T^{(m)}_{\mu \nu } = (\rho _m+p_{\perp m})u_\mu u_{\nu } + p_{\perp m} g_{\mu \nu } + (p_{rm}-p_{\perp m})\chi _{\mu }\chi _{\nu },\nonumber \!\!\!\!\!\!\!\\ \end{aligned}$$
(11)

\(u_{\mu }\) being the four-velocity, normalized as \(u_{\mu }u^{\mu }=-1\), \(\chi ^{\mu }\) the unit spacelike vector in the radial direction, \(\rho _m\), \(p_{rm}\) and \(p_{\perp m}\) the energy density, radial pressure and tangential pressure, respectively.

Note that this present theory considers standard General Relativity for the gravitational Lagrangian with the novelty being the introduction of a 3-form field fluid in the matter sector which modifies the overall dynamics through the intimate relation of the gravitational field equations (9). Thus, we will informally refer to this model as three-form gravity.

3 Spherically symmetric and static background

This work aims at finding solutions for three-form supported stars. To this effect, let us consider the following static and spherically symmetric line element [51]:

$$\begin{aligned} {\text {d}}s^2 =- & {} \exp \left[ -2\int _r^{\infty }h(\tilde{r}){\text {d}}\tilde{r}\right] {\text {d}}t^2 + \frac{{\text {d}}r^2}{1-2\, G \,m(r)/r} \nonumber \\+ & {} r^2 \left( {\text {d}}\theta ^2 + \sin ^2\theta \,{\text {d}}\phi ^2 \right) , \end{aligned}$$
(12)

where m(r) is the mass function, and h(r) is the gravity profile, both functions of the radial coordinate, r, only.

Similar to previous studies [47, 49] on static and spherically symmetric backgrounds, we will assume an ansatz for our three-form field through the aid of its Hodge dual 1-form \({\mathbf{B}}=\star {\mathbf{A}}\), which fully characterize the components of the field, i.e.,

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

Hence, we choose a convenient parametrization of \({\mathbf{B}}\) as the radially directed vector with components,

$$\begin{aligned} B^{\delta }=\left( 0,\sqrt{1-\frac{2Gm(r)}{r}}\zeta (r),0,0 \right) \,. \end{aligned}$$
(14)

dependent on a suitable scalar function \(\zeta (r)\,\). Accordingly, the dynamical equations for the three-form will be expressed entirely in terms of \(\zeta (r)\). Since this scalar function fully determines the components of the three-form, by way of Eq. (13), throughout this manuscript we may naively refer to \(\zeta \) simply as the three-form, although, formally, it is the scalar function determining the components of the three-form, \(A_{\alpha \beta \gamma }\). We refer the reader to Appendix A of [36] to examine how the action Eq. (3) may be written solely in terms of the dual vector B, highlighting how this model can be recast into a vector-tensor theory at the background level.

In this gravitational setting, the contraction \(A^2\) becomes

$$\begin{aligned} A^2 = -6B^{\delta }B_{\delta } = -6 \,\zeta (r)^2, \end{aligned}$$
(15)

and the kinetic invariant

$$\begin{aligned} F^2:= & {} F^{\alpha \beta \gamma \delta }F_{\alpha \beta \gamma \delta } = -24\,(\nabla _{\delta }B^{\delta })^2 \nonumber \\= & {} -24\left( 1-\frac{2Gm}{r} \right) \left[ \zeta \left( \frac{2}{r} - h \right) + \zeta ' \right] ^2, \end{aligned}$$
(16)

where derivatives with respect to the radial coordinate r are denoted by a prime and \(\zeta = \zeta (r)\). We may now express the dynamics of the three-form, Eq. (6), using the metric, Eq. (12), solely in terms of \(\zeta \) and the metric functions, through

$$\begin{aligned}&r^2\zeta ^{\prime \prime }\left( \frac{2Gm}{r}-1 \right) +r\,\zeta ^{\prime } \left[ G m^{\prime }+\frac{3Gm}{r}-2 \right. \nonumber \\&\quad \left. +rh\left( 1-\frac{2Gm}{r} \right) \right] +\zeta \left[ 2+2G \left( m^{\prime }- \frac{3m}{r}\right) \right. \nonumber \\&\quad \left. +rhG\left( \frac{m}{r}-m^{\prime }\right) +r^2 h^{\prime }\left( 1-\frac{2Gm}{r} \right) \right] -r^2V_{,\zeta }=0,\nonumber \\ \end{aligned}$$
(17)

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

The components of the energy-momentum tensor of the three-form can be computed by plugging Eq. (12) into Eq. (7), which leads to

$$\begin{aligned} T^{(A)\,\,t}{}_{t}= & {} -\rho _A = \frac{1}{48}F^2 - V + \zeta V_{,\zeta }, \end{aligned}$$
(18)
$$\begin{aligned} T^{(A)\,\,r}{}_{r}= & {} p_{rA} = \frac{1}{48}F^2 - V, \end{aligned}$$
(19)
$$\begin{aligned} T^{(A)\,\,\theta }{}_{\theta }= & {} T^{(A)\,\,\phi }{}_{\phi } = p_{\perp A} = T^{(A)\,\,t}{}_{t}, \end{aligned}$$
(20)

with the object \(F^2\) given by Eq. (16).

The components of the modified field equations, Eqs. (9), under such a gravitational framework yield

$$\begin{aligned} \rho _A + \rho _m= & {} \frac{m^{\prime }}{4\pi r^2}, \end{aligned}$$
(21)
$$\begin{aligned} p_{rA}+p_{rm}= & {} - \frac{m}{4\pi r^3}-\frac{h}{4\pi G r}\left( 1- \frac{2G m}{r}\right) , \end{aligned}$$
(22)
$$\begin{aligned} p_{\perp A}+p_{\perp m}= & {} \frac{1}{8\pi G}\left( 1-\frac{2G m}{r} \right) \left[ \left( h-\frac{1}{r} \right) \right. \nonumber \\&\left. \times \left( h+\frac{Gm'r-Gm}{r(r-2Gm)} \right) -h'\right] . \end{aligned}$$
(23)

In addition to this, the conservation equation (10) takes the form

$$\begin{aligned} p_{rm}^{\prime }=\frac{2}{r}\left( p_{\perp m}-p_{rm}\right) +h\left( \rho _m+p_{rm}\right) . \end{aligned}$$
(24)

By eliminating h between Eqs. (22) and (24) we obtain the generalized Tolman–Oppenheimer–Volkoff (TOV) equation describing the structure of massive compact objects in three-form gravity, given by

$$\begin{aligned} p^{\prime }_{rm}= & {} -\frac{G\left( \rho _m+p_{rm}\right) \left[ \left( p_{rA}+p_{rm}\right) 4\pi r^3+m\right] }{r^2\left( 1-2Gm/r\right) } \nonumber \\&+\frac{2}{r}\left( p_{\perp m}-p_{rm}\right) . \end{aligned}$$
(25)

The TOV equation, together with the mass continuity equation (21), and the equation of motion for \(\zeta \), form a system of differential nonlinear ordinary equations, whose solution fully characterize star-like objects. The system must be considered with the initial conditions \(m(0)=0\), \(p_{rm}(0)=p_{rm}^{(0)}\), and \(\zeta (0)=\zeta _0\), respectively. In order to close the system of equations the equations of state of the dense matter and the functional form of the potential \(V\left( A^2\right) \) must also be specified.

We want to express the results in units of km and \(M_{\odot }\), and hence we use a set of dimensionless quantities defined as

$$\begin{aligned}&\rho _m= \epsilon _0 {\bar{\rho }}_m, \quad p_{rm}= \epsilon _0 {\bar{p}}_{rm}, \quad p_{\perp m}= \epsilon _0 {\bar{p}}_{\perp m},~~ h={\bar{h}}/R_0, \nonumber \\&m=M_\odot \,{\bar{m}}, \quad \zeta ={\bar{\zeta }}\left( M_\odot /R_0\right) ^{1/2},\quad r= R_0 \,\eta , \end{aligned}$$
(26)

where \(\epsilon _0\) is an arbitrary energy-density scale, \(M_\odot \) is solar mass, and \(R_0=G M_\odot /c^2=1.477 ~\mathrm{km}\). In the following we set \(\epsilon _0=M_\odot c^2/80 \pi R_0^3=2.45\times 10^{15} \, \mathrm{g/cm}^3\). One should note that according to the above definitions we have

$$\begin{aligned}&\rho _A= \epsilon _0 {\bar{\rho }}_A, \qquad p_{rA}= \epsilon _0 {\bar{p}}_{rA}, \qquad p_{\perp A}= \epsilon _0 {\bar{p}}_{\perp A}. \end{aligned}$$
(27)

In the dimensionless variables defined above the system of the structure equations of the three-form stars take the form

$$\begin{aligned} \frac{{\text {d}}{\bar{m}}}{{\text {d}}\eta }= & {} \frac{1}{20}\left( {\bar{\rho }}_{A}+{\bar{\rho }}_{m}\right) \eta ^{2}, \end{aligned}$$
(28)
$$\begin{aligned} \frac{{\text {d}}{\bar{p}}_{rm}}{{\text {d}}\eta }= & {} -\frac{\left( {\bar{\rho }}_{m}+{\bar{p}} _{rm}\right) \left[ \left( {\bar{p}}_{rA}+{\bar{p}}_{rm}\right) \eta ^{3}/2+ 20 \,\bar{m}\right] }{20\, \eta ^{2}\left( 1-2{\bar{m}}/\eta \right) } \nonumber \\&+\frac{2}{\eta }\left( {\bar{p}}_{\perp m}-{\bar{p}}_{rm}\right) , \end{aligned}$$
(29)
$$\begin{aligned}&\!\!\!\eta ^{2}\frac{{\text {d}}^{2}{\bar{\zeta }} }{{\text {d}}\eta ^{2}}\left( \frac{2{\bar{m}}}{\eta } -1\right) \nonumber \\&\quad +\eta \frac{{\text {d}}{\bar{\zeta }} }{{\text {d}}\eta }\,\left[ \frac{{\text {d}}{\bar{m}}}{d\eta }+\frac{3{\bar{m}}}{\eta }-2+\eta {\bar{h}}\left( 1-\frac{2{\bar{m}}}{\eta }\right) \right] \nonumber \\&\quad +{\bar{\zeta }} \left[ 2+2\frac{{\text {d}}{\bar{m}}}{{\text {d}}\eta }-\frac{6{\bar{m}}}{\eta }+\eta {\bar{h}} \left( \frac{{\bar{m}}}{\eta }-\frac{{\text {d}}{\bar{m}}}{{\text {d}}\eta }\right) \right. \nonumber \\&\quad \left. +\eta ^{2}\frac{{\text {d}}{\bar{h}}}{{\text {d}}\eta }\left( 1-\frac{2{\bar{m}}}{\eta }\right) \right] -\frac{1}{80 \pi }\frac{{\text {d}}{\bar{V}}\left( {\bar{\zeta }}\right) }{{\text {d}}{\bar{\zeta }}} \eta ^{2}=0. \end{aligned}$$
(30)

The system of Eqs. (28)–(30) must be integrated with the boundary conditions \({\bar{m}}(0)=0\), \({\bar{p}}_{rm}(0)=p_{rm}(0)/\epsilon _0\), \({\bar{p}}_{\perp m}(0)=p_{\perp m}(0)/\epsilon _0\) etc., \({\bar{\zeta }}(0)={\bar{\zeta }} _0\), \({\bar{\zeta }}'(0)={\bar{\zeta }} ^{\prime }_0\), and \({\bar{p}}_{rm}(R)=0\), respectively, with the latter condition determining the radius of the compact stellar objet in the three-form fields gravity theory.

In the following we are going to construct specific models of three-form field stars under the assumption that the three-form field vanishes at the center of the star, and increases towards the surface, where it reaches its maximum value. Therefore, we impose the condition \(\zeta (0)=0\), indicating that the field becomes negligibly small near the center. Hence, from a physical point of view we assume that the three-form field, as well as the associated physical parameters, do follow, at least qualitatively, the variation of the mass distribution inside the star, and that the maximum values of the field are related to the maximum values of the mass of the ordinary matter. Moreover, we assume a slow increase of the field near the center of the star, with \(\zeta '(0)\ne 0\), with the field derivative also taking a small value near the center \(r=0\), so that \(\zeta '(0) \ll 1\). Indeed, different types of models in which the field is, for example, maximum at the center of the star, or it has some arbitrary large values at \(r=0\), can also be constructed.

4 Structure and astrophysical properties of compact objects

In the present section, we will investigate the basic astrophysical properties of high density compact objects in the three-form field gravity. In our analysis we will not impose any specific restrictions on the functional form of the three-form field \(\zeta \), which is assumed to obey Eq.  (30). In order to close the system of field equations describing the interior of the stellar objects we need to adopt an equation of state for the ordinary matter. Despite intensive theoretical efforts the equation of state of high density baryonic matter is poorly known, and a realistic description of the physical properties of matter at densities exceeding with several orders of magnitude the nuclear density \(\rho _n=2\times 10^{14}\) g/cm\(^3\) is still missing. For this reason we consider some simple approximations of the matter equation of state, by considering four specific cases, given by: (i) the stiff-fluid equation of state, in which the pressure equals the energy density, \(P=\rho \); (ii) the radiation fluid equation of state corresponding to \(P = \rho /3\); (iii) the important quark matter equation of state \(P= \left( \rho -4B\right) /3\), where B is the bag constant, and (iv) the superfluid neutron matter Bose–Einstein condensate equation of state, which is given by a polytropic equation of state with the polytropic index \(n=1\), so that \(P\propto \rho ^2\).

An important physical quantity with a profound influence on the stellar structure is the three-form potential \(V\left( A^2\right) \). In the following we will restrict our investigations to two functional forms of \(V\left( A^2\right) \). We will assume that it is either a constant, or it has a Higgs-type form, with \(V\left( A^2\right) =aA^2+bA^4\), respectively, where a and b are constants. From a physical point of view the constant \(a<0\) is related to the mass of the three-form field via the relation \(m_{\zeta }^2= bv^2/2=-2a\), where \(v^2=-a/2b\) gives the minimum of the Higgs-type potential.

4.1 Stellar structures with stiff-matter fluid

As a first example of a stellar model in the three-form field gravity we consider the case of an isotropic distribution of matter, with \(p_{rm}=p_{\perp m}\), with the matter pressure satisfying the stiff-fluid (Zeldovich) equation of state [52, 53]

$$\begin{aligned} p_{rm}=\rho _m. \end{aligned}$$
(31)

From a physical point of view such an equation of state may describe the matter behavior at densities ten times higher than the nuclear density, i.e., at densities greater than or equal to \(10^{17}\) g/cm\(^3\), corresponding to temperatures \(T=(\rho /\sigma )^{ 1/4} >10^{13}\) K, where by \(\sigma \) we have denoted the radiation constant [53]. An important characteristic of the Zeldovich equation of state is that for stiff matter the speed of sound is equal to the speed of light, \(c_s^2=\partial p_{rm}/\partial \rho _m=1\). One reason why such an equation of state was proposed is that for stiff matter the matter perturbations cannot move at speeds greater than the speed of light. The stiff-matter equation has found many applications in astrophysics. In [54] it was shown that for the equilibrium configuration of a stellar object with a very high density the maximum mass cannot surpass the upper limit value of \(3.2M_{\odot }\). To obtain this fundamental result it was assumed that the stellar interior is described by the spherically symmetric static Einstein field equations, and that the principle of causality and Le Chatelier’s principle both hold. Moreover, it was assumed that at very high densities the neutron matter satisfies the stiff-fluid equation of state \(p=\rho \). This important result on the numerical value of the maximum mass of stable compact objects offers a fundamental criterion for observationally distinguishing ordinary neutron or other types of compact stars from black holes, and it represents an important prediction of the theory of general relativity.

In considering the properties of three-form field stars described by the stiff-matter equation of state we assume that the central density of the compact object varies between the values \(3.1\times 10^{14}\;\mathrm{g/cm}^3\) and \(2.2\times 10^{15}\; \mathrm{g/cm}^3\). We stop the numerical integration when the density reaches the surface value \(\rho _{m}= 10^{14}\; \mathrm{g/cm}^3\), a density lower than the nuclear density. The initial value of the three-form field and of its derivative is \({\bar{\zeta }}_0=0\) and \({\bar{\zeta }}^\prime _0=7\times 10^{-3}\) in all considered cases. In the following, we present the results of the numerical integration for stiff-fluid stars obtained by considering two different choices of the three-form field potential to obtain the interior solutions.

4.1.1 The constant potential case: \(V(A^2)=\lambda \)

As a first case of a stellar type structure constructed in the three-form field gravity we assume that the field potential is a constant, \(V(A^2)=\lambda =\mathrm{constant}\). To describe the properties of the potential we introduce the dimensionless parameter \({\bar{\lambda }}\), defined as

$$\begin{aligned} \lambda =\epsilon _0 {\bar{\lambda }}. \end{aligned}$$

The variations with respect to the dimensionless coordinate \(\eta \) of the energy density of the matter and of the mass of the three-form stars are represented, for different values of \({\bar{\lambda }}\), in Fig. 1. As one can see from the plots, the matter energy density is a monotonically decreasing function of the radial coordinate, and it tends towards the zero value at the star surface. The dependence on the numerical values of the potential \({\bar{\lambda }}\) is weak, at least for small values of the radial coordinate. The interior mass profile of the star is a monotonically increasing function of \(\eta \), and it shows a significant dependence on \({\bar{\lambda }}\).

Fig. 1
figure 1

The interior mass–energy-density profile \({\bar{\rho }}_m\) (left panel) and the mass profile \({\bar{m}}\) (right figure) as a function of the distance from the center of the three-form stiff-fluid stars \(\eta \), in the presence of a constant field potential, for three different values of the constant \(\bar{\lambda }\): \({\bar{\lambda }}=0\) (dashed curve), \({\bar{\lambda }}=0.02\) (dotted curve), and \(\bar{\lambda }=-0.02\) (dot-dashed curve), respectively. For the central density of the star we have adopted the value \(\rho _{{\text {mc}}}=3.7\times 10^{14} \mathrm{g/cm}^3\). The solid curve represents the standard general relativistic density and mass profiles for stiff-fluid stars

Fig. 2
figure 2

Variation of the non-zero component \({\bar{\zeta }}\) (left panel) and of the effective energy density \({\bar{\rho }}_A\) (right figure) of the three-form field as a function of the radial distance from the center of the stiff-fluid star \(\eta \) with constant potential, for three different values of the constant \({\bar{\lambda }}\): \({\bar{\lambda }}=0\) (dashed curve), \({\bar{\lambda }}=0.02\) (dotted curve), and \(\bar{\lambda }=-0.02\) (dot-dashed curve), respectively. For the central density of the star we have adopted the value \(\rho _{{\text {mc}}}=3.7\times 10^{14} \mathrm{g/cm}^3\), while \({\bar{\zeta }}_0=0\) and \({\bar{\zeta }}^\prime _0=7\times 10^{-3}\), respectively

The variations of the three-form field \({\bar{\zeta }}\) and of the energy density of the three-form fields for the stiff-fluid stars in the presence of a constant potential are depicted in Fig. 2. The three-form field \({\bar{\zeta }}\) is a monotonically increasing function inside the star, and it reaches its maximum value on the star surface. On the other hand, inside the star the energy density of the three-form field is a constant, independently of the adopted values of the parameters. Moreover, depending on the numerical values of the three-form field potential, \({\bar{\rho }}_A\) can take both negative and positive values.

The mass–radius relations of the three-form stiff-fluid stars with constant potential are shown in Fig. 3. The general relativistic case is also represented for comparison. As one can see from the figure, the mass–radius relation for compact objects shows significant differences as compared to the standard general relativistic case. Much higher maximum masses, of the order of \(4.1M_{\odot }\), exceeding the maximum mass limit of \(3.2M_{\odot }\) of general relativity, can be achieved for positive values of the three-form field potential. On the other hand negative values of \({\bar{\lambda }}\) lead to a decrease of the maximum mass of the neutron stars to a value of the order of \(2.5M_{\odot }\), and generally of the values of the stellar masses.

Some specific numerical values of the maximum masses of stiff-fluid stars in three-form field gravity are presented, for different values of the field potential, in Table 1. In general relativity for stiff-fluid stars we have \({M_{\max }=3.28\,M_{\odot }}\), \(R=18.53\,\mathrm{km} \), and \({\rho _{{\text {mc}}}=1.21\times 10^{15}\, \mathrm{g/cm}^3}\), respectively.

4.1.2 Higgs-type potential: \(V(A^2)=a A^2 +b A^4\)

For the case of stiff-fluid stars with a Higgs-type three-form field potential given by \(V(A^2)=a A^2 +b A^4\), we rescale the potential parameters to a dimensionless form so that

$$\begin{aligned} {\bar{a}}=R_0^2\, a, \quad {\bar{b}}= R_0 M_\odot b. \end{aligned}$$

In the following, we consider the cases where \({\bar{a}}=0.004\), and \({\bar{b}}=0,\pm 0.01\), respectively. The variations of the interior density and mass profiles are represented in Fig. 4. As required by physical considerations, the energy density of the stiff-fluid three-form star in the presence of the Higgs potential is a monotonically decreasing function of the radial coordinate. The variation of the energy density is practically independent on the values of the parameters of the Higgs potential and, at least for the considered range of values, it is very similar to the general relativistic case. However, there is a significant effect on the variation of the mass profile of the potential parameters, with higher mass values obtained for higher values of \({\bar{b}}\).

Fig. 3
figure 3

Variation of \(M/M_{\odot }\) as a function of the radius of the star R (km) for three-form stiff-fluid stars, for three different values of the constant \(\bar{\lambda }\) for three different values of the constant \(\bar{\lambda }\): \({\bar{\lambda }}=0\) (dashed curve), \({\bar{\lambda }}=0.02\) (dotted curve), and \(\bar{\lambda }=-0.02\) (dot-dashed curve), respectively. The solid curve represents the mass–radius relation for general relativistic stiff-fluid stars

Table 1 The maximum mass and corresponding radius for three-form field stiff-fluid stars with constant potential
Fig. 4
figure 4

Variations of \({\bar{\rho }}_m\) (left panel), and of \({\bar{m}}\) (right panel) as a function of \(\eta \) for three-form stiff-fluid stars in the presence of a Higgs-type potential, for \({\bar{a}}=0.003\) and for three different values of the constant \({\bar{b}}\): \({\bar{b}}=0\) (dashed curve), \({\bar{b}}=-0.033\) (dotted curve), and \({\bar{b}}=0.033\) (dot-dashed curve). For the central density of the star we adopt the value \(\rho _{{\text {mc}}}=3.7\times 10^{14} \mathrm{g/cm}^3\). The solid curve represents the standard general relativistic density and mass profiles for stiff-fluid stars

Fig. 5
figure 5

Variations of \({\bar{\zeta }}\) (left panel) and of \({\bar{V}}\left( {\bar{\zeta }}\right) \) (right panel) as functions of \(\eta \) for stiff-fluid three-form stars in the presence of the Higgs potential, for \({\bar{a}} =0.003\) and three different values of the Higgs potential parameter \({\bar{b}}\): \({\bar{b}}=0\) (dashed curve), \({\bar{b}}=-0.033\) (dotted curve), and \({\bar{b}}=0.033\) (dot-dashed curve). For the central density of the star we adopt the value \(\rho _{{\text {mc}}}=3.7\times 10^{14} \mathrm{g/cm}^3\), while \({\bar{\zeta }}_0=0\) and \({\bar{\zeta }}^\prime _0=7\times 10^{-3}\), respectively

Similarly to the constant potential case, the three-form field component is a monotonically increasing function inside the star, reaching its maximum value near the star surface (see Fig. 5 for details). There is a significant dependence of \({\bar{\zeta }}\) on the parameters of the Higgs potential near the vacuum boundary. The Higgs-type potential is a monotonically decreasing function of \(\eta =r/R_0\), as is transparent from Fig. 5 and, at least for the considered set of parameters, it has only negative values, reaching its minimum on the star surface. The variation of the effective energy density of the three-form field is represented in Fig. 6. As one can see from the figure, \({\bar{\rho }}_A\) is a monotonically increasing function of \(\eta \) inside the star, reaching its maximal value at the vacuum boundary.

The mass–radius relation for three-form field stiff-fluid stars in the presence of a Higgs-type potential is shown in Fig. 7. As one can see from Fig. 7, three-form field compact objects obeying the stiff-fluid equation of state can have much higher masses than their general relativistic counterparts, with the maximum mass reaching values as high as \(5.10M_{\odot }\). Generally, for all considered numerical values the stiff-fluid three-form stars in the presence of a Higgs potential have higher masses as compared to the standard general relativistic approach. The general shape of the mass–radius dependence curve is similar to the one in general relativity, with a displacement of the curve towards higher mass values.

Exact numerical values of the maximum mass and of the corresponding radius for stiff-fluid stars in the presence of a Higgs-type potential are given in Table 2.

4.2 Photon stars

The radiation fluid equation of state \(P=\rho /3\) plays a critical role in physics and astrophysics. By using this equation of state we can describe the physical properties of the cores of neutron stars, assumed to consist of cold, non-interacting and degenerate fermions [55, 56]. The interesting possibility of the existence of stars described by the radiation equation of state, and therefore consisting of a radiation fluid, was also analyzed from different perspectives [57,58,59,60]. Einstein’s field equation describing static, spherically symmetric stars made of a radiation fluid, were investigated numerically in [58]. The way the thermodynamical parameters (entropy, temperature, baryon number, mass–energy, etc.) scale with the size of a photon fluid star were investigated in [60], and an unusual behavior due to the non-extensivity of the system was found. These scaling laws have some similarities with the area scaling law of the black hole entropy.

Another interesting class of stellar-type objects described by the radiation fluid equation of state are represented by the so called “Radiation Pressure Supported Stars” (RPSS), which have the intriguing property that they can exist even in classical Newtonian gravity [61]. They have been generalized in order to incorporate the effects of standard general relativity. The corresponding classes of stellar-type objects are called “Relativistic Radiation Pressure Supported Stars” (RRPSS) [61]. On the other hand, it was already suggested in [62] that the formation of RRPSSs may occur during the gravitational collapse of extremely massive baryonic matter clouds, a process that ends in a very high density state. It turns out that independently of the details of the collapse event, at sufficiently large cosmological redshifts \(z\gg 1\), the radiation flux of the collapsed object (star or black hole) always reaches the maximal Eddington luminosity. The properties of the radiation fluid stars in a non-minimally coupled gravity model of the form \(Y (R)F^2\), where \(F^2\) is the Maxwell invariant and Y(R) is a function of the Ricci scalar R, were considered in [63].

Fig. 6
figure 6

Variations of \({\bar{\rho }}_A\) as a function of \(\eta \) for stiff-fluid three-form stars in the presence of the Higgs potential, for \({\bar{a}} =0.003\) and three different values of the Higgs potential parameter \({\bar{b}}\): \({\bar{b}}=0\) (dashed curve), \({\bar{b}}=-0.033\) (dotted curve), and \({\bar{b}}=0.033\) (dot-dashed curve). For the central density of the star we adopt the value \(\rho _{{\text {mc}}}=3.7\times 10^{14} \mathrm{g/cm}^3\), while \({\bar{\zeta }}_0=0\) and \({\bar{\zeta }}^\prime _0=7\times 10^{-3}\), respectively

Fig. 7
figure 7

Variation of the mass \(M/M_{\odot }\) as a function of the radius of the star R for a three-form stiff-fluid star in the presence of a Higgs-type potential, for \({\bar{a}}=0.003\) and three different values of the constant \({\bar{b}}\): \({\bar{b}}=0\) (dashed curve), \({\bar{b}}=-0.033\) (dotted curve), and \({\bar{b}}=0.033\) (dot-dashed curve). The solid curve represents the mass–radius relation for general relativistic stiff-fluid stars

Table 2 The maximum mass and corresponding radius for three-form field stiff-fluid stars in the presence of a Higgs-type potential for \({\bar{a}}=0.003\), and different values of \({\bar{b}}\)

In the following we will consider the properties of radiation fluid (photon) stars in three-form field gravity. We will assume again an isotropic pressure mater distribution, with \(p_{\perp m}=p_{rm}\), and we consider that the equation of state of the baryonic matter inside the star can be approximated by the photon gas equation of state

$$\begin{aligned} p_{rm}=\frac{\rho _m}{3}. \end{aligned}$$
(32)

We will solve numerically the gravitational field equations for central densities that vary in the range \({3.1\times 10^{14} \mathrm{g/cm}^3}\) and \({2.9\times 10^{15} \mathrm{g/cm}^3}\), respectively. We stop the integration when the density reaches the surface value of \(\rho _{{\text {mc}}}= 2\times 10^{14} \mathrm{g/cm}^3\) (the nuclear density). In the following, to obtain the interior solutions we consider two different choices of the three-form field potential.

Fig. 8
figure 8

Variation of the energy density \({\bar{\rho }}_{m}\) (left panel) and of the mass \({\bar{m}}\) (right panel) in terms of \(\eta \) for three-form field photon stars in the presence of a constant field potential, for three different values of the constant \(\bar{\lambda }\): \({\bar{\lambda }}=0\) (dashed curve), \({\bar{\lambda }}=0.02\) (dotted curve), and \(\bar{\lambda }=-0.02\) (dot-dashed curve), respectively. The central density of the photon star is taken as \(\rho _{{\text {mc}}}=6.4\times 10^{14} \mathrm{g/cm}^3\). The solid curve represents the standard general relativistic density and mass profiles for photon stars

Fig. 9
figure 9

Variation of \({\bar{\zeta }}\)(left panel) and \({\bar{\rho }}_A\) (right panel) as a function of \(\eta \) in three-form field photon stars, for three different values of the constant \({\bar{\lambda }}\): \({\bar{\lambda }}=0\) (dashed curve), \({\bar{\lambda }}=0.02\) (dotted curve), and \(\bar{\lambda }=-0.02\) (dot-dashed curve), respectively. The central density of the photon star is taken as \(\rho _{{\text {mc}}}=6.4\times 10^{14} \mathrm{g/cm}^3\)

4.2.1 Constant potential: \(V(A^2)=\lambda \)

As a first case in the investigation of the photon stars in three-form field gravity we assume that the field potential is a constant, \(V(A^2)=\lambda =\mathrm{constant}\). To numerically integrate the gravitational field equations we consider three different values of \({\bar{\lambda }}\), \({\bar{\lambda }}=0,\pm 0.02\). For the central density of the photon star we adopt the value \(\rho _{{\text {mc}}}=6.4\times 10^{14} \mathrm{g/cm}^3\).

The behaviors of the interior profiles of the matter energy density and of the mass are presented in Fig. 8. The matter energy density is a monotonically decreasing function of the dimensionless radial coordinate, and it reaches the nuclear density at much smaller radial coordinate values, as compared to the stiff-fluid case. The matter energy density is also much sensitive to the variations of the potential. The interior mass profile of the star is a monotonically increasing function of \(\eta \), and it strongly depends on the numerical values of \({\bar{\lambda }}\), with the mass of the star increasing for positive values, and decreasing for negative values of \({\bar{\lambda }}\), respectively.

The variations of the three-form field \({\bar{\zeta }}\) and of the effective energy density of the three-form field \({\bar{\rho }}_A\) in terms of the distance from the center to the surface of the star are shown, for different values of the constant field potential, in Fig. 9. Similarly to the stiff-fluid case, the three-form field is a monotonically increasing function inside the high density object. Near the surface of the star the values of \({\bar{\zeta }}\) are dependent on the constant values of the field potential. On the other hand, the energy density of the three-form field takes constant values inside the star, and its numerical values are strongly dependent on the model parameters.

The mass–radius relation for photon stars is plotted in Fig. 10. The presence of the three-form field inside the star generates a complicated pattern of behaviors for photon stars. Depending on the sign and numerical value of the constant potential, both higher- and lower-mass stars do exist in this model. If in standard general relativity the maximum mass of the photon star is close to \(2M_{\odot }\), three-form field photon stars can reach maximum masses of around \(2.5M_{\odot }\). For negative three-form field potentials smaller masses than in general relativity are also possible. The general shape of the mass–radius function for three-form field photon stars is similar to the standard general relativistic one.

Several numerical values of the maximum mass of photon stars in three-form field gravity are presented in Table 3. In general relativity for photon stars we have \(M_{\max }=2.03\,M_{\odot }\), \(R=12.71\,\mathrm{km} \), values corresponding to a central density \(\rho _{{\text {mc}}}=1.87\times 10^{15}\, \mathrm{g/cm}^3\).

4.2.2 Higgs-type potential: \(V(A^2)=a A^2 +b A^4\)

We consider now the case in which the three-form field photon star also contains a Higgs-type potential, of the form \(V(A^2)=a A^2 +b A^4\). For the parameters of the potential we adopt the numerical values \({\bar{a}}=0.007\) and \({\bar{b}}=0,\pm 0.11\), respectively, that is, in the following numerical investigations we keep \({\bar{a}}\) as constant, and we vary \({\bar{b}}\). Moreover, we fix the central density of the three-form field photon star as \(\rho _{{\text {mc}}}=7.12\times 10^{14}\; \mathrm{g/cm}^3\).

Fig. 10
figure 10

Variation of the total mass \(M/M_{\odot }\) as a function of the radius of the star R for three-form field photon stars, for three different values of the constant \({\bar{\lambda }}\): \({\bar{\lambda }}=0\) (dashed curve), \({\bar{\lambda }}=0.02\) (dotted curve), and \(\bar{\lambda }=-0.02\) (dot-dashed curve), respectively. The solid curve represents the mass–radius relation for general relativistic photon stars

Table 3 The maximum mass and the corresponding radius for three-form field photon stars with constant potential

The interior energy density and mass profiles of the three-form field photon stars are depicted in Fig. 11. The energy density is almost insensitive to the variations of the parameters of the Higgs potential. On the other hand, an explicit dependence on the values of \({\bar{b}}\) exist in the case of the interior mass distribution, which indicates a significant dependence of the mass of the star on the three-form field potential. We refer the reader to Fig.  12 for the behavior of the three-form field component, which is a monotonically increasing function inside the star, reaching its maximum value near the star surface. Note that the Higgs-type potential is a monotonically decreasing function of \(\eta =r/R_0\).

The variation of the energy density of the three-form field is presented in Fig. 13, indicating that \({\bar{\rho }}_A\) is a monotonically increasing function inside the star.

The mass–radius relation for three-form field photon stars in the presence of a Higgs-type potential is presented in Fig. 14. As one can see from the \(M=M(R)\) relation of three-form field photon stars in the presence of a Higgs field potential, there is a significant increase of the masses of compact objects that essentially depends on the parameters of the potential. If for standard general relativistic photon stars the maximum mass is of the order of \(1.9M_{\odot }\), the presence of the Higgs potential leads to an increase of the maximum mass to values of the order of \(3M_{\odot }\). There is also a slight increase in the radii of the maximum mass stars, as compared to the standard general relativistic case, which are of the order of 14 km.

Fig. 11
figure 11

Variation of the energy density \({\bar{\rho }}_m\) (left panel) and of the interior mass profile \({\bar{m}}\) (right panel) for three-form field photon stars in the presence of a Higgs-type potential as functions of \(\eta \), for \(\bar{a}=0.007\), and for three different values of the constant \({\bar{b}}\): \({\bar{b}}=0\) (dashed curve), \({\bar{b}}=-0.11\) (dotted curve), and \({\bar{b}}=0.11\) (dot-dashed curve). The central density of the photon star is \(\rho _{{\text {mc}}}=7.12\times 10^{14}\; \mathrm{g/cm}^3\). The solid curve represents the standard general relativistic density and mass profiles for photon stars

Fig. 12
figure 12

Variation of \({\bar{\zeta }}\) (left panel) and of the Higgs-type potential \({\bar{V}}\left( \zeta \right) \) (right panel) for three-form photon stars as functions of \(\eta \), for \({\bar{a}}=0.007\) and three different values of the constant \({\bar{b}}\), and for three different values of the constant \({\bar{b}}\): \({\bar{b}}=0\) (dashed curve), \({\bar{b}}=-0.11\) (dotted curve), and \({\bar{b}}=0.11\) (dot-dashed curve), respectively. The central density of the photon star is \(\rho _{{\text {mc}}}=7.12\times 10^{14}\; \mathrm{g/cm}^3\), while \({\bar{\zeta }}_0=0\) and \({\bar{\zeta }}^\prime _0=7\times 10^{-3}\), respectively

A selected set of numerical values of the maximum masses of three-form field photon stars are presented, for different central densities, in Table 4.

4.3 Quark stars

According to the present day view of astrophysics, in neutron stars nuclear matter, assumed to be in beta equilibrium, is present from very low densities to several times the nuclear saturation density, given by \(\rho _n=0.16\;\mathrm{fm}^{-3}\) [53, 62, 64]. However, neutrons are not really elementary particles, and they can be regarded as bound states of their valence quarks and antiquarks. Quarks are spin half particles, and thus they are fermions. In the simple phenomenological MIT bag model, the quarks inside the hadrons are confined by a bag, with the quarks constrained to move freely and independently inside the hadrons by an infinite potential well. In a more general model, called the potential model, quarks are confined inside the bag by a phenomenological confinement potential, which usually is taken as a harmonic oscillator potential. At the high nuclear densities reached in the interiors of compact neutron stars, a hadron-quark phase transition may occur, leading to the transformation of the neutron star into a quark star [65,66,67]. The equation of state of the quark matter can be derived from the fundamental Lagrangian of Quantum Chromodynamics (QCD) [68, 69],

$$\begin{aligned} L_{QCD}= & {} \frac{1}{4}\sum _{a}F_{\mu \nu }^{a}F^{a\mu \nu } \nonumber \\&+\sum _{f=1}^{N_{f}}{\bar{\psi }}\left( i\gamma ^{\mu }\partial _{\mu }-g\gamma ^{\mu }A_{\mu }^{a}\frac{\lambda ^{a}}{2}-m_{f}\right) \psi , \end{aligned}$$
(33)

where the subscript f indicates the differen(left panel) and of the Higgs-type potential \({\bar{V}}\left( \zeta \right) \) (right panel)t quark flavors u, d, s, c etc., while the nonlinear gluon field strength is defined according to the Yang–Mills prescription as \( F_{\mu \nu }^{a}=\partial _{\mu }A_{\nu }^{a}-\partial _{\nu }A_{\mu }^{a}+gf_{abc}A_{\mu }^{b}A_{\nu }^{c} \).

Fig. 13
figure 13

Variation of \({\bar{\rho }}_A\) for three-form field photon stars as function of \(\eta \), for \({\bar{a}}=0.007\) and three different values of the constant \({\bar{b}}\), and for three different values of the constant \({\bar{b}}\): \({\bar{b}}=0\) (dashed curve), \({\bar{b}}=-0.11\) (dotted curve), and \({\bar{b}}=0.11\) (dot-dashed curve), respectively. The central density of the photon star is \(\rho _{{\text {mc}}}=7.12\times 10^{14}\; \mathrm{g/cm}^3\), while \({\bar{\zeta }}_0=0\) and \({\bar{\zeta }}^\prime _0=7\times 10^{-3}\), respectively

Fig. 14
figure 14

Variation of the stellar mass \(M/M_{\odot }\) as a function of the radius of the star R for the three-form field photon stars in the presence of a Higgs-type potential, for \({\bar{a}}=0.007\) and three different values of the constant \({\bar{b}}\): \({\bar{b}}=0\) (dashed curve), \({\bar{b}}=-0.11\) (dotted curve), and \({\bar{b}}=0.11\) (dot-dashed curve), respectively. The solid curve shows the mass–radius relation for photon stars in standard general relativity

Table 4 The maximum masses and the corresponding radii for three-form field photon stars in the presence of a Higgs-type potential, for \({\bar{a}}=0.007\), and for different values of \({\bar{b}}\)

An important prediction of QCD is the weakening at short distances of the quark–quark interaction. Under the assumption that interactions of the quarks and gluons are weak, for a quark–gluon plasma at temperature T and having the chemical potential \(\mu _{f}\), one can obtain the energy density \(\rho \) and the pressure p by using thermal quantum field theory [68]. By neglecting the quark masses, in the first order perturbation theory, the equation of state of the quark-gluon plasma is given by [62, 68]

$$\begin{aligned} \varepsilon =\sum _{i=u,d,s,c;e^{-},\mu ^{-}}\varepsilon _{i}+B,p+B=\sum _{i=u,d,s,c;e^{-}, \mu ^{-}}p_{i}, \end{aligned}$$
(34)

where by B (the bag constant) we have denoted the difference between the energy densities of the non-perturbative and perturbative QCD vacuum. Hence, in the MIT bag model the equation of state for quark matter can be obtained as [62, 68, 69]

$$\begin{aligned} p_{rm}=\frac{1}{3}\left( \rho _{m} -4B\right) . \end{aligned}$$
(35)
Fig. 15
figure 15

Variation of the energy density \({\bar{\rho }}_{m}\) (left panel) and of the mass \({\bar{m}}\) (right panel) in terms of \(\eta \) for three-form field quark stars in the presence of a constant field potential, for three different values of the constant \(\bar{\lambda }\): \({\bar{\lambda }}=0\) (dashed curve), \({\bar{\lambda }}=0.02\) (dotted curve), and \(\bar{\lambda }=-0.02\) (dot-dashed curve), respectively. The central density of the quark star is taken as \(\rho _{{\text {mc}}}=1.23\times 10^{15} \mathrm{g/cm}^3\). The solid curve represents the standard general relativistic density and mass profiles for quark stars. The horizontal line in the left panel shows the surface density \(\left. {\bar{\rho }} _m\right| _{\mathrm{surface}}=4B\) of the star

Fig. 16
figure 16

Variation of \({\bar{\zeta }}\)(left panel) and \({\bar{\rho }}_A\) (right panel) as a function of \(\eta \) in three-form field quark stars, for three different values of the constant \({\bar{\lambda }}\): \({\bar{\lambda }}=0\) (dashed curve), \({\bar{\lambda }}=0.02\) (dotted curve), and \(\bar{\lambda }=-0.02\) (dot-dashed curve), respectively. The central density of the quark star is taken as \(\rho _{{\text {mc}}}=1.23\times 10^{15} \mathrm{g/cm}^3\)

Equation (35) corresponds to the equation of state of a system of massless particles, containing the important corrections originating from the QCD trace anomaly, and from perturbative interactions in the system. All these effects are described by the bag constant B. The corrections to the equation of state of the free photon gas are always negative, and they have the important consequence of reducing, at a given temperature, the energy density of the quark-gluon plasma by a factor of about two [68]. Quark stars have the significant property that their surface density is non-zero, reaching the value \(\rho _m=4B\). On the other hand, the quark pressure vanishes on the quark star’s surface, thus allowing the definition of the radius of the compact object.

In the following numerical simulations of the structure of isotropic quark stars, with \(p_{rm}=p_{\perp m}\), in three-form field gravity we have adopted for the bag constant the numerical value \(B=10^{14}\;\mathrm{g/cm}^3\). We consider central densities in the range \(4.2\times 10^{14}\;\mathrm{g/cm}^3\) and \(2.2\times 10^{16}\;\mathrm{g/cm}^3\), respectively. In all cases, the numerical integration stops at \(p_m=0\), which corresponds \(\rho _m=4 B\), according to the equation of state Eq. (35).

4.3.1 Constant potential case: \(V(A^2)=\lambda \)

We investigate first quark stars in the presence of a constant three-form field potentia, with \(V(A^2)=\lambda \). We consider three different values of \({\bar{\lambda }}\) as \({\bar{\lambda }}=0,~\pm 0.02\). The density and mass profiles for these cases are shown in Fig. 15. The density is a monotonically decreasing function of the radial distance, and on the surface of the star it reaches the value 4B, while the pressure identically vanishes. The density variation is not influenced significantly by the presence of the three-form field, and it is similar to the general relativistic case. The mass profile is monotonically increasing, and indicates that the field potential can have a major influence on the mass of the quark stars.

The variations of the three-form field and of \({\bar{\rho }}_A\) are depicted in Fig. 16. Similarly to the previous stellar models, the three-form field is a monotonically increasing function of the radial coordinate \(\eta \), reaching its maximum value at the star surface. The variation of the field is slightly dependent on the values of the potential, and this dependence appears mostly near the surface of the star. Similarly to the other constant potential cases, the effective energy density of the three-form field is constant inside the quark star.

The mass–radius relation for quark stars are shown in Fig. 17. There is a significant effect of the three-form field, and of its potential, on the maximum masses, and stability regions of quark stars. Increasing the value of the constant field potential leads to a large increase in the maximum mass. However, quark stars with masses smaller than the general relativistic ones are also possible. The general form of the mass–radius relation is also preserved.

The maximum masses of a selected sample of three-form field quark stars are presented, for different central densities, in Table 5.

In general relativity for quark stars we have \({M_{\max }=2.03\,M_{\odot }}\), \({R=11.07\,\mathrm{km}}\), values corresponding to a central density of \({\rho _{{\text {mc}}}=1.93\times 10^{15}\, \mathrm{g/cm}^3}\).

4.3.2 Higgs-type potential: \(V(A^2)=a A^2 +b A^4\)

In the case of three-form field quark stars in the presence of a Higgs-type potential for the model parameters we adopt the values \({\bar{a}}=0.01\), and we consider three different values of \({\bar{b}}\), as \({\bar{b}}=0,~\pm 0.11\). The density and mass profiles for these forms of the potential are shown in Fig. 18. The density monotonically decreases from the center towards the surface of the star, where the pressure vanishes. The radius of the star is not significantly influenced by the parameters of the Higgs potential. However, a significant influence does appear in the case of the mass distribution, with the mass reaching much higher values than in the general relativistic case.

Table 5 The maximum mass and the corresponding radius for three-form field quark stars with constant potential
Fig. 17
figure 17

Variation of the total mass \(M/M_{\odot }\) as a function of the radius of the star R for three-form field quark stars, for three different values of the constant \({\bar{\lambda }}\): \({\bar{\lambda }}=0\) (dashed curve), \({\bar{\lambda }}=0.02\) (dotted curve), and \(\bar{\lambda }=-0.02\) (dot-dashed curve), respectively. The solid curve represents the mass–radius relation for general relativistic quark stars

The variation of the three-form field and its Higgs-type potential are depicted in Fig. 19. The three-form field \({\bar{\zeta }}\) is a monotonically increasing function inside the quark star, reaching its maximum at the surface. Also near the surface the dependence of the field on the parameters of the potential becomes stronger. On the other hand, the Higgs-type potential, is a monotonically decreasing curve inside the star, and it takes negative values in all range of distances from the center to the surface. When approaching the surface the variation of the potential shows a significant dependence on the parameters \({\bar{a}}\) and \({\bar{b}}\).

The variation of \({\bar{\rho }}_A\) for three-form quarks stars in the presence of a Higgs-type potential is shown in Fig. 20, indicating an increase of the field energy density inside the quark star.

The mass–radius relation for three-form field quark stars in the presence of a Higgs-type potential are shown in Fig. 21. There is also a significant increase of the masses of the quark stars in three-form field gravity, as compared to the general relativistic case. The radii corresponding to the maximum masses also do increase.

Fig. 18
figure 18

Variation of the energy density \({\bar{\rho }}_m\) (left panel) and of the interior mass profile \({\bar{m}}\) (right panel) for three-form field quark stars in the presence of a Higgs-type potential as functions of \(\eta \), for \(\bar{a}=0.01\), and for three different values of the constant \({\bar{b}}\): \({\bar{b}}=0\) (dashed curve), \({\bar{b}}=-0.11\) (dotted curve), and \({\bar{b}}=0.11\) (dot-dashed curve). The central density of the quark star is \(\rho _{{\text {mc}}}=1.28\times 10^{15}\; \mathrm{g/cm}^3\). The solid curves correspond to the standard general relativistic case. The horizontal line in the left panel shows the surface density \(\left. {\bar{\rho }} _m\right| _{\mathrm{surface}}=4B\) of the star

Fig. 19
figure 19

Variation of \({\bar{\zeta }}\) (left panel) and of the Higgs-type potential \({\bar{V}}\left( \zeta \right) \) (right panel) for three-form quark stars as functions of \(\eta \), for \({\bar{a}}=0.01\) and for three different values of the constant \({\bar{b}}\): \({\bar{b}}=0\) (dashed curve), \({\bar{b}}=-0.11\) (dotted curve), and \({\bar{b}}=0.11\) (dot-dashed curve), respectively. The central density of the photon star is \(\rho _{{\text {mc}}}=1.28\times 10^{15}\; \mathrm{g/cm}^3\), while \({\bar{\zeta }}_0=0\) and \({\bar{\zeta }}^\prime _0=7\times 10^{-3}\), respectively

Some specific values of the maximum masses of the three-form field quark stars in the presence of the Higgs-type potential are presented in Table 6.

Table 6 The maximum masses and the corresponding radii for three-form field quark stars in the presence of a Higgs-type potential for \({\bar{a}}=0.01\) and different values of \({\bar{b}}\)

4.4 Bose–Einstein condensate stars

If the temperature T of the boson gas drops below the critical value \(T_{cr}\), given by \({T_{cr}=2\pi \hbar ^2\rho _{cr}^{2/3}/ \zeta ^{2/3}(3/2)m^{5/3}k_B}\), where m is the particle mass in the condensate, \(\rho _{cr}\) is the critical transition density, \(k_B\) is Boltzmann’s constant, and \(\zeta \) is the Riemann zeta function, respectively, then the gas undergoes a phase transition becoming a Bose–Einstein condensate, with all particles in the system quantum mechanically correlated [70]. It is assumed that due to the superfluid properties of the neutron matter a significant amount of mass of astrophysical compact objects may be in the form of a Bose–Einstein condensate [13, 62]. In the non-relativistic and the Newtonian theoretical limits a Bose–Einstein condensate confined gravitationally can be described as a gas, with the pressure and density related by a barotropic equation of state of the form \(p=p(\rho )\). The equation of state of the condensate is characterized by two important physical parameters, the mass of the condensate particle m, and the scattering length a [71]. For a condensate with quartic non-linearity, the equation of state is polytropic with index \(n = 1\), and it is given by \(p (\rho )=K\rho ^2\) [71,72,73]. In the case of a neutron condensate fluid, in which the nuclear matter underwent a phase transition after the neutrons formed Cooper pairs, the constant K is given by [13]

$$\begin{aligned} K=\frac{2\pi \hbar ^{2}a}{m_c^{3}} =0.1856\times 10^5 \left( \frac{a}{1\;\mathrm {fm}}\right) \left( \frac{m_c}{2m_n}\right) ^{-3}, \end{aligned}$$
(36)

where \(m_n=1.6749\times 10^{-24}\) g denotes the mass of the neutron. If the particles at the core of compact high density stellar objects form Cooper pairs having masses of the order of two neutron masses, and by assuming a scattering length of the order of 10–20 fm, Bose–Einstein condensate stars can achieve maximum central densities of the order of 0.1–0.3 \(\times 10^{16}\) g/cm\(^3\), minimum radii in the range of 10–20 km, and maximum masses of the order of 2\(M_{\odot }\) [13].

Fig. 20
figure 20

Variation of \({\bar{\rho }}_A\) for three-form quark stars as functions of \(\eta \), for \({\bar{a}} =0.01\) and for three different values of the constant \({\bar{b}}\): \({\bar{b}}=0\) (dashed curve), \({\bar{b}}=-0.11\) (dotted curve), and \({\bar{b}}=0.11\) (dot-dashed curve), respectively. The central density of the photon star is \(\rho _{{\text {mc}}}=1.28\times 10^{15}\; \mathrm{g/cm}^3\), while \({\bar{\zeta }}_0=0\) and \({\bar{\zeta }}^\prime _0=7\times 10^{-3}\), respectively

Fig. 21
figure 21

Variation of the stellar mass \(M/M_{\odot }\) as a function of the radius of the star R for the three-form field quark stars in the presence of a Higgs-type potential, for \({\bar{a}}=0.01\) and three different values of the constant \({\bar{b}}\): \({\bar{b}}=0\) (dashed curve), \({\bar{b}}=-0.11\) (dotted curve), and \({\bar{b}}=0.11\) (dot-dashed curve), respectively. The solid curve shows the mass–radius relation for quark stars in standard general relativity

Fig. 22
figure 22

Variations of the density profiles \({\bar{\rho }}_m \) of the superfluid matter (left panel) and of the normalized mass profile \({\bar{m}}\) (right panel) as a function of the radial distance \(\eta \) for the three-form field Bose–Einstein condensate stars, for three different values of the constant \({\bar{\lambda }}\): \(\bar{\lambda }=0\) (dashed curve), \({\bar{\lambda }}=0.0013\) (dotted curve), and \({\bar{\lambda }}=-0.01\) (dot-dashed curve). The solid curves represent the corresponding standard general relativistic quantities

Fig. 23
figure 23

Variation of \({\bar{\xi }}\) (left panel) and \({\bar{\rho }}_A\) (right panel) as a function of distance from the center of the star \(\eta \) for the three-form field Bose–Einstein condensate stars, for three different values of the constant field potential \({\bar{\lambda }}\): \(\bar{\lambda }=0\) (dashed curve), \({\bar{\lambda }}=0.0013\) (dotted curve), and \({\bar{\lambda }}=-0.01\) (dot-dashed curve)

In the following we will investigate the properties of Bose–Einstein condensate/superfluid stars in three-form field gravity theory. We assume that matter inside the star can be described by the equation of state

$$\begin{aligned} p_{rm}=K \rho _m^\gamma , \end{aligned}$$
(37)

where \(\gamma \) is a constant. We consider again an isotropic pressure distribution with \(p_{\perp m}=p_{rm}\). We consider that the superfluid interior of the star can be described by a Bose–Einstein condensate, and hence we adopt for \(\gamma \) the value \(\gamma =2\). As for the constant K in the following we fix it to the rescaled value \({\bar{K}}=\epsilon _0 K=0.4\). We will numerically investigate three-form field Bose–Einstein condensate stars for two different choices of the field potential. In our analysis we assume that the central density varies between the values \(2.44\times 10^{14}\; \mathrm{g/cm}^3\) and \(6.43\times 10^{15}\; \mathrm{g/cm}^3\), respectively.

4.4.1 Constant potential: \(V(A^2)=\lambda \)

As a first example of a Bose–Einstein condensate star in three-form field gravity we consider the case of the constant field potential, \(V(A^2)=\lambda =\mathrm{constant}\). We will adopt three different choices for \({\bar{\lambda }}\) as \({\bar{\lambda }}=0, \,0.0013\) and \(-0.01\), respectively. The variation of the matter density and of the mass profile inside the three-form field Bose–Einstein condensate star are represented in Fig. 22.

As one can see from the plots, the density profile is a monotonically decreasing function of the radial coordinate, while the mass profile is a monotonically increasing function of \(\eta =r/R_0\). The density and the pressure vanish on the vacuum boundary of the star, and hence the radius of the star is uniquely defined. Both the density and the mass profiles show a significant dependence on the numerical values of the potential, and a significant increase in the numerical values of the mass of the condensate star can be obtained by varying the values of \({\bar{\lambda }}\).

The variations of the non-zero component of the three-form field and of the energy density \({\bar{\rho }}_A\) in the interior of the star are depicted in Fig. 23. The three-form field is a monotonically increasing function of the radial distance from the center of the star, and it reaches its maximum on the star surface. For the adopted set of potential values the variation of \({\bar{\zeta }}\) is almost independent on \({\bar{\lambda }}\). The energy density of the field remains a constant inside the star.

The mass–radius relation of the three-form field Bose–Einstein condensate stars is depicted in Fig. 24. The three-form field Bose–Einstein condensate stars with a constant field potential may have higher masses than similar general relativistic stars. The mass increase/decrease depends on the sign of the constant potential. The increase of the mass is also associated to a slight increase of the radius of the star. A selected number of numerical values of the maximum masses of the three-form field Bose–Einstein condensate stars in the presence of a constant potential are presented in Table 7.

In general relativity for Bose–Einstein condensate stars we have \(M_{\max }=2.00\,M_{\odot }\), \(R=11.15\,\mathrm{km} \), values corresponding to a central density \({\rho _{{\text {mc}}}=2.58\times 10^{15}\, \mathrm{g/cm}^3}\).

4.4.2 Higgs-type potential: \(V(A^2)=a A^2 +b A^4\)

We consider now three-form field Bose–Einstein condensate stars in the presence of a Higgs-type potential \(V(A^2)=a A^2 +b A^4\). For our numerical investigations we fix the value of \({\bar{a}}\) as \({\bar{a}}=0.001\), and we consider three different choices of \({\bar{b}}\) as \({\bar{b}}=0,\, 0.6\) and \(-0.1\). In the cases \({\bar{b}}=0,\, 0.6\) the density profile has a local minimum, where its value is still greater than zero. Since the density profile should be a decreasing function of the distance from the center of the compact star, we stop integration around the minimum density \( 2.40\times 10^{13}\, \mathrm{g/cm}^3\), which is smaller than the nuclear density. However, in the case \({\bar{b}}=-0.1\), the density vanishes at a finite \(r/R_0\), and hence the radius of the star can be determined exactly.

The variations of the densities and masses of the three-form field Bose–Einstein condensate stars in the presence of a Higgs potential are plotted in Fig. 25. The density of the condensate stars is a monotonically decreasing function of \(rR_0\), but its behavior near the star surface strongly depends on the model parameters. The radius of the star can be defined either exactly, as the radial distance at which the density vanishes, or as the point corresponding to the local minimum of the density. Once the cut-off point of the density is obtained, one can also uniquely find the maximum mass of the star.

Fig. 24
figure 24

Variations of the mass \(M/M_\odot \) of the three-form field Bose–Einstein condensate stars as functions of the radius of the star R, for three different values of the constant field potential \({\bar{\lambda }}\): \(\bar{\lambda }=0\) (dashed curve), \({\bar{\lambda }}=0.0013\) (dotted curve), and \({\bar{\lambda }}=-0.01\) (dot-dashed curve). The solid curve represents the mass–radius relation for Bose–Einstein condensate stars in standard general relativity

Table 7 The maximum masses and the corresponding radii for three-form field Bose–Einstein stars with constant potential

The variations of the three-form field component \({\bar{\zeta }}\) and of the Higgs-type potential \({\bar{V}}\left( {\bar{\zeta }}\right) \) are represented in Fig. 26. While for smaller values of \({\bar{b}}\), \({\bar{\zeta }}\) is a monotonically decreasing function, reaching its maximum on the stellar surface, for higher values of \({\bar{b}}\), \({\bar{\zeta }}\) reaches its maximum below the surface, and tends to decrease towards the vacuum boundary. The dependence of the behavior of the potential on \({\bar{b}}\) is even stronger. For large values of \({\bar{b}}\) the potential becomes a positive function with the maximum beyond the star’s surface, while for smaller values it is a monotonically decreasing function of the radial coordinate, taking negative values inside the star.

Fig. 25
figure 25

Variations of the matter density profile \({\bar{\rho }}_m\) (left panel) and of the mass profile \({\bar{m}}\) as functions of \(\eta \) for the three-form field Bose–Einstein condensate stars in the presence of a Higgs-type potential, for \({\bar{a}}=0.001\) and three different values of the constant \({\bar{b}}\): \({\bar{b}}=0\) (dashed curve), \({\bar{b}}=-0.1\) (dotted curve), and \({\bar{b}}=0.6\) (dot-dashed curve), respectively. The solid represents the same quantities in standard general relativity

Fig. 26
figure 26

Variation of the non-zero component of the three-form field \({\bar{\zeta }}\) (left panel) and of the Higgs-type potential (right panel) as a function of the radial coordinate \(\eta \) for three-form field Bose–Einstein condensate stars in the presence of a Higgs-type potential for \({\bar{a}}=0.001\), and three different values of the constant \({\bar{b}}\): \({\bar{b}}=0\) (dashed curve), \({\bar{b}}=-0.1\) (dotted curve), and \({\bar{b}}=0.6\) (dot-dashed curve), respectively

The variation of \({\bar{\rho }}_A\) inside the Bose–Einstein condensate star is represented in Fig. 27. The energy density of the field has a complex behavior, which is strongly dependent on the model parameters, and it can take both positive and negative values inside the star.

The mass–radius relations for three-form field Bose–Einstein condensate stars in the presence of a Higgs-type potential are depicted in Fig. 28. Three-form field Bose–Einstein condensate stars in the presence of a Higgs potential can reach higher maximum masses as compared to the general relativistic case, which can exceed \(2M_{\odot }\). However, no significant variation in the radii of the stars does appear due to the presence of the three-form field. Specific values of the maximum masses of three-form field Bose–Einstein condensate stars in the presence of a Higgs-type potential are presented in Table 8.

5 Discussions and final remarks

In the present paper, we have investigated the possible existence of stellar-type massive compact astrophysical objects, described, together with their baryonic content, by the three-form field gravitational theory, in which the standard Hilbert–Einstein action of general relativity is non-trivially extended by the addition of the Lagrangian of a three-form field to the total action. In order to investigate the gravitational properties of this theoretical model we have considered the simplest case, corresponding to a stellar interior described by a static and spherically symmetric geometry. In this case the system of the gravitational field equations in three-form field gravity also depends on the scalar radial function \(\zeta \), the radial component of the dual vector \(B^{\delta }\) of the three-form \(A_{\alpha \beta \gamma }\), and on the arbitrary potential \(V (\zeta )\) of the three-form field. Thus, more degrees of freedom appear than in the case of standard general relativity, in which all the properties of compact objects are determined only by baryons.

As a first step in our study, we have derived the basic equations describing the structure of static spherically symmetric compact objects in three-form field gravity, namely, the mass continuity equation, the generalized hydrostatic equilibrium equation (Tolman–Oppenheimer–Volkoff equation), and the field equation describing the evolution of the three-form field \(\zeta \) in the given static geometry. An important physical parameter determining the properties of the stars in three-form field gravity is the self-interaction potential \(V(\zeta )\) of the three-form field. In the present study we have adopted only two functional forms for \(V(\zeta )\), by assuming that it is either constant inside the stars, or it is of the Higgs type, a second choice which is supported by the role such potentials play in elementary particle physics. The case of the constant potential can also be interpreted by assuming that the three-form field is in the minimum of the Higgs-type potential, and therefore \(V(\zeta )\) assumes a constant value. Indeed, more general forms of the potential (exponential, hyperbolic, trigonometric, power-law etc.) can also be considered, and we expect that the modification of the analytical form of the potential will have a significant impact on the stellar properties, leading to compact objects having different global properties. Other choices of the potential would definitely generate different stellar sequences, corresponding to lighter or perhaps even more massive stars. Unfortunately, presently there are no firm theoretical results establishing a clear relationship between the magnitude of the mass of the star and the field potential, and therefore the role of \(V(A^2)\) can only be understood from the detailed numerical analysis of different stellar models under different physical and initial conditions.

Even in the framework of the simple spherically symmetric static model the field equations of the three-form field theory become extremely complicated. In order to close the system of field equations we must specify the three-form field potential, as well as the equation of state of the baryonic matter. An alternative approach would be to specify the functional forms of the three-form field, together with the form of the baryonic matter density. However, in our present study we have adopted a very general approach to the stellar structure problem, in which we have specified only the three-form field potential (as a constant, or Higgs type), and the equation of state of the baryonic matter. We have considered four types of stellar models, corresponding to four choices of the baryonic matter equation of state.

Fig. 27
figure 27

Variation of the energy density of the three-form field \({\bar{\rho }}_A\) as a function of the radial coordinate \(\eta \) for three-form field Bose–Einstein condensate stars in the presence of a Higgs-type potential for \({\bar{a}}=0.001\), and three different values of the constant \({\bar{b}}\): \({\bar{b}}=0\) (dashed curve), \({\bar{b}}=-0.1\) (dotted curve), and \({\bar{b}}=0.6\) (dot-dashed curve), respectively

Fig. 28
figure 28

Variation of \(M/M_{\odot }\) as a function of the radius R for three-form field Bose–Einstein condensate stars in the presence of a Higgs-type potential, for \({\bar{a}} =0.001\) and three different values of the constant \({\bar{b}}\): \({\bar{b}}=0\) (dashed curve), \({\bar{b}}=-0.1\) (dotted curve), and \({\bar{b}}=0.6\) (dot-dashed curve), respectively. The solid curve represents the mass–radius relation for general relativistic Bose–Einstein condensate stars

Table 8 The maximum masses and the corresponding radii for three-form field Bose–Einstein condensate stars in the presence of a Higgs-type potential for \({\bar{a}}=0.001\), and different values of \({\bar{b}}\)

The first of these equations, the stiff-fluid equation of state [52], gives the limiting case of the causality condition, and guarantees that the speed of sound cannot exceed the speed of light in the baryonic matter. The radiation fluid equation of state plays an important role in astrophysics and physics, and it may be used to model the dense cores of neutron stars [55, 56]. Of particular interest are the quark matter [65, 66] and Bose–Einstein condensate matter [71, 72] equations of state. If the quark model of hadrons is correct, then at high densities (not uncommon inside neutron stars) a deconfinement phase transition must take place, freeing the quarks from the bag, and leading to the formation of a quark star. Quark stars may be born during the collapse of the core of a massive neutron star after the supernova explosion, which may trigger a first or second order phase transition, by the conversion of neutron matter to quark matter in the core of a neutron star, or by the accretion of matter by neutron stars in low-mass X-ray binaries [62, 69]. On the other hand we may conjecture that in three-form field gravity a hadronic matter–quark matter phase transition can take place in extreme astrophysical and gravitational conditions, with the phase transition induced by the presence of the three-form field with Higgs-type self-interaction potential. Some possible cosmic environments in which the phase transition could take place are gamma ray bursts, supernova explosions, or accretion by neutron stars. After such a phase transition the three-form field star reaches a stable state in the minimum of the Higgs potential.

Superfluidity is also assumed to be a common occurrence in high density compact objects. The neutrons inside the star can form Cooper pairs, and a phase transition to a Bose–Einstein condensate can also occur inside the star. Such a bosonic phase transition can also significantly affect the properties of compact objects [13].

By numerically integrating the structure equations of the star, for the set of four equations of state we have considered, we have constructed four classes of three-form field gravity stellar models, corresponding to the stiff-fluid, radiation fluid, quark matter and Bose–Einstein condensate superfluid phase, respectively. In all of these cases we have thoroughly investigated the astrophysical properties of the stars (density, pressure and mass distributions, three-form field behavior), and compared them to similar quantities obtained for standard general relativistic stars. In some of our numerical simulations the three-form potential attains negative values due to the specific choices of \(V(A^2)\) and the relation between \(A^2\) and \(\zeta \), given by Eq. (15). It is also important to note that during our analysis, the specific cases where the three-form potential is given by a constant, the kinetic term of the form field, i.e. \(F^2\) is also constant, rendering the three-form field energy density to behave in the same manner, that is \(\rho _A=\)cte. This is ascribed to the very well known fact [34] that a massless three-form behaves exactly as a cosmological constant. In such cases, the theory reduces to pure GR plus a cosmological constant. Our investigations indicate that for all these four equations of state the three-form field gravity stars are generally more massive than their standard general relativistic counterparts. For example, for the stiff-fluid equation of state, in the presence of the Higgs potential, three-form field stars can reach maximum masses of the order of \(5.2M_{\odot }\), while the stiff-fluid general relativistic stars may reach masses of the order of \(3.2M_{\odot }\) [53]. Three-form field stiff-fluid stars in the presence of a constant potential could have masses as high as \(4.2M_{\odot }\), which are much heavier than the general relativistic stars.

A similar trend can be detected in the case of photon stars, where in the presence of the Higgs-type potential the maximum mass of the star can reach values as high as \(3M_{\odot }\). In the case of quark stars, whose maximum mass is of the order of \(2M_{\odot }\) in standard general relativity, three-form field gravity can induce a significant increase of the mass, which, in the presence of the Higgs-type potential, can reach values as high as \(2.7M_{\odot }\). Thus, for the same central density three-form field quark stars have significantly larger masses. The mass range of the considered superfluid Bose–Einstein condensate three-form field stars is around 1.97–2.08\(M_{\odot }\), which is still higher than the standard general relativistic mass of about \(1.8M_{\odot }\). The mass of the three-form field star is also strongly influenced by its central density, with high central density stars having smaller gravitational masses.

Hence, it turns out that three-form field stars can have a very large mass spectrum. Recently, high precision determinations of the neutron star mass distribution have also confirmed the intriguing fact of the existence of neutron stars with masses in the range of 2–2.20\(M_{\odot }\) [7, 74,75,76,77,78,79]. Such a neutron star having an unexpectedly high mass is the eclipsing binary millisecond Black Widow Pulsar B1957+20, with the mass assumed to be in the range 1.6–2.4\(M_{\odot }\) [77]. The recently discovered J0740+6620 pulsar has an estimated mass of around \(2.14M_{\odot }\) [7]. And of course there is the problem of the very high mass (\(\sim 2.6M_{\odot }\)) of the assumed neutron star in the GW190814 event [5]. A range of 2–2.4\(M_{\odot }\) is very difficult to explain by the standard hadronic matter models, constructed by using the present day knowledge of the nuclear equations of state, and in the framework of general relativity. Even the consideration of exotic models like quark or kaon stars cannot lead to a better understanding of the observations.

While explaining this mass range is problematic in standard general relativity, it can be easily explained in three-form field gravity, where one can construct stellar models having these numerical values of the masses, by using the standard knowledge about the equation of state of dense matter. The large mass spectrum of the three-form field gravity stars also could point towards the possibility that high-mass objects, usually interpreted as stellar mass black holes, and having masses in the range of \(3.8M_{\odot }\) and \(6M_{\odot }\), respectively, could be in fact three-form field stars. Such a possibility has been already investigated in [9] for the case of the quarks stars in the color–flavor locked phase. On the other hand, based on the present results on the structure and properties of three-form field stars, the possibility that at least some stellar mass black holes are in fact ordinary three-form field stars in the presence of a field potential cannot be rejected a priori. As we have shown in detail, three-form field gravity stars could have much higher masses than standard neutron stars, and thus they represent some ordinary star alternatives to low-mass black hole candidates.

A basic problem in theoretical astrophysics is finding a convincing way that could differentiate between different types of stellar models, and which also put some constraints on the equation of state of the high density matter. One of these possibilities may be related to the determination of the surface redshift of different types of compact objects. The surface redshift is defined as

$$\begin{aligned} z=\left( 1-\frac{2GM}{c^2R}\right) ^{-1/2}-1. \end{aligned}$$
(38)

In standard general relativity, where the mass–radius ratio satisfies the Buchdahl inequality \(2GM/c^2R\le 8/9\) [80], the surface redshift satisfies the constraint \(z\le 2\). A selected sample of surface redshifts of the maximum mass stars of all considered three-form field models is presented in Table 9, with the corresponding standard general relativistic values also presented. As one can see from the table, there are significant differences between the general relativistic and the three-form field stars from the point of view of the redshift. Important differences also appear between different types of stars, indicating the important role the equation of state of dense matter plays in the description of the global astrophysical parameters of relativistic stars. Negative values of the constant three-form potential lead to smaller redshift values as compared to the GR case, while positive values of \(\lambda \), as well as the Higgs potential case lead to higher redshift values. An interesting case is represented by the Bose–Einstein condensate three-form field stars, having almost the same surface redshift as their general relativistic counterparts. However, distinguishing different types of stellar objects based solely on their surface redshift may prove to be an extremely difficult observational task, which would require a significant increase in the observational accuracy.

Table 9 The surface redshift of different types of three-form field stars with maximum mass, for different model parameters

Another important observational possibility that may allow one to distinguish between three-form field stars and standard general relativistic stars or stellar mass black holes could be represented by the detailed study of the physical and astrophysical properties of the thin accretion disks that form around compact objects. The radiation properties of the accretion disks around general relativistic stars and black holes and three-form field stars are different due to the differences in both internal and external geometry. Therefore the electromagnetic emission properties of the accretion disk (energy flux, temperature, luminosity), and of the compact central objects, may provide the basic signature that could lead to the convincing differentiation between modified gravity stars and black holes from their standard general relativistic counterparts [81,82,83,84,85,86]. Multimessenger astronomy and gravitational wave observations will also certainly play an important role in this direction.

Three-form field gravity stars manifest a very complex internal structure, which also determines a complex stellar dynamics. The presence of the three-form field can lead to a number of distinctive astrophysical signatures, which could help in the observational detection of these types of objects. However, this may prove to be an extremely difficult task. Further issues related to the astrophysical/observational importance of the three-form field gravity stars will be considered in a future publication.