Stellar structure of quark stars in a modified Starobinsky gravity

We propose a form of gravity–matter interaction given by ωRT\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega RT$$\end{document} in the framework of f(R, T) gravity and examine the effect of such interaction in spherically symmetric compact stars. Treating the gravity–matter coupling as a perturbative term on the background of Starobinsky gravity, we develop a perturbation theory for equilibrium configurations. For illustration, we take the case of quark stars and explore their various stellar properties. We find that the gravity–matter coupling causes an increase in the stable maximal mass which is relevant for recent observations on binary pulsars.


Introduction
Modern day scenarios such as inflation [1,2], late-time cosmic acceleration [3][4][5], flat rotation curves [6][7][8][9] etc. are incompatible with the standard prescription of general relativity (GR). Although the predictions of GR in the weakfield regime are precise, it falls short in the higher curvature regime in the sense that it predicts singularities such as the big bang and the black hole singularities. It has been shown that quantum corrections generate higher order self-coupling curvature in addition to the original scalar curvature [10,11]. This motivates one to consider non-linear curvature theories to see if they provide a better descriptions of gravitation phenomena.
A nonlinear curvature theory of gravity was proposed by Starobinsky [12] in order to address the issue of the bigbang singularity. He considered the Einstein field equations G μν = κ T μν where the right hand side gives quantum mechanical contributions due to coupling between quantum matter fields (having different spins) with classical gravitational field, with the assumption of isotropy and homogeneity a e-mail: a.mathew@iitg.ac.in b e-mail: m.shafeeque@iitg.ac.in c e-mail: mknandy@iitg.ac.in (corresponding author) and absence of radiation field. In one-loop approximation, and upon regularization, T μν was found to be a function of the Riemann geometric quantities. Based on these findings, Starobinsky exhibited the existence of a one-parameter family of non-singular solutions of the de-Sitter type which could be analytically continued into the region t < 0. The de-Sitter phase naturally explains the inflation scenario without having to include any inflaton field.
However, another approach involves a generalization of the Einstein-Hilbert action where an arbitrary function f (R) represents the Lagrangian density [13]. In the Starobinsky model, namely, f (R) = R + α R 2 , and its other generalisations, inflation has been explained to obtain increasingly better fits the observational data [14,15]. Moreover, various forms of f (R) gravity have been able to explain the late-time cosmic acceleration [15,16]. In addition, a simple powerlaw form of f (R) gravity is able to explain [17,18] rotation curves in the spiral galaxies. The power-law form has been explored [19][20][21] to find a basis for the modified Newtonian dynamics (MOND) which is the most successful scenario in explaining rotation curves in many different types of galaxies [22][23][24][25].
Inclusion of the effect of classical matter with f (R) gravity came in two different forms, namely, f (R, L m ) and f (R, T ), where L m is the matter Lagrangian and T is the trace of energy-momentum tensor. While f (R, L m ) gravity has been studied extensively in various contexts [26,27], f (R, T ) gravity entered the literature somewhat recently [28]. It was noted that the T dependence may arise due to exotic imperfect fluid or quantum effects. Thus it is natural to expect that f (R, T ) gravity may be a suitable candidate for compact objects such as neutron stars and quark stars where quantum effects are expected to play a significant role.
Models of extended gravity have been employed to study the stellar structure of different compact objects. Starobinsky gravity with f (R) = R + α R 2 has been applied to neutron stars treating α R 2 as perturbation with well-known models for the equation of state [29]. It was found for some cases that the maximal mass could approach ∼ 2 M only for negative values of α. Moreover the logarithmic model, f (R) = R + α R 2 + αγ R 2 ln(R/μ 2 ) [30], was also studied perturbatively for neutron and quark stars that exhibited similar trends for different γ values. The Starobinsky model was further explored non-perturbatively for neutron stars [31]. They observed that, for positive values of α, GR yeilded higher maximum mass values than the Starobinsky case. They also studied the model f (R) = R+α R 2 (1+γ R) which exhibited low sensitivity on the γ value. On the other hand, for their model R 1− , GR gave the lowest maximum mass and the mass value increased to very high values approaching 2.5 to 3 M .
Yazadjiev et al. [32] solved for the stable configurations of neutron stars in the Starobinsky model f (R) = R + α R 2 for increasing values of the parameter α. By constructing an equivalent scalar-tensor theory, they obtained the stellar structure non-perturbatively and compared their results with perturbative estimates. While the perturbative result was unphysical because it gave a decreasing mass with respect to the radial distance in a region interior to the star [33], no such unphysical behaviour was observed in the non-perturbative framework. Staykov et al. [34] included a slow rotation in neutron and strange stars in a non-perturbative framework of Starobinky gravity. While the slow rotation does not affect the mass and radius with respect to the static Starobinky case, they found a measurable increase in the moment of inertia with respect to GR.
The Starobinsky model R + α R 2 was further studied for neutron and quark stars non-perturbatively [35]. For positive and non-zero values of α, they observed that the scalar curvature does not decrease to zero at the surface (unlike the perturbative results) and it exponentially falls off outside the star. The stellar mass contribution until the surface plus the gravitational mass contribution outside the star constitute the total mass which is actually observed by a distant observer. The gravitational redshift for the distant observer will be determined by the total (stellar + gravitational) mass. The gravitational mass contribution from the outside of the star is remarkably in distinction with the perturbative approaches where the exterior solution is assumed to be Schwarzschild. For negative values of α, they [35] found that the Ricci scalar executes a damped oscillation beyond the surface of the star and the gravitational mass contribution increases indefinitely. In an earlier paper, the same authors [36] compared the prediction of the Starobinsky model and the corresponding scalar-tensor theory. Their non-perturbative analysis indicated that the star is surrounded by a dilaton sphere whose contribution to the mass is negligible.
Models of f (R, T ) gravity and its generalisations were studied for equilibrium configurations of compact stars. Carvalho et al. [37] considered the model f (R, T ) = R − 2λT to find the equilibrium configurations for white dwarfs. The maximum mass limit obtained was slightly above the Chandrasekhar limit. In comparison to GR and f (R) predictions, the white dwarfs were found to have larger radii as the λ value was increased from zero. Deb et al. [38] considered the same model to obtain the equilibrium stellar structure of quark stars. They demonstrated that the M-R curves are different for positive, negative and zero λ values.
It is important to note that, in the Starobinsky model R + α R 2 , a maximum value of 2 M or beyond is reached only when the α value is chosen to be negative [35]. However, this leads to an issue, namely, the Ricci scalar executes a damped oscillation and the gravitational mass contribution increases indefinitely in the exterior region. On the other hand, the Ricci scalar smoothly decreases to zero at infinity for positive α values, for which the star can support a maximum mass lower than 2 M . Thus a physical theory based on Starobinsky model requires a positive α value whence the Ricci scalar behaves properly everywhere. However, in order to reach 2 M or beyond, the Starobinsky model requires modification. We therefore consider the model f (R, T ) = R + α R 2 + ω RT with α > 0. This modification implies inclusion of gravity-matter interaction in the description via the term ω RT . It would be sufficient to show that the maximum mass attainable is greater than the Starobynsky prediction even if we take a simple form R(1 + α R + ωT ), and treat ωT perturbatively on the background of non-perturbative Starobinsky solution.
In this paper, we obtain the field equations for spherically symmetric distribution of matter for f (R, T ) = R + α R 2 + ω RT to O(ω). With this, we solve for equilibrium configurations of quark stars with the equation of state given by the bag model, namely, p = k(ε − 4B), where B = 60 MeV/fm 3 is the bag constant, and we take the physical value k = 0.28 which is valid for strange quark mass m s = 250 MeV/c 2 . For the pure Starobinsky case, a maximum mass of 1.832 M is obtained for α = 10r 2 g , whereas GR gives a maximum mass of 1.764 M [36]. On the other hand, the present model yields a maximum mass ∼ 2 M , which is consistent with different observations of binary millisecond pulsars, namely, J0348+0432, J1614-2230, and J0740+6620, with pulsar masses 2.01 ± 0.04 M [39], 1.93 ± 0.017 M [40,41], and 2.14 +0. 20 −0.18 M [42], respectively. The paper is organised as follows. In Sect. 2, we lay out the preliminary details for the field equations and energy conservation in f (R, T ) gravity. In Sect. 3, we present the details of calculation for the proposed model with gravity-matter interaction. There, we also develop a perturbative treatment as the gravity-matter interaction is expected to be small. We thus obtain the stellar equation for equilibrium configurations in spherically symmetric stars. We apply these equations to quark stars in Sect. 4 and obtain the stellar properties. Sect. 5 contains a discussion on the obtained results and the main conclusions are given in Sect. 6.

