1 Introduction

It is possible to merger geometry and matter in the same action. An enlightening discussion regarding this question was presented in [1], in which an action like

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

was proposed, with \(\kappa =8\pi G/c^{4}\), G the gravitational constant, c the speed of light, g the metric determinant, L the matter lagrangianFootnote 1 and R the Ricci scalar.

Interestingly, the dynamics in this theory can only exist in the presence of matter, which, indeed, suggests a deeper link between space-time and matter. The inexistence of dynamics in the absence of matter in a theory of gravity fulfills Einstein’s initial proposal of having a gravity theory satisfying Mach’s principle [2].

It was also shown that the theory in [1] reduces to a special case of the scalar-tensor pressuron theory [3, 4].

The theory described from action (1) can be seen as a particular case of the well-known f(RL) theory [5], proposed by T. Harko and F.S.N. Lobo. In [5], the authors generalized the f(R)-type gravity models [6,7,8] by assuming that the gravitational Lagrangian is given by an arbitrary function of both R and L.

The viability of f(RL) candidates for dark energy was analysed from a dynamical system approach in [9]. Some constraints were put to f(RL) theories using the COBE-FIRAS measurement of the spectral radiance of the cosmic microwave background [10]. Constraints on f(RL) gravity were also put via energy conditions in [11, 12]. Wormhole solutions were also investigated in the f(RL) gravity context as one can check [13, 14].

Some f(RL) models do not conserve the energy-momentum tensor and the mechanism responsible for that was said to be a gravitational induced particle production, as it was carefully discussed in [15, 16].

It is also worth to quote that in Ref. [17], it was indicated that the f(RL) theories of gravity may be regarded as a subclass of the also well known f(RT) gravity theories [18], for which T is the trace of the energy-momentum tensor. Compact stars have been intensively studied in modified theories of gravity, in particular, the f(RT) gravity has attracted a lot of researcher’s attention [19,20,21,22,23,24,25,26]. The f(RT) gravity can be considered a generalization of the well known f(R) theory of gravity. f(R) gravity is one of the most famous modified theory of gravity being much explored in past literature [27,28,29,30,31,32]. Among the motivations of using modified theories we have: the possibility of surpassing the standard maximum mass limits of compact stars and, thus, get in touch with recent observations of massive pulsars and the indirect evidence of super-Chandrasekhar white dwarfs, and, higher order curvature corrections or curvature-matter coupling may occur in the extreme gravity regimes inside a neutron star. However, f(R) and f(RT) gravity lacks the possibility of getting a very small neutron star radius, and thus, the Buchdal and Schwarzschild radius limits are not surpassed. A non-minimal GMC model is a particular choice within the f(RL) gravity. The f(RL) gravity has two major advances in comparison with the other quoted modified theories of gravity: (1) The energy-momentum tensor of the f(RL) gravity is covariantly conserved by using natural choices for the matter lagrangian and for any choice of the f(RL) functional, this is frequently a problem within the f(RT) theory of gravity, (2) the junction with the exterior Schwarzschild solution is easily satisfied within the GMC model in f(RL) gravity, which in turn is often an issue in f(R) models.

Moreover, in Reference [33], the f(RL) gravity action was generalized by inserting on it a scalar field and a kinetic term, constructed from the gradients of the scalar field. A further model with geometry-matter coupling (from now on referred to as GMC) was proposed by Harko in [34]. For a review on generalized GMC theories, one can check [35].

GMC models have shown to be able to provide great outcomes when applied to fundamental issues of standard gravity, such as dark matter and dark energy, as one can check Refs. [36, 37]. The first GMC model was proposed by Nojiri and Odintsov in [38], where they apply a GMC theory to explain the cosmic acceleration. A GMC model was also considered for solar-like stars in [39, 40], where authors used the Newtonian limit to obtain numerical results for hydrostatic equilibrium of this type of stars and also to derive constraints to the GMC model.

Here in this work, instead, we will be concerned to the outcomes of applying a GMC model for obtaining the hydrostatic equilibrium configurations of neutron stars (NSs). NSs are supernova remnants known for their high density, strong gravitational field and rapid rotation rate. NS binary systems are among the leading gravitational wave sources [41, 42] and have, indeed, been already detected by Advanced LIGO and Advanced Virgo detectors [43].

