Non-trivial class of anisotropic compact stellar model in Rastall gravity

We investigated Rastall gravity, for an anisotropic star with a static spherical symmetry, whereas the matter-geometry coupling as assumed in Rastall Theory (RT) is expected to play a crucial role in differentiating RT from General Relativity (GR). Indeed, all the obtained results confirm that RT is not equivalent to GR, however, it produces the same amount of anisotropy as GR for static spherically symmetric stellar models. We used the observational constraints on the mass and the radius of the pulsar \textit{Her X-1} to determine the model parameters confirming the physical viability of the model. We found that the matter-geometry coupling in RT allows slightly less size than GR for a given mass. We confirmed the model viability via other twenty pulsars' observations. Utilizing the strong energy condition we determined an upper bound on compactness $U_\text{max}\sim 0.603$, in agreement with the Buchdahl limit, whereas Rastall parameter $\epsilon=-0.1$. For a surface density compatible with a neutron core at nuclear saturation density, the mass-radius curve allows masses up to $3.53 M_\odot$. We note that there is no equation of state is assumed, however, the model fits well with a linear behavior. We split the twenty pulsars into four groups according to the boundary densities. Three groups are compatible with neutron cores while one group fits perfectly with higher boundary density $8\times 10^{14}$ g/cm$^3$ which suggests that those pulsars may have quark-gluon cores.