Preliminary details
In this section we briefly present the preliminary details of f (R, T ) gravity needed for our later developments. The action of the most general f (R, T ) gravity is given by [28] where L m is the Lagrangian density of matter. The stressenergy tensor T μν is obtained from the matter Lagrangian L m as Field equations following from Eq. (1) are Assuming the matter to be a perfect fluid, the stress-energy tensor T μν can be obtained from going over to the proper frame and then switching back to the gravitational frame [43], yielding where the energy density ε and pressure p are the proper values and u μ is the macroscopic four-velocity. It can be shown from Eq. (4) that the above form of stressenergy tensor can be obtained from the choice L m = p [44]. Consequently one obtains Thus the field equation (3) become which can be re-written as where G μν = R μν − 1 2 g μν R is the Einstein tensor. It is shown that the covariant divergence of the field equations give the identity [45] which in turn gives Substituting for Θ μν from Eq. (5), we obtain for a perfect fluid.

Present model
In this section we define the present model of gravity-matter interaction in f (R, T ) gravity. We derive the corresponding field equations, the modified TOV equations and also discuss the far-field solution.

Gravity-matter interaction
We consider gravity-matter interaction in a modified gravity represented by where the last term represents the gravity-matter interaction. This form reduces the field equation (7) to The corresponding trace equation is Since the gravity-matter interaction is expected to be small, we shall take a perturbative approach about the exact solutions of R + α R 2 by assuming |ωT | 1. To the first order in ω, we get where the subscript "0" indicates unperturbed quantities when ω = 0, so that φ 0 = 1 + 2α R 0 . The corresponding trace equation (13) is obtained as Since we are interested in the spherically symmetric and static case, we assume the metric ds 2 = −e ν(r ) c 2 dt 2 + e λ(r ) dr 2 + r 2 dΩ 2 (16) along with φ = φ(r ) and T = T (r ), and dΩ 2 = dθ 2 + sin 2 θ dϕ 2 . Thus the above trace equation reduces to The tt-component of the field equation (14) yields The rr-component yields the following equation