Although most NSs have masses \(\sim 1.3-1.4 {\hbox {M}}_\odot \) [44], there is ample observational support for NSs with greater masses (\(\sim 2M_\odot \)) [45, 46], which causes some controversy regarding the NSs origin, equation of state (EoS) or even the underlying gravitational theory.

It is worth to remark that in a recent article [47], it was concluded that the GW 170817 event, coming from a binary system with \(2.74^{+0.04}_{-0.01}M_\odot \), may have resulted in a super-massive magnetar. Also, in [48] a pulsar with \(2.27^{+0.17}_{-0.15}M_\odot \) was reported, being the most massive NS already detected, named PSR J2215 + 5135. If NSs in f(RL) gravity can attain this value for the mass it will be the first time that such an object is theoretically predicted in alternative gravity. It also worth to cite that hydrostatic equilibrium configurations of compact were never performed before for the GMC theory considered in this work.

In this work we will investigate the possibility of predicting the above pulsar by altering the underlying gravitational theory. Particularly, we will investigate the hydrostatic equilibrium configurations of NSs, with particular equations of state, from a non-minimal GMC model which shall be presented in Sect. 2. In Sect. 3 we will derive the hydrostatic equilibrium equations for the concerned theory. In Sect. 4 we will present the EoS that we shall consider for numerically solving the hydrostatic equilibrium equations in Sect. 5. We will discuss our results in Sect. 6.

2 A non-minimal geometry-matter coupling

Here we will work with a GMC theory named f(RL) gravity [5], whose action reads

$$\begin{aligned} S=\int d^{4}x\sqrt{-g}f(R,L), \end{aligned}$$
(2)

with f(RL) being a function of R and L and \(8\pi G\) and c are taken as 1 from now on. One can note from (2) that when \(f(R,L)=R/2+L\), the usual Einstein–Hilbert action is retrieved, such that the variational principle application implies the usual Einstein’s field equations \(G_{\mu \nu }=T_{\mu \nu }\), with \(G_{\mu \nu }\) being the Einstein tensor and \(T_{\mu \nu }\) the energy-momentum tensor.

By following Refs. [13, 14], we will consider the case \(f(R,L)=f_1(R)/2+[1+\sigma f_2(R)]L\), with \(f_1(R)\) and \(f_2(R)\) being functions of R only and the parameter \(\sigma \) can be said to characterize the strength of the coupling. For the sake of simplicity we shall take \(f_1(R)=f_2(R)=R\). Also, we assume \(L=-p\) [35], with p being the pressure of the fluid considered. Taking all these considerations into account, the variational principle applied to (2) yields the following field equations [5]

$$\begin{aligned}&(1-2\sigma p)G_{\mu \nu }+\frac{1}{3}Rg_{\mu \nu }-\frac{\sigma p}{3}Rg_{\mu \nu }\nonumber \\&\quad =(1+\sigma R)\left( T_{\mu \nu }-\frac{1}{3}Tg_{\mu \nu }\right) -2\sigma \nabla _\mu \nabla _\nu p. \end{aligned}$$
(3)

Moreover, the covariant derivative of the energy-momentum tensor reads [5]

$$\begin{aligned} \nabla ^{\mu }T_{\mu \nu }=(-pg_{\mu \nu }-T_{\mu \nu })\nabla ^{\mu }\ln (\sigma R). \end{aligned}$$
(4)

One may wonder, in another perspective, about the cosmological consequences of Eq. (3). In fact, a cosmological model derived from the substitution of the Friedmann–Lemâitre–Robertson–Walker metric as well as the energy-momentum tensor of a perfect fluid in Eq. (3) has not been reported in the literature so far, but shall be soon. For now, we can mention that other GMC models have already been applied to cosmology, yielding well behaved scenarios [37, 49,50,51].

3 The hydrostatic equilibrium equations in a non-minimal geometry-matter coupling model

The hydrostatic equilibrium equations in the concerned GMC theory will be obtained from the substitution of the static spherically symmetric metric

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

