Tidal Love numbers of neutron stars in f(R) gravity

The recent detection of gravitational waves from a neutron star merger was a significant step towards constraining the nuclear matter equation of state by using the tidal Love numbers (TLNs) of the merging neutron stars. Measuring or constraining the neutron star TLNs allows us in principle to exclude or constraint many equations of state. This approach, however, has the drawback that many modified theories of gravity could produce deviations from General Relativity similar to the deviations coming from the uncertainties in the equation of state. The first and the most natural step in resolving the mentioned problem is to quantify the effects on the TLNs from the modifications of General Relativity. With this motivation in mind, in the present paper we calculate the TLNs of (non-rotating) neutron stars in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R^2$$\end{document}R2-gravity. More precisely, by solving numerically the perturbation equations, we calculate explicitly the polar and the axial \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$l=2$$\end{document}l=2 TLNs for three characteristic realistic equations of state and compare the results to General Relativity. Our results show that while the polar TLNs are slightly influenced by the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R^2$$\end{document}R2 modification of General Relativity, the axial TLNs can be several times larger (in terms of the absolute value) compared to the general relativistic case.


Introduction
The first detection of binary neutron star mergers [1] contributed in several ways to the efforts of constraining the nuclear matter equation of state (EOS) [2][3][4][5][6][7][8][9]. Perhaps one of the most elegant and straightforward constraint comes from the measurement of the tidal Love numbers (TLNs) of the merging neutron stars. On its basis one can already exclude a large number of modern equations of state. This approach, a e-mail: yazad@phys.uni-sofia.bg b e-mail: daniela.doneva@uni-tuebingen.de c e-mail: kostas.kokkotas@uni-tuebingen.de though, have the drawback that many alternative theories of gravity would produce deviations from pure general relativity (GR) similar in magnitude and characteristics to the uncertainties in the EOS. 1 Thus, it is difficult to determine the EOS from the current and forthcoming gravitational wave observations in a theory of gravity independent way. The first step in solving this problem is to quantify the effects from modification of general relativity on the TLNs.
The TLNs characterize the response (deformability) of a body to an external tidal force [17,18]. They encode information about the internal structure of the body and the strong field regime of gravity and most importantly -the tidal Love numbers can be determined through the gravitational wave emission of merging neutron stars [19][20][21][22][23][24][25][26][27][28]. TLNs in alternative theories of gravity were studied only in few cases [29]: the dynamical Chern-Simons (dCS) gravity [10,11], the Eddington-Inspired Born-Infeld (EiBI) gravity [30] and in massless scalar-tensor theories (STT) of gravity [31]. One should note that in the dCS gravity the polar TLNs, that give the dominant contribution to the gravitational wave signal, are the same as in GR up to a leading order and for EiBI gravity the standard GR framework for calculating the TLNs, using an apparent EOS formulation, can be employed. On the other hand in [31] the perturbation equations for general STT were derived and the TLNs were calculated for the case of massless scalar field. The TLNs of different exotic compact objects and for BHs in some alternative theories of gravity (having nonzero TLNs in contrast to pure GR) were examined in [32].
The interest in different modifications of GR on the other hand is growing in the last decade. The reason steams from the fact that from one hand there are phenomena such as the accelerated expansion of the Universe, that can not be explained well within Einstein's theory of gravity without requiring fine-tuning or other problems appearing. From the other hand there are purely theoretical arguments coming from the theories trying to unify all the interactions, the quantum corrections to the GR Lagrangian, or the attempt to quantize gravity [33].
Thus a natural step is to calculate the TLN of neutron stars in a broader class of alternative theories of gravity. A significant complication comes from the fact that the exterior solution of the perturbations equations is not known analytically in the general case, in contrast with pure general relativity. A way to circumvent this problem is to consider an alternative theory of gravity with finite range scalar forces (see for example [32]). In this way the analytical GR solution can be used in the far region from the star where the effective scalar field drops off exponentially and is practically zero. As a particular example of such generalized theories of gravity we consider the f (R) theories and more precisely -R 2 -gravity having a Lagrangian of the form f (R) = R + a R 2 where a is a free parameter. The motivation comes on one hand from the fact that this is one of the most popular and widely used theory possessing an effective finite range scalar force (the R 2 gravity is mathematically equivalent to a particular class of scalar-tensor theories with massive scalar filed). Moreover, the natural application of TLN is connected with the inspiral phase of binary neutron star merger. Currently the only studies of mergers in a theory with effective finite range scalar force are done exactly in R 2 gravity [34] and that is why it is natural to consider the same case.
The f (R) theories, though, are explored mainly on cosmological scales because of the connection to the dark energy problem and compact objects in these theories are more scarcely studied. Non-perturbative models of neutron stars were constructed in the static case in these theories in [35][36][37][38] and they were later extended to the slowly [39] and rapidly rotating [40] cases. Binary neutron star mergers in R 2 gravity were examined in [34]. The results in this paper show that the future gravitational wave observations of merging neutron stars can impose strong constraints on R 2 -gravity.
The goal of the present paper is to calculate the TLN of neutron stars in R 2 -gravity. The current observational constraints on astrophysical scales impose the following upper bound a 10 11 m 2 [41] based on the satellite mission Gravity Probe B, which leaves space for significant deviations from pure GR. As far as merger observations are concerned, only the neutron star mergers can impose serious constraint on f (R) theories because of the presence of no-hair theorems for black holes. Due to the still limited accuracy in these observations and the fact that only one neutron star merger is detected [1], no serious constraint on the parameter a are imposed yet.
Even though we concentrate on f (R) theories, the general framework for calculating TLN developed in the paper is valid for a much larger class of alternative theories of grav-ity, the massive scalar-tensor theories (as we have already commented the TLNs in massless scalar-tensor theories have already been calculated in [31]). The reason comes from the fact that f (R) theories are mathematically equivalent to a particular class of scalar-tensor theory with nonzero scalar field potential. Even more, we use this equivalence explicitly in the present paper in order to simplify the calculations.
The paper is organized as follows. The mathematical framework behind the R 2 -gravity and the way of constructing equilibrium neutron star solutions is examined in Sect. 1. In Sect. 2 the formulas for the calculation of the TLN, both polar and axial, are derived. Section 3 is devoted on the numerical results and the comparison with pure GR. The paper ends with Conclusions.