Extended TOV equation
Covariant divergence of the field equation (7) yields Substituting f T = ω R, we obtain For the spherically symmetric static metric (16), we obtain from Eq. (4) Since Ricci scalar R and pressure p are functions of r alone, the conservation equation (21) becomes To the first order in ω, we obtain For the case of vanishing ω, we recover the original TOV equation. For ω = 0, we designate Eq. (24) as the extended TOV (ETOV) equation, where the pressure gradient depends on the value of ω as well as the Ricci scalar R. It is thus evident from Eq. (24) that the pressure gradient inside a spherically symmetric star will change as compared to the standard GR case. However, similarly to GR, the cumulative mass m(r ) is related to the metric potential λ(r ) as

Far field solution
In the region exterior to the star, the trace equation (15) takes the form which is identical to that obtained for f (R) = R + α R 2 gravity in vacuum. This suggests that the exterior solution has an identical form in both Starobinsky gravity and for the given particular form of f (R, T ). For spherically symmetric static metric (16), this equation takes the form We note that the Starobinsky correction α R 2 is a very weak contribution as one approaches infinity. This is also immediately obvious from the above equation because the last term dominates in the limit α → 0 giving us back R = 0. Thus the choice of the Starobinsky form R + α R 2 has to coincide with the solution of Einstein gravity at infinity. In order to see how the Einstein limit is approached at infinity, we must do an approximate analysis of Eq. (27). Since ν and λ and their first derivatives are expected to approach zero on approaching infinity (also confirmed by exact numerical calculations), we can approximate Eq. (27) to the form Solution of this equation is given by Since R → 0 as r → ∞, we have to set the integration constant c 2 = 0, giving which approaches zero faster than r −1 as r → ∞ for positive value of α. However, for negative values of α, the far field solution given by (29) is oscillatory in nature implying that negative α values are unphysical.