The GR theory has been proven to be a successful theory of gravity on solar system scales (spacetime curvature is weak) by many observational tests, c.f. [38], and also on black hole scales (spacetime curvature is extremely strong) using black hole shadows observations by Event Horizon Telescope [39] not to mention the perfect predictions of the gravitational waves due to binary black holes merging [40][41][42]. On the cosmological scales, however, the GR does not provide answers for explaining the late accelerated expansion [43][44][45][46]. Even in presence of a cosmological constant Λ, the discrepancy of the current Hubble parameter H 0 value, between early universe observations by Planck satellite and late universe measurements by distance ladder or strong lensing, may point out the need to modify the GR theory [47], see also [48][49][50].
Many efforts have been done to generalize GR theory by using general function in Einstein-Hilbert action instead of the Ricci invariant, e.g. f (R), f (G), f (T ) and mimetic gravity [48,49,. In fact, these modified theories kept the fundamental assumption that the covariant divergence of the energy-momentum vanishes, i.e. T α β;α = 0 where the semicolon denotes the Levi-Civita covariant derivative. On the contrary, Rastall attempted to modify GR by dropping this assumption replacing it by setting T α β;α = a β where a β vanishes in flat spacetime (vacuum) and recovers GR, otherwise it does not [72,73]. Rastall showed that a β ∝ ∂ β R is a reasonable choice which reflects the non-minimal coupling between matter and geometry. Interestingly, some cosmological models have been constructed using RT [74][75][76] as well as black hole solutions [77][78][79][80][81]. Moreover, shadows of rotating black holes have been investigated within RT by Kumar et al. [82] and Heydarzade and Darabi [77]. Furthermore, the effects of the non-minimal coupling between matter and geometry on thermodynamics of black holes in RT has been investigated in comparison with the GR theory [79,[83][84][85]. a e-mail: nashed@bue.edu.eg b e-mail: waleed.elhanafy@bue.edu.eg Recently, Vissar claimed that RT is completely equivalent to GR [86]. On the contrary, Darabi et al. investigated Visser's claim but they concluded that Visser misinterpreted the matter-geometry coupling term which led him to wrong conclusion [87]. In addition, they showed that by applying Visser's approach to f (R) theory one may conclude that it is equivalent to GR as well which is not true. Different studies have proven that RT is not equivalent to GR [88][89][90][91]. Visser's conclusion is correct when Ricci scalar vanishes for black holes in general, otherwise the claim is incorrect and both theories are not equivalent. One of the good examples which may reveal the contribution of the matter-geometry coupling in RT in contrast to GR is the stellar models when the presence of matter plays a crucial role. It is the aim of the present study to derive a anisotropic static spherically symmetric interior solution using RT and confront it with pulsars observations. The arrangement of this study is as follows: In Section 2 we briefly review the main assumptions and features of the nonconservative theory of Rastall's gravity. In Section 3, we apply Rastall's field equations to a spherically symmetrical object in presence of anisotropic matter source. The system consists of three non-linear differential equations in five unknown functions (two metric potentials in addition to the energy-density, radial and tangential pressures). Therefore, two constraints have been imposed to close the system without using equation of state (EoS): We assumed a specific form of one of the metric potential, g rr as usually done in interior solutions, additionally we made the metric potential g tt contribution to the anisotropy to vanish. Consequently, we derive the analytic form of energy-density, radial, and tangential pressures which satisfy Rastall field equations. In Section 4, we list the necessary and sufficient physical conditions that any stellar model must fulfill to be compatible with a realistic compact star. In Section 5, we investigated the viability of the obtained solution with the conditions listed in Section 4. In Section 6, we matched our model with the exterior vacuum solution (Schwarzschild solution) fixing the model parameters and constants. The pulsar Her X-1 , with mass M = 0.85 ± 0.15 M ⊙ , and radius R = 8.1 ± 0.41 km [92] is used to determine the constants of the model. In Section 7, we investigated the stability of the model using a modified version of Tolman-Oppenheimer-Volkoff (TOV) equation of hydrostatic equilibrium and the adiabatic indices. In Section 8 we use other observational mass-radius constraints of twenty pulsars' data to test the viability of the model. Also we plot the mass-radius profile of the model for different surface densities compatible with neutron core pulsars to show the maximum masses allowed by the model. Finally, we summarize and conclude in Section 9.

Rastall gravitational theory
In Riemann geometry, by making use of the contracted Bianchi Identity on one hand and the minimal coupling procedure on the other hand, this led Einstein to formulate the consistent field equations of GR where χ = 8πG N /c 4 where G N is the Newtonian gravitational constant and c is the speed of light, G αβ denotes Einstein tensor, R αβ denotes Ricci tensor and R = g αβ R αβ denotes Ricci invariant. Rastall, however, dropped the minimal coupling procedure assuming non-divergence-free energy-momentum in curved spacetime [72,73] T α where the constant of proportionalityǫ measures how much the conservation law is locally violated. According to this assumption Rastall obtained a consistent set of field equations Obviously Rastall gravity reflects non-minimal coupling features of matter and geometry whereas Einstein's field equations are recovered when the spacetime is flat. Alternatively, Eq. (4) can be rewritten as Contracting the above equation gives where T = g αβ T αβ is the trace of the energy-momentum tensor. Thus the field equations of RT read where It proves convenient to use a dimensionless Rastall's parameter ǫ =ǫχ, c.f. [93]. Then, Eq. (5) becomes and the tensor T αβ in Eq. (7) reads For ǫ = 0 case, the conservation law is restored and the GR version of gravity is recovered. In this sense, RT generalizes Einstein's one by assuming a local violation of conservation law in curved spacetime due to non-minimal coupling between matter and geometry. Otherwise, flat spacetime, both theories are equivalent. Therefore, one of the important applications, which differentiate both theories, is stellar structure models when presence of the matter sector plays a crucial role in interior solutions.

spherically symmetric interior solution
Providing that the static spherically symmetrical spacetime is given by the following metric 1 where F(r) and G(r) are unknown functions. For the spacetime metric (11), we obtain the non-vanishing components of Levi-Civita connection where prime (double prime) denotes first (second) derivative with respect to the radial coordinate. We write the Ricci invariant We assume the energy-momentum tensor for a anisotropic fluid with spherical symmetry, i.e.
where ρ = ρ(r) is the fluid energy density, p r = p r (r) its radial pressure (in the direction the time-like four-velocity u α ), p t = p t (r) its tangential pressure (perpendicular to u α ) and ζ α is the unit space-like vector in the radial direction. Then, the energy-momentum tensor takes the diagonal form T α β = diag(−ρ, p r , p t , p t ). Applying Rastall's field equations (7) to the spacetime (11) where the matter sector is as given by (14) we obtain, respectively, the components t t, r r and θ θ (= φ φ) as follows: Additionally, we define the anisotropy of the system (15) using the parameter We note that the Rastall parameter quantitatively contributes in the radial and the tangential pressures by equal amount, i.e. Rastall parameter has no contribution in the anisotropy parameter, and this conclusion is valid for any spherically symmetric interior solution. For the ǫ = 0 case, the differential equations (15) coincide with Einstein field equations of an interior spherically symmetrical spacetime [94,95]. The system (15) consists of three independent non-linear differential equations in five unknowns G, F, ρ, p r and p t . Therefore, two conditions need to be imposed to close the above system. We follow the ansatz given in [96] (also used in [95]). First, we assume the metric potential G to have the form with a 2 is a dimensionless constant to be determined by boundary condition and R is the radius at the star boundary. We note that the above ansatz is regular everywhere inside the star, i.e. 0 ≤ r ≤ R, where |a 2 | < 1. Substituting (17) in the anisotropy parameter (16), we get Now, we impose the second condition by assuming that the component g tt has no contribution on the anisotropy parameter, i.e.
This choice clearly gives no anisotropy at the center, r = 0, which is physically a reasonable feature [96], see the physical conditions in Sec. 4. Using Eqs. (18) and (19) and by solving for the metric potential F, we obtain: where the constants of integration a 0 and a 1 are dimensionless to be fixed by matching conditions. Up to this step the obtained results are the same as given by [96]. However, the matching conditions with the exterior solution give the set of constants {a 0 , a 1 , a 2 } in terms of Rastall parameter ǫ. Therefore, any deviation from the GR in the following results should be related to the assumed matter-geometry coupling in Rastall field equations (15), i.e. Rastall parameter ǫ. Substituting the metric potentials (17) and (20) into the system (15), we get the energy-density, radial and tangential pressures in the form Equations (21) coincide with the GR version when Rastall parameter ǫ vanishes [94]. It is to be mentioned that the anisotropic force, F a = 2∆ r , becomes attractive if p t − p r < 0 and repulsive if p t − p r > 0. The mass contained within a radius r of the sphere is defined as Using the energy-density as defined in Eqs. (21) and the above equation (22), we get where ℵ = (a 0 + 2a 1 a 2 2 )a 1 . It proves convenient to use the compactness parameter [97,98] of a spherically symmetric source with radius r, to study the stability of compact objects. Similarly we use the gravitational red-shift parameter Z which is related to the metric potential as These parameters will be used later in Sec 6. In the next section, however, we are going to discuss the physical requirements to derive viable stellar structures in order to test the viability of the obtained model.

Physical conditions for a stellar model
For a stellar model to be physically well behaved, it needs to satisfy the following conditions: (i) For the geometric sector, the metric potentials F and G should be free from coordinate and physical singularities within the interior region of the star 0 ≤ r ≤ R, where the center (boundary) is at r = 0 (r = R) respectively.
(ii) The metric potentials of the interior solution and the exterior 2 should match smoothly at the boundary.
(iii) For the matter sector, the fluid density, radial and the tangential pressures should be free from coordinate or physical singularities within the interior region of the star. In addition, they should be maximum at the center of the star and monotonically decrease towards the boundary of the star. i.e.
(v) At the boundary of the star (r = R), the radial pressure should vanish, i.e. p r (r = R) = 0. However, the tangential pressure at the boundary should not necessarily vanish.
(vii) The fluid density, radial and tangential pressures should fulfill the following energy conditions: a. Null energy condition (NEC): ρc 2 + p t > 0, ρ > 0, b. Weak energy condition (WEC): ρc 2 + p r > 0, ρ > 0, c. Dominant energy conditions (DEC): ρc 2 ≥ |p r | and ρc 2 ≥ |p t |, d. Strong energy condition (SEC): (viii) The causality condition should be satisfied, that is the speed of sound should be smaller than unity everywhere inside the star and monotonically decrease toward the boundary, i.e. for the radial velocity 0 ≤ v r /c = 1 c dp r dρ ≤ 1 and v ′ r 2 < 0, and for the tangential (ix) The stability condition should be satisfied, i.e. −1 < (v 2 t − v 2 r )/c 2 < 0 within the star [99]. (x) The gravitational red-shift should be finite and positive everywhere inside the star and decreases monotonically toward the boundary, i.e. Z > 0 and Z ′ < 0. (xi) The adiabatic index stability condition for anisotropic star should be fulfilled, i.e. the adiabatic index Γ > γ where γ = 4/3 is the adiabatic index corresponds to the isotropic case.
We note that the stellar model which fulfills the above mentioned conditions is physically viable and well behaved. In the following sections we are going to examine the model at hand with these conditions investigating possible roles of Rastall parameter.

Physical properties of the model
We test the model showing that the metric potentials (17) and (20), and the solution (21) are free from physical or geometric singularities. Then we derive some physical quantities which will be used later in Sec. 6 to examine the viability of the present model. In addition, we use the matching conditions at the boundary surface of the star to set the constraining equations on the model parameters {ǫ, a 0 , a 1 , a 2 }.

Non singular model
From Eqs. (17) and (20) one finds that the metric potentials at the center read This ensures that the gravitational potentials are finite at the center of the star. Moreover, the derivatives of these potentials are finite at the center, i.e. F ′ (r = 0) = G ′ (r = 0) = 0. Equation (26) ensures that the metric is regular at the center. This satisfies condition (i) in Sec. 4. From Eqs. (21) one finds that the density, radial and tangential pressures at the center are , .
These ensure that the anisotropy parameter has a vanishing value at the center as stated in condition (iv) on Sec. 4. The density and the pressures should be positive at the center which sets two constraints on the model parameters/constants. Additionally, the Zeldovich condition [100] states that the radial pressure must be less than or equal to the density at the center, i.e. p r (0) Using Eqs. (21) we give the derivative of energy density, radial and tangential pressures, respectively, as follows We use Eqs. (29)- (31) to show that the gradients of the energy-density, radial and tangential pressures are negative later in Sec. 6.

Matching conditions
We note that the exterior spacetime of a static spherically symmetric star is the same for both GR and RT, since the exterior region is vacuum. Thus no reason to expect any solution rather the exterior Schwarzschild one for Rastall's theory, that is where M is the total mass r > 2M. We are going to match the interior spacetime metrics (17) and (20) and the exterior Schwarzschild spacetime metric (34) at the boundary of the star r = R. Therefore, the continuity of the metric functions, as stated by condition (ii) in Sec. 4, across the boundary gives the conditions In addition, the radial pressure (21) approaches zero at the star boundary, p r|r=R = 0, which reads The above constraint ensures that condition (v) in Sec. 4 is fulfilled. From the above conditions, namely (35) and (36), we get the constraints on the set os constants {a 0 , a 1 , a 2 } in terms of the start mass M, radius R in addition to the Rastall parameter ǫ. Using observational pulsars data, knowing the observed values of M and R, we obtain the corresponding numerical values for a particular choice of ǫ.   (19), anisotropic force F a = 2∆/r. We note that the Rastall parameter has no contribution in the anisotropy, therefore GR and RT predicts same anisotropy in the case of spherical symmetry as discussed after (16). For ǫ = 0 and ǫ 0, the gradients of the density, tangential and radial pressures given by Eqs. (29)-(31) versus the radial coordinate r in km of the pulsar Her X-1.

Astrophysical observational constraints
We convert back to physical standard units in order to correctly determine the numerical values of the model parameters. Using masses and radii of observed pulsars along with the physical conditions mentioned in the previous section, one can determine the constant parameters of the obtained model (21) and test its viability. We use the observational constraints of the particular pulsar Her X-1, whose mass M = 0.85 ± 0.15M ⊙ and radius R = 8.1 ± 0.41 km [101], where M ⊙ (= 1.989 × 10 30 kg) denotes the solar mass. Then, the boundary conditions (35) and (36)  Noting that we select a 2 < 1 which is required by the regularity condition of ansatz (17). Substituting the above expressions into Zeldovich condition (28), keeping in mind that the RT predictions are not expected to be far from GR ones (i.e. ǫ should be small), we obtain the following constraints on Rastall parameter −1.880 ǫ 0.259.
Next, we use the above expressions to plot some physical quantities for ǫ = 0 (GR case) and ǫ 0 (RT case) which enable us to test the model viability and investigate the role of Rastall parameter. We note that one should be careful for the choice of Rastall parameter in order to satisfy the physical conditions in Sec. 4 with a particular compact object. In agreement with Oliveira et al. [93] we confirm that the present model is physically stable whereas the Rastall parameter ǫ is negative 3 . For the particular pulsar Her  In Figs. 1(a)-(c), recalling condition (iii) in Sec. 4, we represent the behavior of energy-density, radial and tangential pressures for the pulsar Her X-1 from its center to the boundary surface. This shows that the density, radial and tangential pressures are positive as required for realistic stellar configuration with maximal finite values at the center decrease monotonically towards the surface. Also, as Fig. 1 (a) shows, the surface density is above 7 × 10 14 g/cm 3 while the central density is comparable to the nuclear density ∼ 10 15 g/cm 3 which indicate that the pulsar Her X-1 may has no neutron core but quark-gluon one. This will be revisited in Sec. 8. In addition, Fig. 1(b) shows that the radial pressure at the surface is null in agreement with condition (v) in Sec. 4.
In Fig. 2(a), recalling condition (iv) in Sec. 4, we show that the anisotropy parameter (19) is null at the center and increases monotonically toward the surface of the star. Also, it shows that the anisotropic force F a = 2∆ r is positive reflecting the repulsive behaviour of the force as p t − p r > 0. Recalling Eq. (16) we verify that the Rastall parameter has no contribution to the anisotropy parameter so both GR (ǫ = 0) and RT (ǫ 0) cases coincide. However, this is not the case with other ingredients. For the GR (ǫ = 0) and RT (ǫ 0) cases, Figs. 2 (b) and (c) respectively confirm that the gradients of energy-density, radial and tangential pressures, given by Eqs. In Figs. 5(a) and (c), for ǫ = 0 and ǫ 0 cases, we plot the evolution of the radial and the tangential EoS parameters, i.e. w r (r) = p r /ρ and w t (r) = p t /ρ, respectively, verses the radial distance from the center. The plots show that the EoS parameters slightly vary from maximum at the center with a monotonically decreasing behavior toward the surface whereas w r drops to zero at the surface as expected. We note that the variation of the EoS parameters is more severe in RT than GR, which should be understood as an effect of the matter-geometry coupling. On the other hand, we remind that there are no EoS are preassumed in the model setup, Figs. 5(b) and (d), however, prove that the EoS p(ρ) in both radial and tangential directions is in well fit with linear patterns, whereas the best fit are given by p r = 0.414 ρ − 27.6 and p t = 0.223 ρ − 11.8 for the pulsar Her X-1. Notably the slopes of the fitted lines dp r /dρ = 0.414 = v 2 r and dp t /dρ = 0.223 = v 2 t are consistent with the corresponding values of Table 2 for the pulsar under investigation. Fig. 6 (a) which shows that it is a monotonically increasing function of the radial coordinate whereas the mass at the center M(r = 0) = 0. The plot shows that theoretical M(r) curve of the pulsar Her X-1 according to RT is compatible with the observed values of the mass and the radius. In Fig. 6 (b) we show the behavior of the compactness parameter of the pulsar. Interestingly the RT curve allows less size than the GR for same mass which may evidently reflects the role of the matter-geometry coupling to obtain slightly higher compactness. Also Fig. 6 (c) shows the redshift of the pulsar is maximum at the center and decreases monotonically toward the surface. Notably the redshift at the surface is found to be Z R ∼ 0.204 which is in an agreement with the upper bound constraints Z R ≤ 2 as given by Buchdahl [102], for anisotropic spheres see [103,104] and for the upper bound on the surface redshift in presence of a cosmological constant see [105].

Stability of the model
In addition to the stability condition (ix) in Sec. 4 which has been shown to be verified we are going to discuss the stability of the obtained stellar model using two different techniques in the present section; that are the modified TOV equations and the adiabatic index.  The mass function plot confirm that the model can predict the a of the pulsar Her X-1 in agreement with observational data. The plot shows that RT predicts compactness values higher than GR which reflects the role of the matter-geometry coupling to allow more compactness values. The redshift is finite everywhere within the pulsar and decreases toward the surface as stated by condition (x) in Sec. 4 and also predict a surface redshift consistent with the upper limit constraints as given by [105].

Equilibrium analysis via Tolman-Oppenheimer-Volkoff equation
We assume hydrostatic equilibrium to be everywhere within the stable compact star. This configuration, then, can be described by the GR based TOV equation [106][107][108] which gives the following stability constraint where M = M g (r) is the gravitational mass within a radius r, which is defined by the Tolman-Whittaker mass formula Inserting Eq. (38) into (37), we get 2 where F g = − F ′ 2F (ρ + p r ) and F h = − dp r dr are the gravitational and the hydrostatic forces respectively, in addition to the anisotropic force F a . We note that the TOV equation should be modified in RT due to the non-minimal coupling constraint, T α β;α = ǫ ∂ β R, to include one more force F R as following where These different forces, for GR (ǫ = 0) and RT (ǫ 0), are plotted in Fig. 7 using the pulsar Her X-1 data. For the RT case, the figure clearly shows that the pulsar is dominated by the negative gravitational force over the positive hydrostatics and anisotropic forces. Therefore, it holds the pulsar in a static equilibrium. In conclusion, we verify the stability of the model via TOV equation using the pulsar Her X-1 data.

Relativistic adiabatic indices
Another verification of the stable equilibrium configuration of a spherically symmetric object can be done via the adiabatic index, that is defined as the ratio of two specific heats and can be given as follows [109][110][111] For the general case of anisotropic spheroid fluid, it has been shown that the object is in a neutral equilibrium if its adiabatic index Γ = γ and in a stable equilibrium if Γ > γ [111], whereas (a) γ in cases ǫ = 0 and ǫ 0 (b) Γ r in cases ǫ = 0 and ǫ 0 (c) Γ t in cases ǫ = 0 and ǫ 0 Fig. 8 Plots of the Adiabatic indices γ, Γ r and Γ t , namely (43)- (45), versus the radius r using the constants constrained from Her X-1. For RT, the adiabatic index γ less than the GR case but still greater than the neutral equilibrium value γ = 4 3 . The radial and tangential adiabatic indices have higher values whereas the stability constraints Γ r > γ and Γ t > γ are fulfilled everywhere within the pulsar.
In Figs. 8(a)-(c), we plot the profiles of the adiabatic indices γ, Γ r and Γ t of the pulsar Her X-1, respectively. It is clear that the value of the adiabatic index γ as predicted by RT is less than the corresponding GR value everywhere within the pulsar. On the other hand, the other adiabatic indices Γ r and Γ t are higher in RT than GR, whereas the stability constraints, namely Γ r > γ and Γ t > γ, as shown in the figures are clearly fulfilled.

More Observational Constraints
In this section we confront the present model with other pulsars' data to examine its validity with wide range of astrophysical observations. Also, we plot the mass-radius profile corresponding to different choices of the surface density compatible with the nuclear saturation density as boundary conditions showing the capability of the model to predict masses in the lower mass gap 2.5 − 5M ⊙ . On the other hand, we confront the compactness parameter as predicted by the model with Buchdahl compactness bound.

Pulsars' data
In addition to Her X-1, a similar analysis is developed for other twenty one pulsars cover a wide range of masses from 0.8M ⊙ to heavy pulsars of mass 2.01M ⊙ . In Table 1, we give the observed masses and radii for each associated with the corresponding model parameters {a 0 , a 1 , a 2 } assuming that the Rastall parameter ǫ = −0.1. We confirm that the model at hand can predict masses of those pulsars compatible with their observed values. In Table 2, we give the calculated values of physical quantities of most interest. As seen the values of the density are consistent with the nuclear density. It is to be mentioned that there are no EoSs are assumed, the obtained results however fit well with the linear behaviour, whereas the slopes dp r /dρ and dp t /dρ of the best fit lines are consistent with the values as obtained in Table 2. On the other hand these values confirm that the model is stable and causal for each case and  Table 2 Calculated physical quantities of the most interest.

Mass-Radius Profile
As is shown in Table 2 the surface densities of the listed pulsars, 2.14 × 10 14 ρ R 1.81 × 10 15 g/cm 3 , are mostly compatible with a neutron core. For four different values of the surface density of the pulsars ρ R = 2.7 × 10 14 g/cm 3 , 4 × 10 14 g/cm 3 , 6 × 10 14 g/cm 3 and 8 × 10 14 g/cm 3 we plot the corresponding compactness-radius curve as in Fig. 9(a). In all cases the maximum compactness values do not exceed unity. However for a compact object to be stable it should satisfy Buchdahl compactness bound U = 2G N M c 2 R ≤ 8/9 (for isotropic sphere [102]). We visualize Buchdahl upper bound on the compactness parameter with the corresponding maximum radii as obtained for the four surface densities. It is convenient to give the model parameters {a 0 , a 1 , a 2 } in terms of the total compactness parameter U. Recalling the matching conditions (35) and (36) we write , Notably for the pulsar Her X-1 with ǫ = −0.1, the compactness parameter U = 0.31, which reproduces the set of constants a 0 ≈ 0.369, a 1 ≈ −0.622 and a 2 ≈ 0.298 as previously obtained in Sec. 6. Using Buchdahl limit U = 8/9, we obtain the maximum radii R = 17.66 km for ρ R = 2.7 × 10 14 g/cm 3 , R = 14.51 km for ρ R = 4 × 10 14 g/cm 3 , R = 11.85 km for ρ R = 6 × 10 14 g/cm 3 and R = 10.26 km for ρ R = 8 × 10 14 g/cm 3 . However, more restricted values on the maximum compactness can be imposed by utilizing the SEC as obtained in the present model. Thus we write the density, radial and tangential pressures, namely (21), in terms of the compactness parameter U, after straight forward calculations we find that the SEC is satisfied when the compactness U ≤ 0.603 (using ǫ = −0.1). Higher U values can be obtained for higher values of ǫ, e.g. U(ǫ = −0.05) 0.639. However, we restrict our discussion to ǫ = −0.1. We visualize the SEC constraint on the upper compactness value as in Fig. 9(a). It can be seen that the lower U value as obtained by the SEC gives almost the same maximum radii as Buchdahl compactness bound or slightly lesser. However, it puts more restricted values on the maximum masses as will be discussed shortly. For a compact object to be physically stable we use U ≤ 0.603 as obtained by the SEC. This determines the maximum radii as allowed by the present model as R = 17.30 km for ρ R = 2.7 × 10 14 g/cm 3 , R = 14.22 km for ρ R = 4 × 10 14 g/cm 3 , R = 11.61 km for ρ R = 6 × 10 14 g/cm 3 and R = 10.05 km for ρ R = 8 × 10 14 g/cm 3 . In Fig. 9(b), we plot the mass-radius relation as obtained by the present model for the four surface densities along with the scatter plot of the observed mass-radius of the pulsars listed in Table 1. Consequently, we determine the maximum masses allowed by the present model according to Buchdahl limit for the surface densities ρ R = 2.7 × 10 14 g/cm 3 , 4 × 10 14 g/cm 3 , 6 × 10 14 g/cm 3 and 8 × 10 14 g/cm 3 , which have been found to be 5.31M ⊙ , 4.37M ⊙ , 3.56M ⊙ and 3.09M ⊙ respectively. In comparison with the GR case as obtained by Das et al. [96], RT allows more masses. Indeed these values provide unstable models, since they do not satisfy the SEC. We use the maximum compactness as obtained by the SEC, U = 0.603, which determines the physically stable pulsars with maximum masses 3.53M ⊙ , 2.90M ⊙ , 2.37M ⊙ and 2.05M ⊙ . Obviously from Fig. 9(b) all pulsars are below the SEC exclusion limit.
In Fig. 9(c), we split the pulsars into four groups according their matching with the mass-radius profiles of the surface densities. Notably, the nuclear saturation density ρ sat = 2.7 × 10 14 g/cm 3 and solidification density ρ solid is expected to be from 2.8 × 10 14 g/cm 3 to 5 × 10 14 g/cm 3 whereas the maximum density of cold matter ρ max 4 × 10 15 g/cm 3 based on an EoS-independent analytic solutions of Einstein field equations. In the present work, we find that the pulsar X7 has a neutron core with density less than ρ sat which matches anisotropic superfluid nuclear density ∼ 1.5 × 10 14 g/cm 3 . The pulsars PSR J0030+0451, PSR J0037-4715 and the LIGO constraints have a neutron core with density matches ρ sat . The pulsars M13, KS 1731-260, EXO 1745-268, GW170817-2, GW170817-1, 4U 1820-30, 4U 1724-207 and SAX J1748.9-2021 have a neutron core with density ∼ ρ solid . Notably, the most heaviest pulsars PSR J1614-2230 and PSR J0348+0432 fit well with the mass-radius curve corresponding to surface density 4 × 10 14 g/cm 3 which suggests that these pulsars are having solidified neutron cores. Similarly, the pulsars Her X-1, LMC X-4, EXO 1785-248, Cen X-3, 4U 1608-52 and Vela X-1 lie beyond neutron core density and match perfectly with mass-radius curve with a surface density of ρ R = 8 × 10 14 g/cm 3 , this suggests that those pulsars may have quark-gluon cores. This claim should be verified by gravitational wave signals and by experimental tests. It is to be mentioned that the four categories presented here fit well with linear EoS as is shown earlier in Fig. 5 whereas their slopes are given in Table 2.
Interestingly, the obtained curves allow for a heavy neutron star to lie in the lower mass gap 2.5 − 5M ⊙ which keeps an open window for the binary secondary companion GW190814 [119] to be a neutron star with an anisotropic core whereas the EoS p r (ρ) and p t (ρ) are obtained by simple linear relations. Clearly RT (with ǫ < 0), using Buchdahl limit, allows for higher maximum masses than the corresponding results as obtained by Das et al. within the GR framework [96]. This result is consistent with our earlier conclusion after Fig. 6(a) and also in agreement with the previous study of isotropic neutron stars in RT by Oliveira et al. [93]. However, in the present work we investigate anisotropic compact objects without assuming any EoS. In the following section we summarize and conclude the present work.

Summary and conclusions
One way to modify GR is to drop one of its main assumptions. Rastall questioned the local conservation in presence of curved spacetime, which led him to a set of field equations consistent with a different assumption, that is a non-minimal coupling between matter and geometry whereas the energy-momentum tensor T α β;α ∝ R ;β . Vissar, recently, claimed that RT is completely equivalent to GR [86]. On the contrary, Darabi et al. investigated Visser's claim but they concluded that Visser misinterpreted the matter-geometry coupling term, and consequently led him to a wrong conclusion [87]. In addition, they showed that by applying Visser's approach to f (R) theory one may conclude that it is equivalent to GR as well, which is not true. One of the good examples which may reveal the contribution of the matter-geometry coupling in RT in contrast to GR is the stellar models when the presence of matter play a crucial role. In this study, we applied RT to anisotropic static star with spherical symmetry seeking for a non-trivial interior solution. We showed that, for static spherically symmetric stellar model, the Rastall parameter in general does not contribute to the anisotropy parameter while it contributes non-trivially to the density and the radial/tangential pressures. Applying appropriate ansatz we derived the corresponding density and pressure forms.
We matched the interior solution to the exterior Schwarzschild one (which is a solution in RT as well) in the case of the pulsar Her X-1 (M = 0.85 ± 0.15 M ⊙ , R = 8.1 ± 0.41 km) and by recalling Zeldovich constraint, we determine the numerical values of the set of constants {a 0 , a 1 and a 2 } for a particular choice of Rastall parameter which has been found as ǫ = −0.1. We verified that the obtained model satisfies the physical constraints as listed in Sec. 4 and their predictions are slightly different from GR ones. We give detailed patterns of all physical quantities in Figs. 1-8. In addition, we confronted the model at hand with other pulsars' observations and we confirm that all physical conditions are fulfilled in all cases. In Table 1, we give the numerical values of the set of constants for the same Rastall parameter value ǫ = −0.1 of twenty additional pulsars as is done in the case of the pulsar Her X-1 . Also, in Table 2, we give the corresponding physical quantities of the most interest for all pulsar at the present study. All results confirm that the model satisfies the physical conditions as listed in Sec. 4.
For four surface densities, we utilized the SEC to set an upper limit of the compactness U ∼ 0.603 (where ǫ = −0.1) less than Buchdahl compactness U = 8/9, see Fig. 9. Consequently we determine the maximum possible masses with a stable configuration. For a boundary density compatible with neutron star core at saturation nuclear density the maximum mass ∼ 3.53M ⊙ . This results have been obtained without assuming a specific EoS, however, the pressure-density relation fits well with linear behaviour whereas the slopes keep a stable and causal structures. We split the twenty pulsars into four groups, three groups fit well with neutron cores with anisotropic superfluid, saturated and solidified cores. We found that the most massive pulsars PSR J1614-2230 and PSR J0348+0432 fit will with the mass-radius curve corresponding to surface density 4 × 10 14 g/cm 3 which suggests that those pulsars are having solidified neutron cores. Furthermore, we suggest that the following group of pulsars to have a quark-gluon cores, Her X-1, LMC X-4, EXO 1785-248, Cen X-3, 4U 1608-52 and Vela X-1. This group fit perfectly with a boundary density 8 × 10 14 g/cm 3 lies beyond the neutron core and suggested to have quark-gluon core. Even for those of relatively small masses we find that corresponding radii are much smaller than what is expected for a typical size of neutron star 10 km, which characterizes low-mass quark star [128,129]. In our case as we showed that the additional Rastall force in TOV equation contributes to resize the compact object to be smaller than GR prediction for the same mass. This in return results in higher core density ∼ 10 15 g/cm 3 . Noting that the nuclear saturation density ρ nuc ≈ 2.7 − 2.8 × 10 14 g/cm 3 , while the typical solidification density of neutron matter ρ solid 4 × 10 14 g/cm 3 .
To conclude: We investigated Rastall gravity, for an anisotropic star with a static spherical symmetry, whereas the mattergeometry coupling as assumed in RT is expected to play a crucial role differentiating RT predictions from GR ones. Indeed, all results confirm that RT is not equivalent to GR. However, Rastall parameter does not contribute to the anisotropy parameter, in general, i.e. RT predicts same amount of anisotropy as GR in the static spherically stellar models. Interestingly RT predicts more mass relative GR within same size which reflects the role of matter-geometry coupling to allow masses in the lower mass gap. This keeps the opportunity that the binary secondary companion GW190814 to be a neutron star with an anisotropic core whereas the EoS are obtained by simple linear relations.