in (3)–(4), with \(\alpha (r)\) and \(\beta (r)\) being metric potentials depending on r only.

From the substitution of (5) into (4) and assuming the energy-momentum tensor of a perfect fluid, such that \(T_{\mu \nu }=\texttt {diag}(e^\alpha \rho ,e^\beta p,r^2p,r^2\sin ^2\theta p)\), with \(\rho \) being the matter-energy density, we have for \(\nu =0\)

$$\begin{aligned} \nabla ^{\mu }T_{\mu 0}= & {} (-pg_{\mu 0}- T_{\mu 0}) \nabla ^\mu \mathrm{ln} (\sigma R),\nonumber \\ \nabla ^{\mu }T_{\mu 0}= & {} (-pg_{0 0}- T_{0 0}) \nabla ^0 \mathrm{ln} (\sigma R), \end{aligned}$$
(6)

and considering that we are concerned to a static case, we have \(\nabla ^0 \mathrm{ln} (\sigma R)=0\). Now, for \(\nu =i\), where \(i=1\), 2 and 3, we have

$$\begin{aligned} \nabla ^{\mu }T_{\mu i}= (-pg_{\mu i}- T_{\mu i}) \nabla ^\mu \mathrm{ln} (\sigma R), \end{aligned}$$
(7)

and the right side of (7) provides the term \((-pg_{i i}- T_{i i})\), which is identically null for all values of i. This implies that the energy-momentum tensor is covariantly conserved, i.e.,

$$\begin{aligned} \nabla ^{\mu }T_{\mu \nu }=0, \end{aligned}$$
(8)

which justifies our choice of matter Lagrangian as \(L=-p\), since a different choice of Lagrangian would give a non-conserved energy-momentum tensor.

The 00 and 11 components of the field equations (3) for metric (5) read, respectively,

$$\begin{aligned}&\frac{(1-2\sigma p)}{r^2}\left( r-re^{-\beta }\right) ' +(1-\sigma p)\frac{R}{3}\nonumber \\&\quad =(1+\sigma R)\left( \frac{2}{3}\rho +p\right) +\sigma e^{-\beta }\alpha ' p',\end{aligned}$$
(9)
$$\begin{aligned}&\frac{(1-2\sigma p)}{r^2}\left( e^{-\beta }-1+e^{-\beta }\alpha ' r\right) +(\sigma p-1)\frac{R}{3}\nonumber \\&\quad =(1+\sigma R)\frac{\rho }{3}-2\sigma e^{-\beta }\left( p''-\frac{\beta '}{2}p'\right) , \end{aligned}$$
(10)

with primes denoting derivatives with respect to the radial coordinate r.

The conservation of the energy-momentum tensor (8) yields to

$$\begin{aligned} p'=-(\rho +p)\frac{\alpha '}{2}. \end{aligned}$$
(11)

As the Ricci scalar becomes a degree of freedom, another equation can be derived from the trace of the field equations as

$$\begin{aligned} (1+2\sigma p)R=-(1+\sigma R)T-6\sigma \Box p, \end{aligned}$$
(12)

where the D’Alambertian operator reads

$$\begin{aligned} \Box =-e^{-\beta }\left[ \frac{d^2}{dr^2}-\frac{\beta '}{2}\frac{d}{dr}+\frac{\alpha '}{2}\frac{d}{dr}+\frac{2}{r}\frac{d}{dr}\right] . \end{aligned}$$
(13)

In order to obtain the hydrostatic equilibrium configurations, Eq. (12) needs to be included in the system of differential equations (9), (10) and (11). The unknowns are R, \(\alpha \), \(\beta \), \(\rho \) and p so that we have five variables and four equations. In this way, we need to define a relation between p and \(\rho \), namely an EoS.

4 The equation of state for nuclear matter inside neutron stars

Fig. 1
figure 1

Equation of state properties. Sound velocity, \(v_s\), must respect the causal limit, which states \(v_s<1\). Bottom panel shows that polytropic EoS does not violate causality, while SLy4 does violate for densities slightly above the star central density for the maximum mass star in GR (see Fig. 3b).