Quark stars with gravity-matter coupling
In this section, we examine in detail the stellar structure of quark stars in the modified gravity model f (R, T ) = R+α R 2 +ω RT that incorporates gravity-matter interaction.
In massive compact stars (such as quark stars and neutron stars), we expect the gravitational field to be strong enough so that the gravity-matter coupling has a appreciable contribution. With the above choice of f (R, T ) gravity, Sects. 3.1 and 3.2 give a perturbative solution on the background of unperturbed Starobynsky gravity given by The field equations (18), (19), (24), together with the trace equation (17), are reduced to a set of five first order differential equations, given by where we have defined a new field variable Ψ = R .
To complete the solution of the above equations, we take the equation of state of the quark star as that of quark-gluon plasma given by the bag model [46,47] Figure 1 in Ref. [48].
Consistency of the perturbation theory requires |ωT | 1 at all densities. This condition is satisfied throughout the star if one require that |ωT | 1 is true at the center. By defining this condition becomes (1 − 3k)β ε c 4B + 3βk 1. For the choice of β ∼ 10 −2 and ε c 4B ∼ 10, |ωT c | ≈ 4.48 × 10 −2 , thus ensuring the validity of the perturbative approach.
We solve field equations (31)-(35) numerically upon making them dimensionless by defining η = r/r g , χ = Rr 2 g , 5 cm, is taken as the the scaling parameter.
The numerical integrations of the above set of differential equations are carried out by requiring that the metric is asymptotically flat at infinity for the initial conditions λ(0) = 0 and ν(0) = ν c . Since the metric potential ν(r ) enters the field equations only through its derivatives, the central value ν c remains arbitrary and fixed by specifying a value that satisfies ν → 0 for r → ∞. The same initial conditions are imposed on the unperturbed metric, that is, λ 0 (0) = 0 and ν 0 (0) = ν 0c . The central value of pressure p(0) = p c is assigned by the equation of state for a given central density ρ c . Correspondingly, the unperturbed pressure takes the value p 0 (0) = p c . The surface of the star is identified at a radial distance r s (stellar radius) for which the pressure p vanishes. Moreover, since the value of the scalar curvature is maximum at the center, and it gradually decreases towards the surface, we have the boundary condition Ψ (0) = 0.
A unique choice for the central value of the scalar curvature in both the unperturbed scenario (Starobinsky gravity) and the actual case requires the knowledge of the exterior solution. This requires one to continue the integration outside the star with boundary conditions λ(r s ) = λ s , ν(r s ) = ν s , R(r s ) = R s and Ψ (r s ) = Ψ s at the surface, obtained from the interior solution for an initial guess for R c . To set the above initial conditions, we first carry out a numerical integration for the unperturbed case, with similar boundary conditions λ 0 (r s ) = λ 0s , ν 0 (r s ) = ν 0s , R 0 (r s ) = R 0s and Ψ 0 (r s ) = Ψ 0s , imposed for an initial guess R 0c .
For convenience, the initial guess in both the cases are taken to be the GR value E R c = κ(ρ c c 2 − 3 p c ) since this value is not too far from the required values. The integration is carried out several times for different initial guesses until the required conditions R → 0 and |Ψ | → 0 as r → ∞ and R 0 → 0 and |Ψ 0 | → 0 as r → ∞ are satisfied. This procedure gives the correct central values R c and R 0c , which are found to be lower than the GR value. For the sake of accuracy, we fine-tune χ c (= R c r 2 g ) and χ 0c (= R 0c r 2 g ) up to twelve decimal figures by taking the exterior solution as far as η max (= r max /r g ), where η max satisfies χ(η max ) ∼ 10 −12 and χ 0 (η max ) ∼ 10 −12 .
In the original framework of general relativity, the Ricci scalar vanishes immediately outside the surface, the trace equation being R = −κ T . In our present case, the vacuum solution is given by Eq. (26) and the Ricci scalar does not vanish but decays exponentially outside the star, as also implied by the far-field solution, Eq. (30). This gives rise to two distinct masses [36], namely, the stellar mass M s = m(r s ), the mass within the stellar radius r s , and the mass M as seen by a sufficiently distant observer, estimated as Since the numerical calculations are sufficiently accurate with a sufficiently large value of r max , this estimate for M is expected to be close to the one for r → ∞.
In the following subsections, we analyze the exact numerical solutions of the field equations given by Eqs. (31)-(35) with the boundary conditions discussed above for quark stars with the equation of state given by the bag model. We also compare these results with the Starobinsky case, ω = 0. We take α = 10r 2 g = 2.1804 × 10 11 cm 2 (that is, √ α = 3.16r g = 4.6694 × 10 5 cm), which is smaller than the estimated upper bound √ α < 7 × 10 7 cm as predicted by binary pulsar data [49].

