{'fraction_non_alphanumeric': 0.04260877979550366, 'fraction_numerical': 0.016224236364222818, 'mean_word_length': 4.488405027438485, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The mass-radius relation for neutron stars in f (R) = R + αR2 gravity: a comparison between purely metric and torsion formulations

P. Feola1, Xisco Jiménez Forteza3, S. Capozziello2,3,4, R. Cianci1, S. Vignolo1 1DIME Sez. Metodi e Modelli Matematici, Università di Genova, Via All’ Opera Pia 15a 16100 Genova (Italy) 2 Dipartimento di Fisica, “E. Pancini” Università “Federico II” di Napoli, Compl. Univ. Monte S. Angelo Ed. G, Via Cinthia, I-80126 Napoli (Italy) 3 INFN Sez. di Napoli, Compl. Univ. Monte S. Angelo Ed. G, Via Cinthia, I-80126 Napoli (Italy) and 4Laboratory for Theoretical Cosmology, Tomsk State University of Control Systems and Radioelectronics (TUSUR), 634050 Tomsk, Russia. (Dated: September 20, 2019)


I. INTRODUCTION
Compact objects, such as Neutron Stars (NS), are astrophysical objects that can be described by General Relativity (GR). These relativistic stars are natural laboratories for studying the behavior of high-density nuclear matter using an appropriate equation of state (EoS), which relates the pressure and density of degenerate matter. This allows to obtain the Mass-Radius relation, M − R, and other macroscopic properties such as the tidal deformability and the stellar momentum of inertia [1].
Since the internal structure of a NS cannot be reproduced in the laboratory because of the extreme conditions in which it operates, only theoretical models can be formulated where there are a very large number of EoS candidates. The astrophysical measurements of the macroscopic properties of NS are very useful because they allow us to understand what can be realistic EoS . In fact, they can provide information on whether the EoS is soft or stiff and what is the pressure several times the density of nuclear saturation [2][3][4][5]. Therefore, measuring the mass value of a NS could help us to describe matter at extreme gravity regimes.
Einstein's theory describe accurately the physical properties that govern the stability of NS where Chandrasekhar, considering degenerate matter, fixed a theoretical upper limit of 1.44M so that the stability of a non-rotating degenerate star is conserved [6]. Instead, as confirmed by several astrophysical observations, there exist binary systems with NS having mass values that violates this limit allowing larger masses [7][8][9][10][11][12].
To study these observational evidences, as already done in some previous works, developed in metric formalism, [13][14][15][16][17][18], Extended Theories of Gravity [19] can be used. In particular f (R) gravity, i.e. a class of Lagrangians considering a generic function of the Ricci curvature scalar. The primary objective is to obtain the M − R relation for a NS that allows, given an EoS, to derive the maximum mass value.
From a cosmological point of view, f (R) theories, beside addressing in a straightforward way the inflationary paradigm [20], could be useful in view of problems like the accelerated expansion of the universe (the dark energy issue), confirmed by several observations [21][22][23][24][25][26], and the problem of the formation of large-scale structures, called dark matter. Unlike the standard Concordance Lambda Cold Dark Matter (ΛCDM) Model [27][28][29], similar results can be obtained without considering dark components but extending the gravitational sector at infrared scales [19,[30][31][32][33][34][35][36][37]. Specifically, f (R) gravity is acquiring a growing interest because it allows a good description of gravitating structures without non-baryonic dark matter: extra degrees of freedom of gravitational field can be dealt as effective scalar fields contributing to the structure formation and stability [38,39]. In this perspective, it is possible to unify the cosmic acceleration [30,40], the early-time inflation [20,41], thus leading to a complete picture of the evolution of the Universe [32,[42][43][44][45][46][47] and large-scale structures therein [48][49][50][51]. However, the dark side and the f (R) descriptions are, in some sense, equivalent at large scale so one needs an experimentum crucis capable of discriminating among the two competing pictures. Discovering new particles out of the Standard Model or addressing gravitational phenomena that escape from the GR description could be an approach to fix this challenging issue. Observing exotic stars modeled by some alternative theory of gravity could be a goal in this perspective.
In this paper, we derive the M − R diagram, using realistic EoS compatible with the LIGO constraints [52] for a f (R) = R+αR 2 Lagrangian, using two different approaches: the purely metric theory and a theory with torsion that allows to introduce the spin degrees of freedom in GR [53]. In our specific model, the torsion field is due to the non-linearity of f (R). Here the massenergy is the source of curvature and the spin is the source of torsion. In this way, torsion contributions could provide additional information for compact stars in extreme gravity regimes. The goal of this paper is to obtain realistic M − R relation by solving numerically a modified system of equations, derived from Tolman-Oppenheimer-Volkoff (TOV) [54] equations, and compare results with the LIGO constraints. Specifically, we shall consider quadratic corrections to the Ricci scalar and discuss models with and without torsion comparing them with GR.
The paper is organized as follows. In Section II we derive the TOV equations for f (R) gravity in the metric and torsion formalism. Section III is devoted to the problems related to the numerical aspects of TOV equations in f (R) gravity. In Section IV we derive the numerical solutions of stellar structure equations and compare the results of the M − R relations. Discussion and conclusions are given in Sec. V.

A. The metric theory
In the metric formulation, the action for f (R) gravity (in units for G = c = 1) is given by where f (R) is a function of the scalar curvature R, g is determinant of the metric tensor g i j and L matter is the matter Lagrangian.
Varying the action (1) with respect to the metric tensor g i j , one gets the field equations: In eqs. (2), R i j is the Ricci tensor, f (R) denotes the derivative of f (R) with respect to the scalar curvature, δg i j is the energy-momentum tensor of matter and = 1 Here we adopt the signature (+, −, −, −).
In order to describe stellar objects, we assume that the metric is static and spherically symmetric of the form: where ψ and λ are functions depending only on the radial coordinate r. We assume that the interior of the star matter is described by a perfect fluid, with energy-momentum tensor Σ i j = diag(e 2ψ ρ, e 2λ p, r 2 p, r 2 p sin 2 θ), where ρ = ρ(r) and p = p(r) are the matter density and pressure respectively. By a direct calculation, it is possible to show that field eqs. (2), evaluated in the metric (3), are equivalent to a set of equations consisting of the Tolmann-Oppenheimer-Volkov (TOV) equations for f (R) gravity and a continuity equation given by the contracted Bianchi identity ∇ i Σ i j = 0. Specifically, the TOV equations for f (R)gravity are while the continuity equation is Here R r and R r,r denote respectively the first and second derivative of R(r) with respect to radial coordinate r. In order to solve numerically the equations (4), (5) and (6), we can consider the scalar curvature R as an independent dynamical field. In doing this, we need an additional equation which is directly obtained form the definition of scalar curvature: Indeed, inserting the content of eqs. (4), (5) and (6) into (7), we get the dynamical equation for R: Finally, the numerical solution of the resulting dynamical equations relies on the assignment of a suitable EoS, p = p(ρ), relating pressure and density inside the star, as well as of initial data (i.e. values of the fields at the center of the star).

B. The theory with torsion
In f (R) gravity with torsion, the gravitational and dynamical fields are pairs (g, Γ) consisting of a pseudo-Riemannian metric g and a metric compatible linear connection Γ. with non-vanishing torsion.
The corresponding field equations are obtained by varying the action functional (1) independently with respect to the metric and the connection. It is worth noticing that now R refers to the scalar curvature associated with the dynamical connection Γ.
Moreover, we recall that any metric compatible linear connection Γ may be decomposed as the sum whereΓ h i j is the Levi-Civita connection associated with the given metric g and K h i j denotes the contorsion tensor, related to the torsion tensor T h i j = Γ h i j − Γ h ji by the relation [55]: The contorsion tensor (10) verifies the antisymmetry property K j h i = −K h j i and, together with the metric tensor g, identifies the actual degrees of freedom of the theory.
Making use of eqs. (9) and (10), we can decompose the Ricci and the scalar curvature of the dynamical connection respectively as: and whereR i j andR are the Ricci and the scalar curvature of the Levi-Civita connection induced by the metric g.
In the absence of matter spin density, variations of (1) yield the field equations [56][57][58][59][60]: and where Σ i j denotes again the energy-momentum tensor of matter, and the non-linearity of the gravitational Lagrangian function f (R) becomes a source of torsion. Now, by inserting eqs. (11) and (14) into eqs. (13), it is possible to show that the whole set of field equations evaluated in the metric (3) is equivalent to the system formed by the following two TOV equations together with the continuity equation which also holds in the present case [61,62]. Also in the torsional case, we consider the scalar curvature R as an independent dynamical variable, introducing a consequent additional equation derived from the very definition of R itself. In fact, inserting eqs. (10) and (14) into (12), evaluating all in the metric (3) and making use of eqs. (15) and (16), we obtain the evolution equation: Again, in order to be solved, the set of dynamical equations (15), (16), (17) and (18) for the unknowns R, λ, ψ, p and ρ must be completed by an EoS and initial data.
We consider here the specific form of f (R): where α is the coupling parameter of the quadratic curvature correction. This model is specially suitable to account for cosmological inflation, where higher-order curvature terms naturally lead to cosmic accelerated expansion. The quadratic term emerges in strong gravity regimes, while at Solar System scales and, more in general, in the weak field regime, the linear term predominates. Since the interior of a NS could present energy conditions in some sense similar to those early universe [15], the model (19) is particularly suitable for our considerations. In this model, eqs. (4), (5) and (8) take the explicit form: while eqs. (15), (16) and (18) become respectively: Clearly the torsion contributions emerge in the second system. In the next Section, we shall discuss numerical solutions for the interior space-time of spherically symmetric NS in both metric and torsional f (R) = R + αR 2 gravity. Our aim is to compare the solutions of the above two systems of differential equations in order to point out the torsion contribution with respect to the purely metric one.
In view of this, it is worth noticing that, in vacuo, f (R) = R + αR 2 gravity with torsion amounts to GR plus, eventually, a cosmological term [56,60]. Therefore, under the assumption of spherical symmetry, in the case with torsion, the spacetime outside the star has to coincide with the Schwarzschild one. In order to compare the two models, it is then reasonable and consistent assuming the external Schwarzschild solution as the space-time outside the star also in the case of purely metric theory. It is worth noticing that the Schwarzschild external solution is actually a vacuum solution for purely metric f (R) = R+αR 2 gravity as demonstrated in [63,64].
Therefore, viable interior solutions, at the boundary, have to match the external Schwarzschild solution. In this regard, we recall that junction conditions for f (R) gravity have been studied in [65] for the purely metric formulation, and in [66] for the theory with torsion. Referring the reader to [65,66] for more details, we assume the following junction conditions at the stellar radius Outside the star λ, ψ and R refer to the corresponding Schwarzschild quantities. Eqs. (26) and (27) are the conditions at the stellar radius to be satisfied by the numerical solutions we shall investigate in the next Sections.
The TOV equations presented in Sec. II, together with an EoS, form a closed system of equations that can be solved numerically once a suitable set of initial conditions are provided. The EoS accounts for the behavior of the matter fields in the NS at nuclear level. However, it also dominates the NS macroscopic properties as the total mass M, radius R S and compactness C = M/R S . The total mass M and the radius R S may vary significantly depending on the state of matter in the NS interior where C ≈ [0.02, 0.25], being C = 0.5 the black hole solution. On the other hand, the knowledge of the macroscopic properties provides a direct insight to understand the particle interactions, energy transport and state of the matter in the NS core. Until recently, there were placed only vague constraints on the EoS of NSs from electromagnetic observations [67]. The recent LIGO-Virgo binary neutron star (BNS) observation has significantly clarified the state of art concerning the EoS physics. The largest accuracy of the gravitational wave (GW) channel in relation to the electromagnetic (EM) observations allowed to rule out stiffer solutions (less compact) thus reducing significantly the number of astrophysically relevant EoS. In this section, we discuss some aspects of the numerical solution of TOV equations in the metric and torsional f (R), formulations described above, for four EoS compatible with the recent LIGO constraints: APR4, MPA1, SLy, WFF1 [68][69][70][71], accurately described the piecewise polytropic fits provided in [72].
Then, to solve numerically the TOV equations, we use a dimensionless version of the them by re-scaling our physical variables as where and M is the mass of the sun, r g is the gravitational radius ( 1.5km), G Newton's Gravitational constant and c the speed of light. The two systems of differential equations shown in Subsection II C take the following form, where the primed variables denote radial derivatives. Therefore, we are left to setup five initial conditions (ICs) for the variables {p(0), λ(0), ψ(0), R(0), R (0)} to complete the numerical scheme. ICs are chosen at the center of the star r = 0 in order to preserve regularity, thus preventing the generation of large gradients that may lead to numerical instabilities. Mathematically, this involves that any expansion around the NS center must have a zero first derivative. In particular, the scalar curvature at the NS center may be expanded as, where regularity involves R (0) = 0. Pressure and density at the center ρ(0) = ρ c and p(0) = p c are given by the EoS so they only depend on the type of fluid under consideration. For the metric potential λ, it is natural to fix λ(0) = 0, analogously to what happens in Newtonian gravity, where the λ(r) and ψ(r) variables are matched to the m(r) mass of the system by, Notice that the variable ψ(r) does not enter directly in our system of differential equations which implies that ψ(0) can be defined up to any arbitrary constant. Therefore we adjust ψ(0) conveniently to match (i) the internal solutions with the external Schwarzschild solution at the stellar radius R S and (ii) to obtain asymptotically the O(r −1 ) profile as, The star radius is ideally defined where the pressure p(R S ) ≈ 0 though, in practice, and for numerical reasons, it is sufficient to set a ground value as p(R S )/p c ≤ ∼ 10 −10 .
The fulfillment of eqs. (33) requires to find an optimal choice for the Ricci scalar R c = R(0). In general, this is achieved by shooting the central value R c within some sufficiently large range [R min c , R max c ], containing the true value R c . Then R c is found by applying bijection root-finding methods until eqs. (33) are satisfied up to numerical tolerance. Unfortunately, the existence of such R c strongly depends on the particular form of the f (R) model, giving rise to ghosts in case of ill-defined configuration of the model parameters. This is true for both metric and torsional (R + αR 2 ) theories discussed in this work. Then we choose the sign of α to be the one that better matches the junction conditions at the surface of the star (26), (27) for the metric and torsional theory respectively. As we evince in the following sections, the only choices that reproduce not blowing up solutions are α > 0 for the metric case and α < 0 for the torsion one. Unfortunately, these choices generate some typical tachyonic oscillations due to a bad behaved f (R) and that we could not remove numerically. This effect was also reported in [73] and it shows an oscillatory behavior, in the form of a damped-sinusoid outside the star, even in the minimally perturbed scenario with α 1. These oscillations grow as the value of α increases and they are as well propagated to our metric potentials λ(r) and ψ(r). This inserts some ambiguity in defining the asymptotic conditions (33) for large r since the oscillations are not totally vanished when the numerical noise begins to dominate the solution (for r ∼ 100). To overcome this issue and to reduce the amplitude of the oscillations, we restrict our analysis to small values of α ∈ [0.001, 0.1]. As this is anyway consistent with current observational tests, in doing so, we are not discarding any relevant astrophysical scenario and this fact allows us to set R c ≈ R GR . This hypothesis is shown to have a minimal impact in the M − R diagrams as we will discuss throughout next sections. Moreover, the assumption of a Schwarzschild-type solution outside the star allows us to smooth out these oscillations and to recover a good fulfillment of the junction conditions. According to the above positions, we justify the choice of α > 0 for the metric theory and α < 0 for the torsional one.
Finally, the two systems of ordinary differential equations (ODE) are solved by using a 8th-order Runge-Kutta with adaptive step-size and high-stiffness control methods implemented in the Wolfram Mathematica package [74]. These methods regulate the discretization step-size by estimating the error of the Runge-Kutta method point by point ensuring the numerical convergence of the solution step by step. The stiffness control methods use polynomial extrapolation on the short regimes where the gradients become too large. We have found these methods essential to ensure the accuracy of the solutions in the torsional formulation.

IV. NUMERICAL SOLUTIONS
We compute the M − R diagrams for metric and torsional formulations of f (R) = R + αR 2 gravity. Due to the numerical limitations found throghout our analysis, we restrict |α| ∈ [0, 0.1] where α is required to be positive for the purely metric theory and negative in the theory with torsion to avoid blowing up solutions [73]. These values are anyway consistent with solar system tests of GR [73,75]. Such tests fix light constraints on the form f (R) 10 −6 rather than on the parameter α, thus being translated as R + |α|R 2 10 −6 . Bearing in mind that curvatures themselves are expected to be small, this leaves the parameter α rather unconstrained. Other tests as Eöt-Wash laboratory experiment set α 10 −10 m 2 . On the contrary, there exist alternative observational space-based constraints coming from the Gravity Probe B experiment [76] or the observation of the binary pulsar PSR J0737-3039 [77,78] that set α [5 × 10 11 , 2.3 × 10 15 ]m 2 . Therefore, the discrepancies among the several experiments do not set tight bounds on the value of α, and our choice seems to be compatible with existing data.

A. Purely metric theory
The solutions of the TOV equations for the purely metric f (R) = R + αR 2 model are illustrated in Figure 1. The pressure at the center of the star p c drops quickly until it eventually gets equal to zero, thus defining the radius of the star R S . This radius is used as our reference point to compute the total mass M by means of eq. (32). The numerical system exhibits some dissipative oscillations about the Ricci scalar R and the metric potential λ. These oscillations naturally arise from the harmonic-form of the Ricci scalar R(r) equation in vacuum [73], for a non optimal choice of the Ricci scalar R c at the center of the star, and where optimal choice is here defined as that matching the Schwarzschild junction conditions at the stellar radius. Unfortunately, such a choice becomes increasingly difficult as α tends to zero since the system of equations become also stiffer [79]. Generally speaking, this may appear to be counterintuitive, since α → 0 should exactly recover the GR space-time. However the asymptotic approach to α → 0 of the Ricci scalar equations (22) (25) are ill-defined. This is clear if, for instance, one re-expresses (22) as, Notice that the numerator of the first term is exactly zero in GR and that ideally approaches to zero faster than linear order in α. However, this is not so exact when dealing with numerical uncertainties, where the same factor may behave as a ∼ 0/0 solution for α << 1 thus requiring much more precision on the estimation of central value R c . To overcome this issue, we have set R(0) = R GR = 8π(3p c − ρ c ) to the GR value. Though this seems apparently an arbitrary choice, we notice that, for α < 1, the solution must be close to t GR so the value cannot be further to that of GR. This is self-evident from Fig. 2, where, in the right plot, we illustrate the variations on the pressure p(r) and the Ricci scalar R(r) for different choices of the central value R c = R GR c , 0.2R GR c , 2R GR c . Then, notice than the effect of varying R c on the radius R for such small values of α is about ∼ 2% considering the maximum and minimum choices of R c . This variation is then compared with the uncertainty arising from the definition of the star radius R S to be the place where the pressure drops by a factor . Then, in the left plot, we show that the impact of relaxing this value to ∼ 10 −9 would generate an uncertainty of about 4%, thus larger than the one from varying R c . ; a Schwarzschild fit (green line) to the numerical data outside the star, that is, with R > 11.6km. We note that, for α smal and averaging out all the oscillations, all physical quantities reproduce rather well the Schwarzschild solution outside the star, while matching as well the junction conditions (26). From the fitted results we get M = 1.40M , thus very close to the theoretical one.

Exact
In Fig. (3), we show the behavior of the metric potentials λ(r) and ψ(r) and the derivatives R (r) and ψ (r) paying special attention to: (i) the junction conditions at the NS boundary and (ii) their profiles as r → ∞. We show the full numerical solution (blue line), its corresponding Schwarzschild solution (orange line) given by eqs. (32) with M = 1.43M and the result of fitting the exterior data to the same Schwarzschild-like ansatz in order to quantify the agreement with the Schwarzschild solution outside the star and which results in a NS with total mass M = 1.40M . The good agreement between the three lines confirms that the solution is well approximated by the Schwarzschild solution right outside the star radius better than ∼ 2%. This good match is also extended to their derivatives thus globally satisfying the necessary junction conditions of eqs. (26) once the oscillations are averaged out. On the other hand, since the oscillations do not appear on ψ(r), we choose this quantity more appropriated to define the NS mass M.
Finally, in Fig. 4, we show the M − R diagrams for the four EoS considered in this work. For each choice of the central density ρ c , we get a different estimate of the radius R S and the total mass M. We loop over ρ c until dM/dR = 0 which defines the unstable branch, i.e. the point at which the NS is expected to collapse to a black hole and that provides the maximum allowed mass M max for the given EoS. Note that for all the EoS considered, the total mass tends to increase with respect to GR as in [18,79]. This is because gravity becomes stronger, thus allowing more massive systems. Indeed, in the f (R) = R + αR 2 scenario, Newton's gravitational constant G is replaced by The combined conditions of α > 0 and R < 0 imply then G e f f > G, thus generating a more attractive gravity.

B. Theory with torsion
We repeat the analysis for the torsional f (R) = R + αR 2 theory. Although further models have been also considered in the literature, the numerical complexity of torsional equations makes difficult a full exploration of other kinds of f (R) functions. This issue becomes more relevant when considering the torsional theory with spin [57], where spin gradients add higher order derivatives to our system of equations that increase the stiffness of the numerical system. We plan to extend our study in the presence of spin matter in a forthcoming paper. Then, in Figure 5, we show the results we obtained for the theory with torsion, using the same range for |α| as in the metric case but choosing α < 0. In this scenario, we see that the general trend predicts a decreasing of the total mass of NS, independently of the EoS considered. This could be related with the fact that the stable branch of the solutions, given by the sign of α, is reversed with respect to the purely metric case to avoid for ghosts. However, estimates for the total mass and radius are still compatible with the astrophysical observations [4], thus not allowing us to rule out any of the models studied here. On the other hand, if we further increase |α| the errors generated by eq. (25) and propagated to the total mass M and the total radius R S become too large. Therefore, we restrict our analysis to |α| ≤ 0.1. In Fig. 6 we repeat the same Schwarzschild-based tests adopted for the metric formalism for α = 0.05. In this case, the total mass M = 1.37M is slightly diminished with respect to the metric case. Notice that the Schwarzschild solution is as well verified at the star radius R s , where the metric λ(r) is clearly C 0 and ψ(r) still preserves the C 1 condition. Outside the star, and once the oscillations are vanished, the metric functions λ and ψ still preserve the 1/r decay.  Finally, in Fig 7 we compare the different predictions obtained in the purely metric and the torsional formulation respectively, for α = 0.1. Note that, in the theory with torsion, though the total mass of the NS decreases, while increases with respect to the metric case, the relative deviations, in absolute value, with respect to GR seem to be larger than in the metric case. This is caused by the effective repulsion generated by the extra torsional terms (see eq. (25)) which induce a partial screening of the gravitational field that prevents to reach NS masses as large as in standard GR. This is explicitly shown in Table I, where we show the variation of the maximum mass M max , radius R max and compactness C for the purely metric and torsional theories respectively, corresponding to the points in the M − R diagrams where dM/dR = 0. Note that whereas the purely metric formulation tends to more massive and compact NSs, the opposite occurs when considering torsion. Specifically, as the quadratic term in the curvature increases, the effects of torsion counterbalance the increase of total mass. This can be intuitively derived by using the same reasoning as in (35) with α < 0 and R < 0, implying G e f f < G and thus generating a less attractive gravity.

V. DISCUSSION AND CONCLUSIONS
In this paper, we have studied the existence of realistic NSs in the context of the f (R) = R + αR 2 theory both in purely metric and torsional formulations. The main results concern the computation of the M − R diagrams resulting from the two different theoretical frameworks considered. Matter fields have been represented by static and spherically symmetric perfect fluids where the EoS have been chosen to agree with the recent LIGO-Virgo constraints [52]. The parameter α has been restricted to be smaller than |α| ≤ 0.1 to avoid unrealistically large oscillations (see e.g. [73]) on our metric potentials and therefore ensuring the(i) fullfillment of junction conditions and (ii) the accurate recovery of the Schwarzschild solution far from the source. These two requirements single out four of the five initial conditions: p(0), λ(0),ψ(0) and R (0), while R(0) remains free. R(0) is ideally defined by choosing this parameter in such a way to match the junction conditions (26), (27). However, the oscillatory behavior of some solutions for r → ∞ prevents from finding a unique value for R(0). To overcome this issue, we have set R(0) = R GR identical to the GR value. This assumption have been shown to be valid for small α, being the estimates of the NS radius only mildly dependent on the R(0) choice, but this is no longer true for α 1. In the purely metric theory, the obtained results show a progressive increasing of the total mass as |α| increases, for all the four EoS considered. This allows for higher masses and more compact NSs than in GR. This absloute increasing of the mass and compactness could be also reproduced by assuming softer EoS in GR, consistent with the recent observations [52]. In the case with torsion, the NS mass tends to decrease for all the EoS considered. This could be related with the fact that the stable branch of the solutions is flipped with respect to the purely metric case to ensure the stability of the numerical system. The physical existence of such solutions could help us to describe NS compact or not, based on astrophysical observations, choosing the appropriate theory by simply constraining whether α is positive or negative. In the torsional framework, the differences in the M − R predictions with respect to GR are larger than those obtained in the purely metric case. As a consequence, the allowed intervals on α are poles apart from the two theories. Moreover, the theory with torsion would seem to describe less compact NS. This would allow to obtain solutions that could be reproduced using EoS with stiff matter in the limit of GR. Unfortunately, this is in disagreement with the recent LIGO-Virgo discoveries [52]. What comes to the rescue is that given the current accuracy of electromagnetic observations, we cannot deny the NS observations yet because the differences with the GR are still too small. However, this issues could be addressed by next generation gravitational wave detectors (3G) [80][81][82], where the opportunity to test results presented in this work could be realistic.