f(R) theories and equivalence to scalar-tensor theories
The essence of f (R) theories is that the Ricci scalar R in the action is substituted by a function of this scalar f (R): Here R is the Ricci scalar with respect to the spacetime metric g μν and S matter is the action of the matter where the matter fields are denoted by χ . If we want the theory to be well posed, i.e. to be free of tachyonic instabilities and ghosts, the following inequalities should be fulfilled We will work in a particular class of f (R) theories, the so-called R 2 gravity where where a is a parameter and in order to satisfy the above given inequalities one should require that a ≥ 0. A very common approach is to work not with the original form of the action but transform it to another one by substi- . In this way we obtain an action that is equivalent to the Brans-Dicke theory [42,43] with a parameter ω B D = 0: In the case of R 2 gravity the potential takes the form and therefore the scalar field is massive with It is mathematically equivalent to work either with the original form of the action (1) or with its scalar-tensor representation (3) and this was explicitly shown also in [44] for the case of neutron star solutions. Moreover, it is evident that f (R) theories belong to the class of modified gravity with finite range scalar forces, that would be very important later when calculating the TLN of neutron stars.
The action (3) is written in the physical Jordan frame where there is no direct coupling between the matter and the scalar field in order to satisfy the weak equivalence principle. One can further simplify the problem by introducing the Einstein frame by making a conformal transformation of the metric and redefining the scalar field and the potential Thus we arrive at the following Einstein frame action where R * is the Ricci scalar curvature with respect to the Einstein frame metric g * μν . As one can see, a direct couping between the matter and the scalar field appears in the Einstein frame through the coupling function A 2 (ϕ) = −1 (ϕ). In the particular case of R 2 gravity the coupling function and the scalar field potential take the following form The field equations in the Einstein frame are much simpler compared to the Jordan frame ones and that is why we will employ this frame. Of course, the final quantities that we obtain have to be transformed back to the physical Jordan frame. In addition, the equation of state of the nuclear matter that we use will be also only in the Jordan frame. A more detailed discussion of the problem can be found in [36,39,40].
We will consider nonrotating stars and thus the following general ansatz for the static and spherically symmetric Einstein frame metric can be used where all the metric functions depend on r only. The reduced field equations take the following form Here the Jordan frame pressure p and energy density ρ are used and they are connected to the Einstein frame ones ( p * and ρ * ) via the following relations p * = A 4 (ϕ) p and ρ * = A 4 (ϕ)ρ. The Jordan frame quantities ρ and p are naturally connected via the equation of state (EOS) for the neutron star matter p = p(ρ). In addition, we have to impose the standard boundary conditions -regularity at the center of the star and asymptotic flatness at infinity. The radius of the star is calculated from the requirement of vanishing of the pressure at the stellar surface and the mass is taken from the asymptotic expansion of the metric functions at infinity. It is important to note that for the considered R 2 gravity the mass in the Einstein and the Jordan frame coincide, while the physical Jordan frame radius of the star R S is connected to the Einstein frame one r s in the following way Here we have presented the problem of calculating the background equilibrium neutron star solutions very briefly. More detailed explanations can be found in [36].
In what follows, we shall use the dimensionless parameter a → a/R 2 0 , where R 0 is one half of the solar gravitational radius R 0 = 1.47664 km.