Interior and exterior solutions
In this sub section we elaborate upon the interior and exterior solutions by looking at the radial profiles for pressure p(r ), mass m(r ) and Ricci scalar R(r ).  Radial profiles of the scalar curvature R(r ) for different values of β with central density ρ c = 2.5414 × 10 15 g /cm 3 are shown in Fig. 2. It may be observed that the choice of the central scalar curvature R c is strongly correlated with β (or equivalently ω) in that the value of R c increases with increasing values of ω. On the other hand, when we fix β = 0 and vary ρ c , the central value R c is found to be higher for a higher value of ρ c as shown in the inset of Fig. 2. Although there is an increase in the value of R c when ρ c is increased, the scalar curvature R falls off rapidly for higher value of ρ c than for the lower one. In the former case (with fixed ρ c and varying β), the Ricci scalar maintains higher values thoughtout the star for higher value of β.This is consistent with the fact that the perturbative ω terms in Eq. (32) add with the term κ 6α T e λ , so that the effective value of T changes which is equivalent to a changed value of matter content with respect to the Starobinsky case. Since these perturbative terms disappear outside the star, they act as if they were an additional matter content in Starobinsky gravity. Figure 3 shows the total mass profile m(r ) with ρ c = 2.5414 × 10 15 g /cm 3 for different values of β. The inset represents the mass profile up to the stellar surface r = r s . In comparison with Starobinsky gravity (β = 0), we observe a slight increase (decrease) in the stellar mass M s and stellar radius r s when β is positive (negative). Further we observe that these changes in stellar mass M s and radius r s are larger for β = ±0.025 than β = ±0.01. On the other hand, we see that the mass M measured by a distant observer increases appreciably with increasing β.
As seen from Fig. 2, the scalar curvature does not decrease to zero at the surface of the star and it falls off outside the star. This fall-off is similar to a Yukawa function (as shown in Sect. 3.3). There is a gravitational mass contribution due to the non-vanishing scalar curvature outside the star. The mass profiles shown in Fig. 3 contains both contributions, stellar plus gravitational. The inset shows only the stellar mass contribution that does not extend beyond the stellar radius r s ∼ 10 km. The main graphs in Fig. 3 show both contributions (stellar plus gravitational) extending beyond r s ∼ 10 km. We see that the combined mass profiles approach asymptotic values for large r (∼ 60 km). A sufficiently distant body experiences the gravitational field of the combined mass.

Mass-radius relations
In this section we study the mass-radius (M − R) relations obtained from the field equations for a continuous range of central density ρ c or equivalently the central Ricci scalar R c . In addition, we verified that all energy conditions are satisfied. Figure 4 represents the relation between stellar mass M s and stellar radius r s for different values of β. We see that for In the same regime, if we fix r s , M s is found to increase with increasing β. This fact signify that the presence of the ω terms strengthens the effect of gravity to balance the increased pressure gradient as inferred from the pressure profiles studied in Sect. 4.1. Figure 5 presents the mass M measured by a distant observer against the stellar radius r s for different values of β. As noted earlier, the mass M consists of contributions from both the stellar mass and non-vanishing scalar curvature extending beyond the stellar radius r s . This mass was calculated up to a sufficiently high radial distance (r max ) until the scalar curvature approached very close to zero with the condition given by Eq. (37). The relationship between the curves in Fig. 5 bear similarity with those in Fig. 4 giving qualitatively similar conclusion. However we note the important fact that maximum observed mass M * are appreciably higher than those in Fig. 4.

