Neutron Stars in Palatini $R+\alpha R^2$ and $R+\alpha R^2+\beta Q$ Theories

We study solutions of the stellar structure equations for spherically symmetric objects in Palatini $f(R)=R+\alpha R^2$ and $f(R,Q)=R+\alpha R^2+\beta Q$ in the mass-radius region associated to neutron stars. We illustrate the potential impact of the $R^2$ and $Q$ terms by studying a range of viable values of $\alpha$ and $\beta$. Similarly, we use different equations of state (SLy, FPS, HS(DD2) and HS(TMA)) as a simple way to account for the equation of state uncertainty. Our results show that for certain combinations of the $\alpha$ and $\beta$ parameters and equation of state, the effect of modifications of general relativity on the properties of stars is sizeable. Therefore, with increasing accuracy in the determination of the equation of state for neutron stars, astrophysical observations may serve as discriminators of modifications of General Relativity.


Introduction
General Relativity (GR) as the simplest realisation of a geometric description of gravity has so far shown, perhaps unexpectedly, total agreement with astrophysical and cosmological observations [1,2], as long as the general picture is accepted that the dark sectors, as needed for the ΛCDM model, are related to the particle and energy content of the Universe. Despite this tremendous success, there are reasons to investigate modifications of GR. First, GR breaks down in the high curvature regime. Second, it is only a classical field a e-mail: g.herzog@campus.unimib.it b e-mail: hsanchisalepuz@gmail.com theory without any quantum effects and which is not (perturbatively) renormalisable as a quantum field theory unless quadratic curvature corrections are added to the Lagrangian [3]. Moreover, it is conceivable that all or part of the effects currently attributed to dark sectors are in reality a manifestation of different gravitational dynamics [4][5][6][7].
Some of the (theoretical) problems of GR can be ameliorated in modified gravitational theories, like f (R) and f (R, Q) theories in which the Einstein-Hilbert Lagrangian is modified by adding polynomial terms in the Ricci scalar R or in the contraction of the Ricci curvature Q = R µν R µν (see e.g. [8][9][10][11][12][13] and references therein). From a different point of view, if the different terms in a gravity action are considered as quantum operators in the spirit of effective quantum field theories, studying their renormalisation group evolution indicates that f (R) and f (R, Q) terms render the theory asymptotically safe (and hence viable as a quantum theory) in the Planckian and super-Planckian regime [14][15][16][17][18]. Finally, classical theories of modified gravity can also serve as effective descriptions of their more fundamental quantum counterparts (see e.g. [19,20]).
Studying the compatibility of modified theories of gravity with phenomenology is a highly non-trivial issue since an interpretation of observations involves a combination of many different and uncertain physical mechanisms. As mentioned above, cosmological observations compatible with GR and dark sectors can be also made compatible with certain modifications of gravity and different contributions from dark sectors [21]. Inflationary scenarios can also be expressed in terms of modified theories of gravity (see e.g. [22][23][24][25][26]). It has been suggested that modifications of GR can be identified via the detection of new modes in gravitational waves [27][28][29][30][31][32][33][34]. Another possibility to investigate deviations from GR is with astrophysical observations. Different gravitational dynamics reflects into different stellar equations of structure [2, [35][36][37][38][39] and hence, in particular, in mass to radius ratios arXiv:2102.05722v3 [gr-qc] 19 Oct 2021 of neutron stars different to those expected from GR (see e.g. [40][41][42][43][44][45][46][47][48][49][50][51][52]). This is particularly interesting with the growing volume of observational data on neutron stars and the discovery of compact objects that are difficult to accomodate in a standard GR picture (see e.g. [53][54][55][56]). These measurements are, unfortunately, riddled with enormous uncertainties coming from the equation of state (EoS) describing the interior of neutron stars and can as well be interpreted as providing a handle on EoS assuming the theory of gravity is GR [57]. However, as we will show in this paper, the changes in the mass-radius relation for neutron stars coming from the uncertainty in the EoS are sometimes comparable to the differences generated by the modifications of the gravity theory, which implies that if, in the future, reliable ab-initio calculations of the EoS were available, astrophysical observations could be used as discriminators for gravity theories.
In this paper we shall investigate the influence of two particular models of f (R) and f (R, Q) theories, namely f (R) = R+αR 2 and f (R, Q) = R+αR 2 +β Q on the mass-radius relation for neutron stars. We study those theories in the Palatini formalism in which the dynamical degrees of freedom are the metric and the affine connection [11]. We will perform our calculations using a small number of representative EoS with the goal of estimating the uncertainty stemming from the EoS and from the theory of gravity (in our case simply encoded in the values of α and β ).
The paper is organised as follows. We sketch the derivation of the equations of stellar structure for Palatini theories in Sec. 2 and discuss several aspects of the EoS chosen in Sec. 3. Our results are shown in Sec. 4 and we conclude a discussion of the implications of our results on the validity of theories studied herein and of possible future work.

Stellar Structure Equations for f (R, Q) Palatini Theories
We discuss here the main aspects of the equations of stellar structure solved in this work. We consider only the simplified case of non-rotating and static stars, whose only observables are their mass and radius. The equations for f (R, Q) theories were derived in [36] and we refer to that paper for further details (see also [35,48]). In GR, the equivalent equations are the well-known Tolman-Oppenheimer-Volkoff (TOV) equations.
In Palatini f (R, Q) theories the action is given by where R is the Ricci scalar, Q = R µν R µν the aforementioned contraction of two Ricci tensors, with R µν = R ρ µρν and the Riemann tensor given by R d ab Γ d ec with Γ the connection coefficients, g µν is the metric, κ ≡ 8πG and S m [g, ψ m ] is the matter action. In the Palatini formalism the equations of motion are obtained after varying the action with respect to both, the metric and the connection One can introduce an auxiliary metric h µν via the (matter-content dependent) mapping √ −g ( f R g µν + 2 f Q R µν ) ≡ √ −hh µν such that the equation for the connection becomes ∇ β √ −hh µν = 0 and the connection Γ can thus be written as the Levi-Civita connection for h µν . Therefore, the Ricci tensor in Eqs. 2 and 3 can be written as the standard metric Ricci tensor of the metric h which we denote as R µν (h). Introducing Σ α ν := ( f R δ ν α + 2 f Q B α ν ), with B α ν given as B α ν = R αβ (h)g β ν , the relation between g µν and h µν is [36] The Ricci tensor is given by Note here that, when calculating R and Q from R(h) in f (R, Q), the Ricci tensor R µν (h) must be contracted with the physical metric g µν , that is, via the matrix B ν α above [36] In this paper we will only consider energy-momentum tensors of the form (7) T µν = (ρ + P) u µ u ν + Pg µν , and Lagrangians of the form f (R) = R + αR 2 and f (R, Q) = R + αR 2 + β Q. In these cases we can write the scalars R and Q in terms of ρ and P as R = −κT , as in GR, and [36] (8) The sign in front of the square root must be chosen so as to obtain the right limit at low curvatures.
The stellar structure equations are the dynamical equations for the star's interior metric, parametrised as , r 2 , r 2 sin 2 θ , where A(r) will be written as A(r) = 1 − 2M(r) r . After a tedious but straightforward derivation (see [36] for details) one arrives at the stellar structure equations for Palatini in the case of f (R, Q) = R + αR 2 + β Q, which read where, as before, a subscript denotes partial derivation, i.e. M r = ∂ M ∂ r and Ω rr = ∂ 2 Ω ∂ r 2 , and we introduced for compactness Here, Ω and S are given by the relation between the auxiliary metric h µν and the physical metric g µν where σ 1 and σ 2 appear upon rewriting the matrix Σ as Σ ν α = diag(σ 1 , σ 2 , σ 2 , σ 2 ) from Eq. 5 and are given by [36] where In order to recover the GR limit one has to use the positive sign in front of the square root in Eq. 12 and for σ 1 but the negative sign in λ [36]. Finally, τ refers to the right-hand side of Eq. 5, namely R ν µ ≡ τ ν µ . Note that, in Eq. 11, we have Ω rr = Ω PP P 2 r + Ω P P rr which involves the second derivative P rr , in contrast to the standard TOV equations. As it turns out, it is possible to write P rr in terms of the first order derivative P r . The result is [36] with α r = α P P r and β r = β P P r and (21) where we introduced Φ ≡ (τ r r + Ω S τ t t ). In this way, Eqs. 10-12 are a closed (after fixing an equation of state) system of equations expressed in terms of r, ρ(P), P, M and their first radial derivatives only.

Equations of State
The essential ingredient to solve the equations of stellar structure is an equation of state that relates the pressure to the energy density inside the star.
In oder to estimate the different effects on the mass and radius of neutron stars coming from the modifications of the gravity lagrangian and from the uncertainty in the EoS, we use a number of different EoS in analytic as well as tabulated form. Specifically, we use two different tabulated and two different analytic EoS. The analytic equations of state we used are the so-called SLy and FPS [58] The SLy and FPS EoS are analytic parametrisations of the results of manybody calculations with unified effective nuclear Hamiltonians [59,60], describing all regions of the neutron star interior from crust to core including its transitions. The analytic paramerisations take care that all thermodynamic conditions involving derivatives of the EoS are fulfilled, an aspect that will become problematic when using tabular data directly, as we will see. The SLy and FPS EoS are both parametrised by the following function where ξ = log(ρ/g cm −3 ) and ζ = log(P/dyn cm −2 ) and f 0 (x) is given by [58] f The fitted parameters a i in Eq. 22 are shown in Table 1. Additionally to the two analytic EoS above, we also studied neutron star solutions based on two other EoS in tabulated form. We choose EoS available in the CompOSE online service [61], namely the so called HS(DD2) [62,63] and HS(TMA) [62,64] EoS. These two EoS are both based on the the statistical model presented in [62], which includes the contribution of nuclei, nucleons, electrons, positrons and photons (excluding neutrino contributions) and which requires as input the masses and binding energies of nuclei and an effective model for the nucleon-nucleon interaction, which is then treated in the relativistic mean-field (RMF) approximation. The two EoS thus differ in the different model used for that dynamical input. The HS(DD2) equation of state uses the density-dependent nuclear model called DD2 [63], with the nuclei properties given by the FRDM model [65]. The HS(TMA) equation of state uses instead the TMA model [64], with the masses of nuclei taken from [66].
The drawback of using tabulated data to solve the equations of stellar structure is that interpolation functions and numerical derivatives thereof must be used. This procedure is ambiguous (as one can use many different forms for the interpolating functions) and, as mentioned above, certain thermodynamic relations must be preserved by the derivatives of the EoS (see e.g. [58]). This is a particularly difficult problem to tackle in the case of modified Palatini theories, since higher derivatives of the EoS are needed as compared to GR. Indeed, in the stellar structure equations for Palatini f (R) and f (R, Q) the first and second derivative of the density ρ with respect to the pressure P appear. For the analytically given EoS above this is relatively straightforward. Remembering that P(ρ) = 10 ζ (ξ (ρ)) = 10 ζ (log 10 (ρ)) we have and in order to arrive at the first and second derivative of the density with respect to the pressure we now have to invert Eqs. 24 and 25. One gets The result for the parametrization of the SLy and FPS EoS together with the first and second derivative of the parametrisations are shown in Fig. 1. For tabulated EoS data one can either apply a similar procedure as above after defining appropriate interpolating functions or simply calculate the derivatives numerically, as we did. When calculating the derivatives numerically, the different ways in which the data can be interpolated and hence derivatives can be taken has an impact on the results, as we discuss below.

Results
The mass and radius of a star are obtained by integrating from inside out the stellar structure equations after fixing a central density and until the pressure reaches a pre-determined threshold. In our calculations, the threshold was set to P = P c · 10 −12 where P c is the initial central pressure (determined by the EoS from the central density) in case of the analytic EoS. For tabulated EoS we stopped the integration when reaching the last entry of the EoS tables. We considered central densities in the range ξ = [14.4, 16] (at r = 0.01mm for numerical stability). The integration was performed using a fourth-order Runge-Kutta method with a fixed stepsize of 1m, since after exploring several adaptive-step methods we concluded that a fixed step was accurate enough and sped up the calculations.
A problem arising whenever a modified theory of gravity is used, is to determine the relative strength of the additional terms of the Lagrangian. Since in the present paper we use theories of the type f (R) = R + αR 2 and f (R, Q) = R + αR 2 + β Q, the problem reduces to finding reasonable values for α and β , consistent with phenomenology. However until now no solid experimental bounds on the values of α and β exist for theories in the Palatini formalism (see [67] for a study in Palatini). There exist, however, more restrictive experimental bounds on α for theories in the metric formalism (where only the metric is considered a dynamical degree of freedom). The Gravity Probe B experiment constrains α to α 5 · 10 15 cm 2 , while the constraint coming the Pulsar B in PSR J0737-3039 is four orders of magnitude higher [68]. In [68] they derived an even more stringent constraint on α from the Eöt-Wash experiment, constraining α to α 10 −6 cm 2 . However, this bound was derived in the low curvature regime of a laboratory on earth. Since possible modifications to GR will likely only be relevant in the high curvature regime, we do not take the bound from this experiment into consideration. Recently another bound for α has been derived by calculating the stability of stars in metric f (R) = R + αR 2 gravity [69]. Using polytropic EoS and looking for the maximum value of α which still fulfills the stability criteria they found α 2.4·10 8 cm 2 . To the best of our knowledge, no bounds on β exist at the moment.
Since the dynamics of metric and Palatini theories can be very different, the bounds discussed above may not apply in our case. Lacking a thorough analysis along the lines of [68,69] for Palatini theories, we used values for α and β for which an appreciable change in the mass-radius curve was obtained and that, for the case of α, are still within the bounds obtained in the metric formalism and in [67].

R + αR 2 case
We begin with the study of f (R) = R + αR 2 theories and show the mass-radius relation for a range of positive and negative values of α. Note that an analogous study was performed in [48] for the SLy and FPS EoS. However, our resuls are in contradiction with those presented in [48]. While in [48] the authors found only minor variations with respect to GR of the mass and radius of neutron stars, our calculations show instead a significant deviation from GR in some cases. We show in Fig. 2 our results for the mass-radius relation using the SLy and FPS EoS, compared to the corresponding GR result. As one can see, significant deviations from GR appear in the low-radius region even though smaller deviations are also observed for less compact stars. Interestingly, for negative values of α, f (R) = R + αR 2 theories generate heavier stars than GR and thus push the mass limit for the SLy and FPS EoS.  The situation is less pronounced in the case of the tabulated HS(DD2) EoS, Fig. 3. In this case, the differences in the mass-radius relation with respect to GR are still visible but are significantly smaller than in the previous case. What remains true is that negative values of α allow for a larger maximum mass. In the case of the HS(TMA) EoS, Fig. 3, the deviations with respect to GR are dramatic for the largest values of α we used. This is certainly an unexpected behaviour that we were able to trace back to a small difference in the HS(TMA) EoS with respect to the other EoS we used 1 , which induces a larger difference in the first and second derivative of the EoS. As can be seen in the first panel of Fig. 9 at around P = 10 32 g cm s the HS(TMA) EoS differs slightly from all other EoS but the difference is amplified for the first and second derivatives as can be seen in the second and third panel of Fig. 9. Note finally, that the range of radii obtained in the case of tabulated EoS is limited by the maximum and minimum values of the data provided, unlike in the case where we used analytic expressions.
In addition to the mass-radius plots, it is instructive to analyse the results as in Figs. 4-7, were we show the difference in percentage between the f (R) and the GR solutions, for the mass and the radius of a star separately. In this way we observe that the effect of the R 2 term on the radius of stars is generally smaller than on the mass, except for the lowest central-density regions. The shift in the mass is thus the main responsible of the change of the mass-radius ratios discussed above. The changes, moreover, are strongly dependent on the EoS; SLy, FPS and HS(TMA) exhibit large deviations with respect to GR, but the latter shows a qualitatively different behaviour as a function of the central density. For HS(DD2), the deviations with respect to GR are significantly smaller than for the other EoS, for both mass and radius.  It is clear from our discussion so far, that as long as the ambiguities in the EoS of neutron stars are not resolved, it is hardly possible to draw conclusions from mass and radius observations on the theory of gravity, as for a fixed f (R) (here, for a given α value) the mass-radius relation changes dramatically for different EoS. However, it is still interesting to speculate how relevant would astrophysical observations be, if we knew the correct EoS of a neutron star. For example, the observation of the 2.2M neutron star MSP J0740+6620 [70] has been used as an argument to rule out certain EoS, including some used in this paper. This is can be seen in Fig. 8, where we show the constraints put by MSP J0740+6620 on EoS, assuming GR is the theory of gravity and also using f (R) = R + αR 2 with α = −4 · 10 9 cm 2 (for which the increase of mass is the largest in our calculations). Clearly, some EoS can be rescued adding an R 2 term to the theory of gravity. Reversing the argument, with increasing certainty in the EoS of neutron stars, observations could be used to constraint deviations of the theory of gravity from GR.
At this point, we would like to discuss in some detail the differences in the results induced by the different methods used to calculate derivatives of the EoS in the case of tabulated data. To this end, we repeated the calculations with the analytic EoS SLy, this time generating tabulated data from the analytic expressions as follows. In one case, the derivatives were calculated analytically as outlined in Sec. 3 and separate data tables were generated from them; these values then had to be interpolated (in the same way as the EoS data is interpolated). As a second method we calculated the derivatives numerically with a finite difference approach, using second order central differences. In Fig. 10 the difference between the derivatives using those methods can be seen. While the difference for the first derivative is rather small and never exceeds 4%, the second derivative shows very large differences, of about 200% which rises up to 400%, at some points. To exemplify the impact of such differences in the mass-radius relations, we repeated the calculation for the SLy EoS with the different methods of calculating derivatives. The results are shown in Fig. 11. In the upper panel we show the calculation where the tabulated data for exact derivatives is interpolated, which clearly leads to the same results as when using the fully analytic version. However, in the lower panel, where the derivatives have been calculated with a finite difference method, shows only a minor difference between f (R) and GR, in fact in agreement with the results published in [48]. This discussion points to a potential problem appearing when calculating with tabulated EoS and theories with high-order derivatives of matter; all calculations that rely on numerical differentiation must be taken with due caution.
We proceed now to discuss the results for the f (R, Q) = R + αR 2 + β Q. As already indicated above, we consider val- Fig. 8 Maximum mass allowed by an f (R) = R + αR 2 theory with α = −4 · 10 9 cm 2 for the different EoS used in this paper, compared with the observed mass of the MSP J0740+6620 neutron star.   ues of β for which an appreciable change in the mass-radius curve is obtained. This means that, even the small changes in the stars masses that we observe, as discussed next, may be artificially large. We show in Figs. 12 and 13 the mass-radius relation for two values of β , with and without the R 2 term (here the R 2 term is excluded by setting α to a very small value, for computational convenience). The effect of the Q term appears to be generally small, the deviation from GR mostly coming from the R 2 term as we have seen in the previous section. The effect of the Q term seems to be more pronounced for analytic EoS than for the tabulated ones, which may point to the problem with numerical derivatives that we discussed above.
As in the f (R) case, it is interesting to plot separately the differences between the f (R, Q) and the GR solutions, for the mass and the radius of a star. We show them in Figs. 14-17. In contrast to the R 2 term, for which the dependence of the mass and radius as a function of the central density showed a very different qualitative behaviour for different EoS, in the case of the Q term these show a similar dependence on ξ c for all studied EoS. In all four cases the effect is somewhere below 5% on the mass, with the exception of HS(DD2) which is up to 8% for low central densities. As for the R 2 case, it is interesting to note that the Q terms influences the mass of the star much more than its radius, again with the exception of lowest values of the central density. In that region, the radius and mass effects compensate in a way such that they are not visible in the mass-radius plots.
To conclude, we wish to mention a peculiar feature of the solutions we obtained for f (R) gravities. In Fig. 18, we show a star profile (that is, the metric function M(r) as a function of the radial coordinate r) for the SLy EoS and ξ = 15.1. The metric function M(r) shows an unexpected cusp near the surface of the star. This feature, absent in GR, was already seen in [48] and predicted in [71] (see also [72]), where it was argued that it invalidates Palatini f (R) theories as physically viable. In [71] it was also speculated that a Q term may ameliorate this odd behaviour of the metric inside the star. As we see in Fig. 18, it is true that an f (Q) term (without R 2 ) does not show a cusp 2 . However, in an f (R, Q) model (i.e. with both R 2 and Q terms) the cusp remains; that is, the Q term is not able to compensate for the unexpected behaviour generated by the R 2 term. Note that in [71] the authors argue that for nonlinear f (R) one always finds divergences at the surface of the star, which may in turn generate curvature divergences. For our numerical results, we have carefully investigated the appeareance of divergences in all terms of the equations of structure (e.g. divergences in P r or P rr which may result, upon integration, in finite values of M(r)) but found none. From our calculations, thus, we can only conclude that the cusp in Fig. 18 is simply a dynamical effect with no clear relation to the arguments in [71]. Note finally that, even though in Fig. 18 we showed one particular solution only , in all cases we have studied in this paper the qualitative behaviour of the metric function M(r) was analogous.

Summary
In this paper we studied neutron stars solutions in Palatini f (R) = R + αR 2 and f (R, Q) = R + αR 2 + β Q gravities. We performed calculations using four different EoS and different values for α and β .
Our study shows that the differences in the masses and radius of stars induced by modifications of the gravity Lagrangian can in fact be as large as the uncertainty induced by the use of different EoS. Notably, it is the mass of the star and not its radius that is mostly affected by modifications of GR. We have also discussed separately the effect of the R 2 and the Q terms. The deviations from GR induced by the R 2 term are qualitatively different for different EoS. The Q term instead shows a similar behaviour for all EoS, and in general smaller than the effect from the R 2 term.
Even though it was not a goal of our investigations, we found that some of the EoS which have previously been ruled out by neutron star observations, might become viable again if the αR 2 correction is indeed part of the right theory of gravity. For negative α the maximum mass allowed by a given EoS tends to be increased. For some EoS, such as SLy or HS(TMA), it turns out that an αR 2 term can produce solutions as heavy as the MSP J0740+6620, the heaviest neutron star found so far [70].
To further highlight the relevance of studies like the present one it is interesting to mention the 2.6M secondary component of GW190814 [73]. Calculations using modified gravity theories indicate that it could be described as a neutron star (replacing MSP J0740+6620 as the heaviest known neutron star) Indeed, as shown in [49,50] for metric f (R) = R + αR 2 gravities, it is possible for certain EoSs to support neutron stars in the 2.6M mass range. Interestingly, the authors showed that studying rotating stars also increases the maximum mass of neutron stars for those modified gravities. In this paper we observe a similar effect for Palatini f (R) = R + αR 2 and f (R, Q) = R + αR 2 + β Q, obtaining neutron star solutions in the 2.6M mass range for the HS(TMA) EoS (see Fig. 3). Alternatively, allowing a sufficiently negative value for α also leads to neutron stars in this mass range when using other EoS, like the analytic SLy EoS.
In this work we focused on the relatively simple problem of calculating the mass, radius and profiles of spherically symmetric stars. In the future, it is necessary to perform more detailed studies. For example, it is crucial to study the stability of a solution in order to understand whether it is viable as a physical solution. For example, it could be that the solutions beyond the GR mass limit are unstable and hence not realised in Nature. This is an even more pressing issue, considering the odd behaviour of the metric function M(r) for f (R) theories, as discussed in last section. Also, realistic neutron stars rotate and the system of equations solved herein provides only an estimate of their properties. Finally, on the technical side, it is necessary to establish reliable and unambiguous methods to calculate high-order derivatives of EoS, as we have seen.
Clearly, no conclusive results can be obtained yet, at least until the ongoing efforts to nail down the EoS for a neutron star from first principles are successful. However, our study indicates that, in the future, astrophysical observations may serve as discriminators of modified gravity.