Tidal love numbers
In order to compute the tidal Love numbers we have to consider the stationary perturbations of the static and spherically symmetric stars in R 2 -gravity . The perturbations of the metric can be separated in polar and axial type. Here we present the two cases separately and derive the tidal Love numbers in both cases.

Polar
For the polar perturbations the peturbed Einstein frame metric in the Regger-Wheeler gauge can be written in the form where Y lm (θ, φ) are the spherical harmonics. The perturbations of the scalar field, energy density and the pressure can be decompose in the form δϕ = δφ(r )Y lm (θ, φ), δρ * = δρ(r )Y lm (θ, φ) and δp * = δp(r )Y lm (θ, φ). After perturbing the Einstein frame field equations of the f (R) gravity coupled to a perfect fluid it can be shown that H 0 = −H 2 and H 1 = 0. Also K , δρ(r ) and δp(r ) can be written as functiond of H 0 and δφ. Finally we obtain two equations for H 0 = −H 2 = H and δφ governing the stationary perturbations of the static and sphereically symmetric stars in f (R) gravity: Here 0 , ψ 0 , ϕ 0 , p * 0 and ρ * 0 are the corresponding unperturbed variables taken from the background neutron star solutions andc s is the Jordan frame sound speed defined byc 2 s = ∂ p/∂ρ. In the considered model the scalar field mass is nonzero which means that both the background scalar field ϕ 0 and its perturbation δφ drop off exponentially at infinity. This means that the corresponding scalar field tidal Love number is zero.
The asymptotic behavior of H at large r on the other hand is The tidal Love number k 2 is connected to the coefficients in the above given expansion c 1 and c 2 in the following way: In pure GR the ratio c 1 /c 2 is usually determined after matching at the stellar surface the numerical solution for H inside the neutron star to the analytical solution outside it [19][20][21][22]. Applying this approach directly to our problem is not possible since the equations for H and δφ are coupled and in the general case no analytical solution exist outside the star. The fact that the scalar field is massive, though, simplifies the problem considerably. Since both the scalar field and its perturbation die out exponentially at distances larger than the Compton wavelength of the scalar field λ ϕ = 2π/m ϕ , far enough from the surface of the star the scalar field and its perturbation are practically zero and we can use the same analytical solution as in pure GR. This requires, of course, a matching of the inner numerical and outer analytical solutions to be done not at the surface of the star but far enough from the stellar surface where the scalar field and its perturbation are negligible. This large distance where we match the two solutions will be denoted by r match . We have also verified that matching the two solutions not at the stellar surface but at r match , works very well and does not lead to any numerical problems.
The perturbation equation for H far away from the center of the star, where the scalar field and its perturbation are negligible, can be obtained straightforward from Eq. (17) after substituting p * 0 = ρ * 0 = ϕ 0 = δφ = 0. This equation is the same as in pure GR and its analytical solution can be expressed in terms of elementary functions [20]. As commented, the value of k 2 can be calculated after matching the numerical and analytical solutions at large enough radial distances r match . Since r match is connected to the Compton wavelength of the scalar field, r match is not a constant but increases with the increase of the parameter a. For the l = 2 case one obtains where y = r H /H , C 1 = M/R S is the compactness of the star and C = M/r match . The value of y is calculated after solving numerically the coupled system of equations (17), (18) from r = 0 to r = r match . It is important to note, that the polar TLNs are the same in the physical Jordan and the Einstein frame since the scalar field drops off exponentially outside the star.