Once provided an EoS for the matter inside the star, Eqs. (9)–(12) can be numerically solved. One of the most simple and oftenly used NS EoS in the literature is the polytropic one. Following Refs. [52,53,54,55], the relation between p and \(\rho \) can be regarded as \(p=k\rho ^{\gamma }\), where k is constant and \(\gamma =1+1/n\), where n is the so-called polytropic index. In this work we consider \(k=1.475\times 10^{-3} (\mathrm{fm^3/MeV})^{2/3}\) and \(n=3/2\). Although this kind of equation of state does not represents the state of the art to describe neutron star microphysics, several recent works has used this type of equation of state to study compact stars and, in particular, neutron stars in modified theories of gravity [26, 56,57,58]. The main reason for this is that the particular interest is to investigate the modified theory imprints itself. Polytropic approximations are also often used to describe neutron star EoS by fitting of its parameters. It is worth to cite that, here we are dealing with a polytropic relation between pressure and energy density [59, 60]. This kind of relation is way more efficient to fit the equation of state of neutron stars than a polytropic relation between pressure and rest mass density, and, thus, it is also more efficient to describe the neutron star macroscopic features, such as mass-radius relation. The so called SLy equation of state, described in [67], is a Skyrme type one based on nuclear interactions, and it will also be considered in this work. Both EoS are depicted in Fig. 1, where panel Fig. 1a shows pressure as a function of energy density and panel Fig. 1b presents sound velocity \( v_s = (v_{sound}/c)\) also as function of energy density. From panel Fig. 1a we can see that the polytropic EoS is softer than SLy one, what is the main reason behind the difference on maximum masses in GR, where the polytropic EoS has a maximum mass of \(\sim 1.4M_\odot \) and the SLy \(\sim 2.0M_\odot \). However, both maximum masses in GR are not capable to describe the PSR J2215+5135 mass even if one consider the lower bound. In panel Fig. 1b, the constrain on sound velocity, \(v_s<1\), is not respected for the SLy EoS for densities larger than \(\sim 1800\) MeV/\(\hbox {fm}^3\) (or \(\approx 12.2 \rho _0\), where \(\rho _0\) represents the saturation nuclear energy density), whereas polytropic equation of state violates causality only for densities larger than \(8200~ {\hbox {MeV}}/{\hbox {fm}}^3\, (55.8\rho _0)\).

5 Numerical procedure, mass-radius and mass-central density relations