Stability and energy conditions
In this section, we study mass versus central density to find the maximal mass from the stability of equilibrium configurations. The stable configuration corresponds to the region in the mass-central density curve where ∂ M ∂ρ c > 0, whereas the unstable region is given by ∂ M ∂ρ c < 0 [50]. The onset of instability is identified as the point where ∂ M/∂ρ c = 0, and the mass corresponding to this point is the maximal.
To study the stability, we first examine stellar mass M s versus central density ρ c in Fig. 6 for different values of β. We see that the maximal stable mass M * s (corresponding to the maximal total mass M * ) increases as β increases from β = −0.025 to β = +0.025.
A similar trend is observed in Fig. 7, where the mass M observed by a distant observer is plotted against the central density ρ c . Here we see that the maximal mass M * (denoted by the open circle in the figure) shift to appreciably higher value with respect to the stellar values M * s . Table 1 displays maximal mass values M * observed by a distant observer for different values of β. The corresponding stellar mass values M * s , stellar radius r * s , central Einstein Ricci scalar E R * c , central Starobynsky Ricci scalar R * 0c , central Ricci scalar R * c , and central density ρ * c are also displayed. We see that the maximal mass value M * increases and approaches ∼ 2 M as β is increased. We note that this increase is appreciable even for very small magnitudes of β, suggesting a measurable effect played by gravity-matter interaction.
We verify the validity of the perturbative results by estimating the maximal value of ωT corresponding to the maximal mass. For β = 0.025 (or ω = 0.025/4B) and central density ρ * c = 2.4173 × 10 15 g/cm 3 , we get |ωT c | = (1 − 3k)βρ c c 2 /4B + 3βk = 4.36 × 10 −2 . Thus the maximum value of ωT is of the order of 10 −2 , giving assurance to the validity of the perturbative results. Figure 8 displays the energy conditions [51], namely, null energy condition (ε ≥ 0), weak energy condition (ε+ p ≥ 0), strong energy condition (ε + 3 p ≥ 0), and dominant energy condition (ε −| p| ≥ 0). We see that all energy conditions are satisfied because they are positive in the entire region of the star. These energy conditions are valid to a good approximation since they are large throughout the star (lying between ∼ 1 and ∼ 10 in the units of 4B) compared to the highest perturbation (∼ 10 −2 at the centre).

Discussion
In the original case of Einstein's gravity, the Ricci scalar is a linear function of the central density, expressed as χ = κr 2 g {(1 − 3k)ρc 2 + 12k B}. One might expect such a linear relationship in the Starobinsky model in the low density regime where α R 1. Besides, the same behaviour is expected in the present model for low central densities where the term R dominates over α R 2 and ω RT . However, the situation is completely different in the high density regime in both Starobinsky and the present model. In the Starobinsky model, we found that at higher central densities, the central Ricci scalar R 0c varies slowly as a function of ρ c , as shown in Fig. 9. On the other hand, in the present model, depending on the sign of β (or ω), we find that the central Ricci scalar R c would either increase or decrease with respect to the Starobinsky case. Figure 9 shows the variation of central Ricci scalar with respect to central density for different values of β. For positive β, we see that the R c value increases compared to the Starobinsky model in the higher density regime. This was expected since the additional terms in the present model contribute at higher densities as previously noted in Sect. 4.1. The opposite is true for the case of negative β val-   ues where we found that the central Ricci scalar R c values lie below the Starobinsky values for higher densities as shown in Fig. 9. This happens because O(ω) term gives a negative contribution in this case. For any positive β, the curve in Fig. 9 lies above the Starobinsky case, so that maximum mass values higher than the Starobinsky case would be obtained. On the other hand, the curve for any negative β lies below the Starobinsky case implying that the maximum mass values lower than the Starobinsky case would be obtained.
It may be recalled from the discussion in Introduction that, in the Starobinsky gravity (ω = 0), the star is surrounded by a gravitational halo since the Ricci scalar is non-vanishing outside the star. In fact we see the same behaviour from Fig. 2 in the present model (ω = 0) as well. The far-field solution (given by Eq. 30) shows that the Ricci scalar (and the gravitational halo) falls off exponentially as r → ∞. Moreover, it is evident from Fig. 2 that this fall-off is faster for higher values of the central density ρ c . Figure 10 Fig. 10 The radius of gravitational halo r H versus the central density ρ c until the onset of gravitational instability at ρ * c . The inset shows that the radial profiles for the scalar curvature χ (corresponding to the points denoted by the open circle in the main graph) falls off faster for a higher density, similar to the Starobinsky case (β = 0), as shown in the inset of Fig. 2 graphs in Fig. 10 that the radius of the halo r H decreases with increasing central density ρ c . At the same time, the value of r Sch = 2G M c 2 increases with increasing central density ρ c , as seen from Fig. 7. These two opposite behaviours (shrinking and expansion) continue as ρ c increases. Thus it is apparent that, when the star collapses to a black hole (with an infinite central density), the gravitational halo would shrink and will be well inside the horizon, leading to a vanishing Ricci scalar outside the horizon. This scenario is consistent with the fact that when the coefficient of R 2 term is positive, the only static spherically symmetric solution of a black hole with a regular horizon is the Schwarzschild solution, as shown in Ref. [52].