Axial
In the axial case the metric perturbations are given by where The perturbations of the scalar field, the energy density and the pressure vanish. Using the perturbations of the field equations one can show that h 1 = 0. We are left with only one equation for the metric perturbation h: In this case we also do not have an analytical solution because of the presence of scalar field terms. Similar to the polar case, though, such solution can be found far outside the star where the scalar field has died out exponentially and the solution is the same as in pure GR. The asymptotic equation at such large distances is and it can be solved analytically for a given l [21,22]. The function h has the following asymptotic: The tidal Love numbers for the axial perturbations k axial l are related to the coefficients c 1 and c 2 in the following way The value of k axial l can be found after matching the numerical solution (obtained after integrating the perturbation equation (24) from the center of the star to some large distance r match were the scalar field is negligible) with the analytical asymptotic solution. Thus, one can obtain the following relation for the l = 2 case 2 k axial where y = rh / h, C 1 = M/R S is the compactness of the star and C = M/r match . Similar to the polar case, the axial TLNs are the same in the physical Jordan and the Einstein frame.

Numerical results
We will work with three modern realistic EOS that allow for models with maximum mass above the two solar mass barrier [45,46] and are in agreement with the constraints coming from the observation of the tidal Love numbers of merging neutron stars [1]. These are the APR4 EOS [47], the SLy4 EOS [48] and the MPA1 EOS [49]. We should note, though, that the MPA1 EOS actually do not fit well in the current constrains coming from the electromagnetic observations [50,51] but we included it in our studies so that we can cover a larger range of stiffness. In order to be able to make a better comparison between the three EOS and judge about the effect of R 2 gravity on the background neutron star models, the mass as a function of the radius is plotted in Fig. 1 for these EOSs and for four indicative values of the parameter a. The chosen values of a are the same as the ones used in the calculations of the tidal Love numbers. Since a is back-proportional to the mass of the scalar field (see eq. (5)), a → ∞ corresponds to the massless scalar field case. On the other hand the mass of the scalar field goes to infinity when a → 0 which means that the Compton wavelength is practically zero and the corresponding solutions tend to the pure GR case. As one can see in Fig. 1 the differences of the neutron star masses and radii with pure GR can reach up to roughly 10% and they are comparable both qualitatively and quantitatively with the deviations due to the EOS uncertainties.
Let us comment in more detail the particular values of a that we have chosen. More precisely, we have worked with a ≤ 100 because of the following reason. The parameter a introduces a length-scale related to the Compton wavelength of the scalar field, above which the scalar field drops off exponentially and thus the scalar field has a finite range. 3 This is a crucial ingredient in our calculations of the TLNs since we use the pure GR solution of the perturbation equations outside this Compton radius. Thus, it is natural to require that the Compton wavelength of the scalar field is smaller than the orbital separation between the merging neutron stars at the time when they can be observed by the ground base detectors and the TLNs can be measured. Assuming that with the current instruments we can detect the emitted gravitational waves when the orbital separation drops down to roughly a few hundreds of kilometers, we have chosen to work with a ≤ 100 which leads to λ ϕ ≤ 226 km.
The polar (top panel) and the axial (bottom panel) TLNs as functions of the neutron star compactness are plotted in Fig. 2 for four values of a and only for the APR4 EOS, in order to have better visibility. The equation of state dependence of the results is presented in Fig. 3. As one can see for a fixed Clearly, the changes in the polar Love numbers due to R 2gravity are within the current equation of state uncertainty. Hence, the current generation of gravitational wave detectors is unlikely to be able to set constraints on the parameter a. The three EOSs, though, were chosen to be the ones allowed by the measurement of the TLNs in the neutron star merger [1]. If we take into account also the constraints from the electromagnetic observations [50,51] the picture might change because then the MPA1 EOS is outside the allowed range of masses and radii. Thus, if we consider only the APR4 and Sly4 EOSs, the maximum deviation in R 2 -gravity is larger than the difference between the two equations of state. A definite answer whether this is an observable effect or not can be given only after a detailed analysis of change in the phase of the signal and such a study is underway.
On the other hand the electromagnetic observations are rapidly advancing and the next generation of gravitational wave detectors is already planned. That is why one can expect that when we know the EOS with a better accuracy in the future from the electromagnetic observations, and have more accurate observations of the gravitational waveforms of merging neutron stars, the R 2 -gravity effect will be important, producing effects larger than the equation of state uncertainties. Thus we should take them into account when extracting the relevant parameters from the gravitational wave signal.
The axial TLNs, on the other hand, are significantly influenced (for a ≤ 100) by the modifications of the theory of gravity. Moreover, the absolute value of the axial tidal Love numbers increase compared to the pure GR case. In pure GR the contribution of the axial TLN to the gravitational wave phase is roughly two orders of magnitude smaller that the polar contribution and the contribution of the higher order (higher l) polar TLN would be more important than the axial one as far as the change of the phase in the gravitational wave signal is concerned [22]. The very recent calculations, though, show that the axial TLNs in pure GR would have negligible effect even for the next generation of gravitational wave detectors [52]. Thus, even though the absolute value of the axial TLNs can increase significantly in R 2 gravity, they could be still hardly observed in practice.

Conclusion
In the present paper we have calculated the tidal Love numbers of neutron stars in alternative theories of gravity with finite range scalar field. We have chosen to work with f (R) theories, and more precisely R 2 -gravity, due to the fact that the only studies of binary neutron stars mergers [34] in theories with finite range scalar field are done exactly in R 2gravity. The study is motivated by the fact that the recent detection of gravitational waves from merging neutron stars allowed us to measure their TLN and thus set constraints on the EOS. Since very often there is a degeneracy between effects coming from modifications of GR and uncertainties of the nuclear matter EOS, an estimation of the influence of alternative theories of gravity on the TLN is very important for the proper interpretation of the observational data. As a matter of fact the problem of calculating the TLN in modified gravity is not well studied, with the exceptions of dCS gravity, EiBI gravity and massless scalar-tensor theories.
We have calculated both the axial and the polar TLN for l = 2. The results show that while the polar TLN is only slightly influenced by the modification of GR (the deviations are within the EOS uncertainty is we consider a broader set of EOSs), the axial TLN can be several times larger (in terms of absolute value) compared to pure GR. These conclusions are for values of the free parameter of the R 2 -gravity that fulfill the requirement that the effective radius of action of the R 2 term in the Lagrangian is smaller than the orbital separation between the merging neutron stars when they enter inside the detector sensitivity.
The question is whether such deviations from Einstein's theory can lead to observable effects. The problem is that the quantity that can actually be measured with the current detectors is the polar TLN, while the axial one is expected to give much smaller (yet unmeasurable) contribution in the change of the phase of the signal. The recent calculations suggest that even taking into account that the absolute value of the axial TLN can increase significantly in R 2 -gravity, they could hardly have any measurable influence on the gravitational wave signal even for the next generation of gravitational wave detectors. On the other hand, when the nuclear matter EOS is better constrained by the electromagnetic observations in the future, the deviations in the polar TLN might become observationally important in order to be able to accurately interpreted the detected signal.