The system of differential equations previously derived are solved here considering the usual initial and boundary conditions, i.e., the pressure and density have an initial value at the center of the star, namely, \(p(0)=p_C\) and \(\rho (0)=\rho _C\), and the equations are integrated until the pressure vanishes, that is, \(p(r=R_{\star })=0\), where \(R_{\star }\) is the radius of the star [61]. It is worth to quote that R is also zero at the surface according to (12) and, hence, the junction with the exterior Schwarzschild solution is satisfied. The initial value of \(\beta \) is considered to be null, in analogy with the interior Schwarzschild solution. The initial value of \(\alpha \) is arbitrary since the differential equations depend only on its derivative \(\alpha '\).

Fig. 2
figure 2

Mass-radius (a) and mass-central density (b) relations for several values of \(\sigma \) for the polytropic EoS. The magenta circles mark the maximum mass on each curve. In panel (a) the Buchdahl and Schwarzschild limits are indicated by an orange-dotted and a gray-dotted line, respectively.

It should also be cited that since the energy-momentum tensor is covariantly conserved (see Eq. (8)), the total stellar mass can be calculated in its standard form as

$$\begin{aligned} M=\int _0^{R_\star } 4\pi r^2 \rho (r)dr. \end{aligned}$$
(14)

In Fig. 2a we present the mass-radius relation for the polytropic neutron stars in the \(f(R,L)=R/2+L+\sigma R L\) theory of gravitation for different values of \(\sigma \), the coupling parameter. Magenta circles mark the maximum masses for each value of \(\sigma \) and we also present the lines of the Buchdahl and Schwarzschild radius limits, and show that these limits can be surpassed in the GMC theory. From this figure we can see that the radius of the neutron star ranges from 7 to 18 km, which is within the expected values of neutron star radii from observational constraints [62, 63]. On the other hand, the maximum mass of the polytropic neutron star is largely affected, being able to reach values up to \(\sim 2.69M_\odot \).

Fig. 3
figure 3

Mass-radius (a) and mass-central density (b) relations for several values of \(\sigma \) and for the SLy4 EoS. The magenta circles mark the maximum mass on each curve. In both panels a horizontal dashed-dotted line indicates the PSR J2215 + 5135 mass. In panel b vertical dotted line represents the causality threshold.

Furthermore, in Fig. 3 we present the mass-radius and mass-central density relation for the Sly EoS [67]. From panel Fig. 3a, we observe that stars become more massive according \(\sigma \) is increased. In particular, maximum mass increases from \(~2M_\odot \) (\(\sigma =0\)) to \(2.6M_\odot \) (\(\sigma =30\)), which represents a \(\sim 30\%\) increasing. In addition, for the SLy equation of state the maximum mass points are attained for larger values of radius. This is due the maximum mass points are being obtained for smaller values of central energy density, see Fig. 3b. This means that the GMC model yields more massive and larger stars for the SLy EoS, and thus, in this case, the Buchdahl and Schwarzschild limits are not surpassed. This opposite behavior from one EoS to another could help to differ results from different EoS in GMC gravity once given observational values of gravitational redshift. The main reason for this opposite scenario is that used equations of state are quite non-similar. For example, as stated before, the Sly EoS is almost non-causal, while polytropic one is far from being causal for the considered density interval. This means that the larger the softness possibly the larger the sensibility of results concerning to changes in parameter \(\sigma \). It is worth to cite that further investigations analyzing different EoSs could help to better understand this phenomena.

It is worth noticing that the behavior of the mass-radius relation for polytropic NS changes according to the value of \(\sigma \) in a way that a change of concavity is present for \(\sigma \ge 15\). Due to this effect, in the curve for \(\sigma =20\) the mass does not decrease with the reduction of the stellar radius, i.e., dM/dR is always negative.

A similar behavior is presented in the mass-central density relation of Fig. 2b above, where the mass does not decrease with the increasing on central energy density. It is worth to cite that, a priori, this behavior does not represent an instability on the NSs and hence a maximum stable mass cannot be set only by means of equilibrium configurations. In this sense, the GMC model can predict very high mass NSs, such as the pulsar PSR J2215 + 5135, for instance, with a very simple and soft EoS. Those stars also present small radius, such as 8 km. Such a combination of high mass and small radius can overcome the Buchdahl limit (\(GM/c^2R<4/9\)) and also the Schwarzschild radius limit for the star to become a black hole, indicating that the Schwarzschild radius and Buchdahl limit are changed in the GMC model (see also Fig. 2a). We stress that in order to establish a maximum mass limit within the GMC model one needs to perform a study about radial oscillations [64,65,66]. We plan to address all those stability issues in a forthcoming work. However, most massive polytropic neutron stars obtained here have a huge central energy density, as we already pointed out before, surpassing the saturation nuclear density in a few dozens, which may be unrealistic as we are working with a simple EoS treatment. This indicates that, even polytropic stars respecting causality and trace positiveness, the huge maximum mass may not be physically feasible since they are obtained for very high values of central energy density. In addition, when we took the SLy equation of state a remarkable increasing on the maximum mass is presented. For the SLy EoS, we also noticed that maximum central energy densities are smaller than that obtained in GR, which means causality is always respected and the mass of the PSR J2215 + 5135 is attainable for \(\sigma \sim 20\) with maximum central energy density of \(\approx 1000~{\hbox {MeV}}/{\hbox {fm}}^3\, (\approx 6.8\rho _0)\), which is a more feasible scenario.

Anyhow, a maximum limit for the GMC parameter \(\sigma \) can be established from observations concerning the pulsar PSR J2215+5135. For the considered range of central pressure, the value that describes the mass of such a pulsar is about \(\sigma =18\) for the polytropic EoS and about \(\sigma \approx 20\) for the SLy one. The maximum mass for \(\sigma =20\) is \(M_{\mathrm{max}}=2.69M_\odot \), which represents a remarkable increasing of \(\sim 90\%\) with respect to maximum mass point of General Relativity (\(\sigma =0\)) concerning the polytropic EoS.

6 Discussion

It has been shown that it is possible to avoid the Big-Bang singularity through a GMC model [68]. In such an approach, the authors used the so-called Eddington-inspired Born-Infeld gravity. The achievement of evading the Big-Bang singularity was expected to be attained only through quantum gravity and, in this sense, GMC gravity models can figure as a great alternative until we derive the ultimate theory of quantum gravity. For more insights on GMC gravity, one can also check [69].

Let us briefly visit some other recent GMC proposals and their applications. The first GMC model was proposed in [38]. In [70], the authors presented an extension of teleparallel gravity [71, 72], in which the non-metricity is non-minimally coupled to the matter lagrangian. Some cosmological models were derived from such a formalism and have featured an accelerated expansion for late times. In [36], it was shown that the effective energy-momentum tensor of a GMC theory is, indeed, more general than the usual perfect fluid energy-momentum tensor of General Relativity, and that the referred extra terms could be related to elastic stresses in the body, or to other forms of internal energy. In fact, they could also be related to fluid imperfections, such as viscosity and anisotropy. The field equations for a GMC model by using the Palatini formalism were derived in [73]. Furthermore, in [74] it was shown that the dark matter effects can simply be a consequence of GMC and in [75] it was shown that a generalized GMC is compatible with Starobinsky inflation [76].

In the present work, we have obtained hydrostatic equilibrium configurations of NSs in a non-minimal GMC gravity model. The underlying gravitational theory was chosen to be the f(RL) theory [5]. It is the first time in the literature that the hydrostatic equilibrium equations are solved for neutron stars in a fully relativistic treatment in such a theory. In [39, 40] a detailed study was performed for solar-like stars using a similar equation of state as the one polytropic one employed here. However, their equation of state is written in terms of the baryonic mass density instead of energy density, which yields a very different relation between pressure and energy density as showed by [59, 60]. We also consider in the present work the case of SLy EoS , a more realistic and stiff equation of state (EoS) derived from a Skyrme type force based on effective nuclear interaction.

The main motivation for doing so is related to some recent observations of massive pulsars [45, 46]. In fact, even the detection of some super-Chandrasekhar white dwarfs [77, 78] provides a field of investigation for alternative gravity [79, 80]. More remarkably, in [48], a pulsar with \(\sim 2.27M_\odot \), named PSR J2215+5135, was reported as the most massive NS already observed. The achievement of such a mass scale with a very simple EoS was the main motivation of our work, that is, can such a mass scale be predicted through a GMC gravity model with a very simple, and often disregarded in GR, EoS? The answer for the this question is yes, as we explain in the following.

Figure 2a shows the mass-radius relation of polytropic NSs in the GMC model for different \(\sigma \), which features the strength of the GMC. It is clear that for stronger couplings (\(\sigma \ge 18\)), the PSR J2215 + 5135 mass scale is attained. The mass scale of other massive pulsars, such as J1614 – 2230 [45] and J0348+0432 [46], which is \(\sim 2M_\odot \), is attained for slightly weaker couplings (\(\sigma \sim 15\)). In addition, the Sly EoS can also achieve the mass of PSR J2215+5135 for \(\sigma \sim 20\), in spite of its maximum in GR (\(M_{\mathrm{max}}=2.05M_\odot \)).

Moreover, the equilibrium configurations of neutron stars in the GMC theory was showed to largely depend on the EoS stiffness. This aspect of GMC may be clarified if further investigations on the GMC theory with other EoS are performed. For now, we can conclude only that for softer EOS the effect of the coupling between background gravity and matter in GMC generates an effective geometry-matter pressure that is smaller when this coupling is absent, and as a consequence a higher central energy density in the matter pressure is needed to compensate gravity and stabilize the star turning it more compact. For stiff EOS, we have the usual feature of f(RT) or f(R) theories, where the geometry-matter coupling produces an effective pressure that is stronger and a smaller central energy density in the matter pressure is enough to stabilize the star making it bigger and less dense.