Conclusion
In this paper, we considered a form of f (R, T ) gravity that includes a coupling between gravity and matter on the background of the Starobinsky model. While the Straobinsky model takes account of quantum fluctuations [12] and f (R, T ) gravity may arise due to quantum effects [28], a coupling between matter and gravity is expected to bring about the features in quark stars where the gravitational field is extremely strong. In particular, the stellar structure of quark stars, with equation of state coming from the bag model, is expected to undergo an measurable change due to this coupling. Moreover, we speculate that the maximum mass limit would change appreciably so that astrophysical observations on binary pulsars could be given a theoretical basis.
The Starobinsky model has been applied to quark stars by other authors [35,36] to find their stellar structure. This produced a different stellar structure from the pure GR case and it was found that the maximum mass limit increased from the GR case due to an additional contribution from gravitational mass enveloping the stellar mass. In the present case, we find that this mass is further increased due to additional contribution from the coupling between gravity and matter (for positive values of ω).
To assert the above features, it was sufficient to treat the gravity-matter coupling as a perturbation keeping in mind that the coupling constant is sufficiently small for the validity of the perturbation treatment. We adopted this perturbation treatment in the background of unperturbed solutions of the Starobinsky case. Remarkably, such a treatment gives physically acceptable solutions for both signs of the coupling constant ω representing the strength of gravity-matter interaction.
The gravity-matter coupling term increases (decreases) the magnitude of the pressure gradient p (r ) for positive (negative) values of ω pushing (pulling) the stellar boundary outward (inward) as compared to the pure Starobinsky case. Moreover the strength ω of the gravity-matter coupling determines the central value of the scalar curvature R c for a given central density ρ c . The scalar curvature maintains higher (lower) values throughout the star compared to the pure Starobinsky case for positive (negative) values of ω as the effective matter content increases (decreases) within the star. It is interesting to see that, although there is small increase (decease) in the stellar mass M s for positive (negative) values of ω, the gravitational mass contribution enveloping the star increases (deceases) appreciably with respect to the Starobinsky case. This is because of the increased (decreased) scalar curvature exterior to the star contributing a greater (lesser) gravitational mass than the Starobinsky case. Consequently the quark star can support higher values of maximal total mass (M * ) than the Starobinsky case for positive values of ω.
Recent observations of binary millisecond pulsars, have yielded the pulsar masses to be ∼ 2 M [39][40][41][42]. Such a high value of mass cannot be explained by models based on hyperon or boson condensate equations of state for neutron stars, leading to their possibility of being quark stars. We see that our present model with gravity-matter coupling, although treated as a perturbation, is capable of supporting high values of masses of quark stars.
Acknowledgements Arun Mathew is indebted to the Indian Institute of Technology Guwahati for financial assistance through a Research and Development fund. Muhammed Shafeeque is indebted to the Ministry of Human Resource Development, Government of India, for financial assistance through a doctoral fellowship. The authors would like to thank the anonymous reviewer for constructive comments and suggestions that improved the presentation in the paper.

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: This article represents a theoretical work and the ensuing numerical results are included in detail in this article. For the purpose of comparison, observational data from other sources are quoted and cited in this article.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .