Exact versus Taylor-expanded energy density in the study of the neutron star crust-core transition

The importance of the fourth and higher order terms in the Taylor series expansion of the energy of the isospin asymmetric nuclear matter in the study of the neutron star crust-core phase transition is investigated using the finite range simple effective interaction. Analytic expressions for the evaluation of the second and fourth order derivative terms in the Taylor series expansion for any general finite range interaction of Yukawa, exponential or Gaussian form have been obtained. The effect of the nuclear matter incompressibility, symmetry energy and slope parameters on the predictions for the crust-core transition density is examined. The crustal moment of inertia is calculated and the prediction for the radius of the Vela pulsar is analyzed using different equations of state.


Introduction
The study of the equation of state (EOS) of a neutron star (NS) is an area of current research interest in nuclear astrophysics. There is a lack of consensus among the various model calculations about the presence of deconfined quarks and hyperons inside the NS core [1]. However, there is agreement in the fact that the structure of NSs consists of a solid crust encompassing the dense homogeneous core, which is maintained in a charge neutral β−stable condition and is in a liquid phase [2]. The structure of the inner crust is inhomogeneous. It is made of clusters of positive charges in a solid lattice immersed in a bath of neutrons in the superfluid phase, along with a background of electrons such that charge neutrality is maintained. The complicated inhomogeneous structures often referred to as "nuclear pasta" [3,4,5,6,7] makes the formulation of the EOS in this region a formidable task and the value of the transition density between the crust and the core, ρ t , still remains uncertain. Although the NS crust is a small fraction of the star mass and radius, it plays an important role in various observed astrophysical phenomena [3,8,9,10,11,12,13,14,15,16,17,18,19,20].
The microscopic formulations of the inhomogeneous phase [21,22,23,24,25] favour a continuous nature of the transition between the crust and the core. The density at which the crust-core transition takes place can be estimated by examining the onset of violation of the stability condition of the homogeneous liquid core against small-amplitude density fluctuations, which indicates the formation of nuclear clusters. The three more common approaches to study the crust-core transition are the dynamical method [3,8,9,26,27,28,29,30], the thermodynamical method [13,31,32,33,34,30,35] and the random phase approximation [36,37]. The dynamical and thermodynamical methods are similar in the long wavelength limit giving similar results, but with a slightly smaller value of the transition density in the former case due to the inclusion of density gradient and Coulomb effects [30]. The inequalities resulting from the stability conditions show a direct connection with the isospin part of the EOS. The isospin dependent part of the energy per particle in isospin asymmetric nuclear matter (ANM) is considered as the symmetry energy. The isospin dependence appears in a complicated way in the analytic formulation of the EOS of ANM in any of the model calculations. In order to have analytical simplifications, often the Taylor series expansion of the energy density around zero isospin asymmetry is adopted and the energy density in ANM is expressed as a combination of even powers of the isospin asymmetry parameter. In the many-body model calculations using the Taylor series expansion of the energy density, it has been found that higher order terms beyond 2 nd order have a very small contribution [38]. This finding leads to the popular parabolic approximation of the energy per particle in ANM. Under this parabolic law, the symmetry energy can be expressed either from the 2 nd order term in the Taylor series expansion of the energy density or by the empirical relation defined as the difference between the energy per particle in pure neutron matter (PNM) and in symmetric nuclear matter (SNM).
Although the parabolic approximation is sufficient for the energy calculations, the 4 th and higher order terms in the isospin asymmetry can have some influence when the convergence of the Taylor series is rather slow, as for example in the calculation of the proton fraction in β−equilibrated n + p + e + µ matter. Thus the contributions of these higher order terms may be significant in deciding as to whether the direct URCA process occurs or not in NSs [39]. The prediction of the crust-core transition density ρ t in NSs is another important area where the role of the higher order terms is quite significant and the results with the parabolic approximation can be misleading. There are a limited number of studies [30,35,40,41], using either non-relativistic or relativistic mean field model calculations, aimed at explicitly demonstrating the important influence of 4 th and higher order terms in the prediction of the critical transition density. A general trend that emerges from these calculations is that the inclusion of the 4 th order term reduces the value of the crust-core transition density. This feature points out to the fact that the complete isospin contribution is crucial in the calculation of the transition density.
In non-relativistic model calculations of the crust-core transition using Skyrme type interactions, the study of the influence of the 4 th and higher order terms in the Taylor series expansion of the energy can be easily performed due to the simple form of the analytic expressions. However, for finite range interactions the complexity increases. In Ref. [30], the influence of higher order terms of the Taylor series expansion on the crust-core transition and the properties of a NS has been studied using the MDI energy density [42] that can result from a finite range interaction having a Yukawa form. In the present work, we have obtained general expressions that can be used to obtain the 2 nd and 4 th order contributions of the Taylor series expansion for any conventional form of finite range interactions. The values of ρ t and the pressure at the transition density, P (ρ t ), are the quantities of crucial importance in the calculation of the crustal fraction of the moment of inertia ∆I/I, where I is the moment of inertia of the NS [10,30]. The quantity ∆I/I can be estimated from the observations of pulsar glitches. Glitches are the intermittent sudden jumps in the rotational frequencies of pulsars, which are believed to arise from the interaction between the solid crust and the liquid core.
In section 2, we outline the stability conditions of the thermodynamical method in terms of the nuclear EOS. Analytic expressions for the 2 nd and 4 th order terms in the Taylor series expansion of the energy density in ANM have been obtained. The explicit expressions for the energy per particle in ANM, neutron (proton) chemical potential and other quantities used in the study have been obtained for the Yukawa form of the finite range simple effective interaction (SEI) and are given in the Appendix. Alternatively, and for practical purposes, these expressions can be evaluated numerically (including thermodynamic derivatives). The procedure adopted to determine the parameters of SEI necessary for the study of ANM is outlined. The approximate expression of the crustal momentum of inertia obtained in Ref. [10,30] is outlined in this section. Section 3 contains results and discussions on the transition density obtained with the Yukawa form of the SEI for the 2 nd and 4 th order Taylor series, as well as for the exact analytic case. We discuss the importance of using exact expressions of the energy density instead of their Taylor expansions to compute the crust-core transition density. The influence of nuclear matter saturation properties such as incompressibility, symmetry energy and slope parameter on the transition density are also discussed. The crustal moment of inertia as well as the mass-radius relation for the Vela pulsar are calculated for the EOSs under consideration. The influence of the functional form of the finite range part of the interaction on the transition density is also examined by considering a Gaussian form factor. Finally, Section 4 contains the final remarks and conclusions.

Formalism
The study of the crust-core phase transition in neutron stars in the thermodynamical model amounts to examine the stability of the β−equilibrated charge neutral dense matter in the liquid phase with respect to the stability conditions [31,32], In these equations, P is the total pressure, v = V /B is the volume per baryon number, q is the charge fraction given by q = Y p − Y e , where Y i = ρ i /ρ, i = p, e, are the proton and electron fractions, with ρ = ρ p + ρ n being the total nucleonic density and ρ n , ρ p and ρ e being the neutron, proton and electron densities, respectively. The βequilibrium chemical potential is µ = µ n − µ p = µ e , µ i , where i = n, p, e denotes the respective chemical potentials. The stability conditions in equations (1) and (2) refer to the mechanical and chemical stabilities of the system, respectively. Here, the total pressure P of the (n,p,e) system comprises the baryonic contribution P N and the leptonic contribution P l . Under the consideration of the relativistic Fermi gas model, the pressure of the electrons is P l = µ 4 e 12π 2 (c ) 3 . The electronic contribution to the pressure becomes a function of µ only when the system is in β-equilibrium and does not contribute to the stability condition in equation (1), which now reads Now using the relations v = 1 ρ and µ = − ∂e(ρ,Yp) ∂Yp , where e(ρ, Y P ) is the energy per nucleon in the β-equilibrated nuclear matter (µ = µ e ), the conditions in equations (3) and (4) identically reduce to [34,30,43] In the second condition one may recognize the last term in the right hand side (rhs) as the screening length. The first term in the rhs can be shown, under the parabolic approximation, to be proportional to nuclear symmetry energy E s (ρ). Although the density behaviour of E s (ρ) is still uncertain, the heavy-ion collision studies favour a monotonically increasing behaviour of E s (ρ) and therefore the condition in equation (6) is usually satisfied. The first two terms in the rhs of equation (5) refer to the pressure and incompressibility of nucleonic matter and are positive for fundamental reason, whereas the third term which comes from the presence of leptons contributes negatively. Thus the stability condition in equation (5) reduces to the form where V thermal = ∂P N ∂ρ µ and e(ρ, Y p ) = H(ρ,Yp) ρ is the energy per particle in ANM, H(ρ, Y p ) being the energy density of ANM. From the relation in equation (7) the importance of the EOS of ANM in determining the transition density for the crust-core transition in a NS is evident. Once the equilibrium proton fraction Y p is ascertained for a given nucleonic density ρ from the solution of the charge neutral beta stability conditions, the nucleonic energy per particle e(ρ, Y p ) can be calculated and V thermal evaluated using equation (7) either numerically or by analytical procedure wherever it is possible. Numerical evaluation of equation (7), however, does not provide further insight into the isospin dependence of the EOS of ANM which is an area that is less understood and of current research interest. The analytical calculation of V thermal from equation (7) can be conveniently made provided the Y p -dependence in the expression for the energy density in ANM is separated out. The isospin dependence in the energy density of ANM resulting from a finite range interaction is complicated and in most cases it is not practically possible to express the energy density by segregating the Y p -dependent parts. Due to this, the Taylor series expansion of the energy density of ANM is usually adopted in the studies involving the isospin dependence aspects. For calculating the exact analytic expression of V thermal , equation (7) can be equivalently expressed in terms of derivatives of the chemical potentials µ q (q = n, p) with respect to the neutron and proton densities, ρ n and ρ p , as given by where µ i = ∂H(ρ,Yp) ∂ρ i for i = n, p. Further, the equality ∂µn ∂ρp = ∂µp ∂ρn is used. The exact isospin dependence of the EOS of ANM is poorly known and a popular way of studying it is by making a Taylor series expansion of the energy density in ANM in even powers of the isospin asymmetry parameter, β = ρn−ρp with, H sym,4 (ρ) = 1 384 where H(ρ) = ρ e(ρ) is the energy density in SNM. The analytic expressions of these 2 nd and 4 th order terms for any general effective interaction, v(r), are given as, and respectively. In these expressions k f = (3π 2 ρ/2) 1/3 is the Fermi momentum in SNM, v l d (r), v ul d (r) and v l ex (r), v ul ex (r) are the direct (d) and exchange (ex) parts of the interaction acting between like (l) and unlike (ul) pairs of nucleons, m is the nucleon mass and j 0 and j 1 are the spherical Bessel functions of zeroth and first order. In all the model calculations it is found that the contributions to the energy from the higher order terms, beyond the 2 nd order, is very small and therefore the series is often truncated at the quadratic term resulting into the parabolic expression of the energy per particle in ANM, where the symmetry energy E s (ρ) is obtained from the second order derivative term as H sym,2 (ρ)/ρ. Another empirical expression of the symmetry energy E s (ρ), often used, is by approximating it as the difference between the energy per particle in PNM, e N (ρ), and in SNM, e(ρ), The justification of this approximation is that at the two extreme values of the isospin asymmetry, Y p =0 and Y p =1/2, the energy per particle expression e(ρ, Y p ) = e(ρ) + (1 − 2Y p ) 2 E s (ρ) of ANM reduces to the respective expressions of PNM and SNM, respectively. For the 2 nd order approximation of the energy density in the Taylor series expansion, the stability condition V thermal in equation (7) becomes [30,40,43], where E sym,2 = H sym,2 ρ . Similarly, for the 4 th order approximation, where the Taylor series expansion in equation (9) is truncated at the fourth order term, the stability condition becomes [40,41], (17) where E sym,4 = H sym,4 ρ . The difference between the neutron and proton chemical potentials in the β− stability condition, for the 2 nd and 4 th order Taylor series approximations of the energy density becomes, and respectively.
The stability condition of V thermal is a signature of the crust-core transition in the β−equilibrated neutron star matter. As we summarize in the next subsection, the transition density and the pressure at the transition density play a critical role in the prediction of the crustal fraction of the moment of inertia of neutron stars, used in the possible explanation of the observed glitches in pulsars.

Crustal fraction of the moment of inertia of neutron stars
Based on the hypothesis that the mechanism for glitches observed in the magnetized rotating neutron stars, is due to the pinning of the vortexes in the superfluid neutrons inside the dense liquid core with the superfluid neutrons of the inner crust in the crustcore transition region [44,45,46,47], the crustal fraction of moment of inertia, ∆I/I, can be calculated from the observed glitches. Glitches are the intermittent disruption in the extremely regular pulses emitted from the magnetized rotating neutron stars. An approximate expression for ∆I/I has been obtained by Xu et. al. [30] using the work of Link et. al. [10]. It contains the mass M and radius R of the NS, and the dependence on the EOS through the pressure and the transition density at the crust-core transition, P (ρ t ) and ρ t , respectively, as given by where ξ = GM Rc 2 , G is the gravitational constant and c is the velocity of light. The NS mass and radius for a given EOS are calculated by solving the Tolman-Oppenheimer-Volkov (TOV) equations. The total moment of inertia, I, of the neutron star rotating slowly with a uniform angular velocity Ω is obtained from I = J/Ω, where the total angular momentum J is calculated from [48], The angular velocity of a point in the star,ω, is obtained from the solution of the relevant equation in Ref. [48] subject to the boundary conditions thatω is regular at r → 0 andω → Ω as r → ∞.

Finite range simple effective interaction
The present study is made using the finite range simple effective interaction (SEI) described in [49,50]: where r = r 1 − r 2 , R = ( r 1 + r 2 )/2 and f (r) is the functional form of the finite range interaction, which may be of Yukawa, Gaussian or exponential form, and contains a single range parameter α. The SEI contains altogether eleven parameters, namely, t 0 , x 0 , t 3 , x 3 , b, γ, α, W , B, H and M. However, for the study of isospin asymmetric nuclear matter only nine parameters are required, namely, b, γ, α, ε l 0 , ε ul 0 , ε l ex , ε ul ex , ε l γ , ε ul γ . The connection of the new parameters to the interaction parameters is given in Ref. [51]. In this work we shall be calculating the results for the Yukawa form of f (r). The determination of the nine parameters of the ANM is discussed in detail in earlier studies [51,52,53]. Here we outline the procedure for the sake of convenience of the reader.
The symmetric nuclear matter requires only the following combinations of the strength parameters in the like and unlike channels: together with the γ, b and α parameters, altogether six parameters. For a given value of the exponent γ, and assuming the standard values for the nucleon mass mc 2 , saturation density ρ 0 and binding energy per particle in SNM at saturation e(ρ 0 ), the remaining five parameters ε 0 , ε γ , ε ex , b and α of symmetric nuclear matter are determined in the following way. The range α and the exchange strength ε ex are determined simultaneously by adopting an optimization procedure using the condition that the nuclear mean field in symmetric nuclear matter at saturation density vanishes for a nucleon kinetic energy of 300 MeV, a result extracted from optical model analysis of nucleon-nucleus data [54,55,56]. This requires only the values of mc 2 , ρ 0 and e(ρ 0 ) as discussed in Refs [49,50] and is independent of the other parameters of symmetric nuclear matter including γ. The parameter b is fixed independently for avoiding the supraluminous behaviour in SNM [57]. The two remaining parameters, namely ε γ and ε 0 , are obtained from the saturation conditions. The stiffness parameter γ is kept as a free parameter and its allowed values are chosen subject to the condition that the pressure-density relation in symmetric matter lies within the region extracted from the analysis of flow data in heavy-ion collision experiments at intermediate and high energies [58]. It is found that the maximum value that fulfills this condition is γ=1, which corresponds to a nuclear matter incompressibility K(ρ 0 )=283 MeV. Therefore, we can study the nuclear matter properties by assuming different values of γ up to the limit γ=1. Now, to describe asymmetric nuclear matter we need to know how the strength parameters ε ex , ε γ and ε 0 split into the like and unlike components. The splitting of ε ex into ε l ex and ε ul ex is decided from the condition that the entropy density in pure neutron matter should not exceed that of the symmetric nuclear matter. This prescribes the critical value for the splitting of the exchange strength parameter to be ε l ex = 2ε ex /3 [59]. The splitting of the remaining two strength parameters ε γ and ε 0 , is obtained from the values of the symmetry energy E s (ρ 0 ) and its derivative E dρ 0 at saturation density ρ 0 . Notice that the usual slope parameter of the symmetry energy is By assuming a value for E s (ρ 0 ) within its accepted range [60,61], we determine E ′ s (ρ 0 ) from the condition that the difference between the energy densities of the nucleonic part in charge neutral beta-stable n + p + e + µ matter, referred to as neutron star matter (NSM), and in symmetric matter at the same density be maximal [52]. The value of E ′ s (ρ 0 ) thus obtained for a given E s (ρ 0 ) predicts a density dependence of the symmetry energy which is neither very stiff nor soft and does not predict the direct URCA process in the calculated NSs. It may noted here that the splitting of the exchange strength parameter ε ex into the like and unlike channels in ANM is solely determined, as discussed in Ref [59], from the thermal evolution of nuclear matter properties in SNM and PNM, and is independent of the splitting of the parameters ε 0 and ε γ . Thus, the parameters associated with the finite range exchange part of the mean field and EOS, namely ε ex , α, ε l ex and ε ul ex are independent of the way in which the remaining parameters of the interaction, including γ, are determined and also of the choice of E s (ρ 0 ).
With the parameters determined in this way, the SEI was able to reproduce the trends of the EOS and the momentum dependence of the mean field properties with similar quality as predicted by microscopic calculations [59,62]. The two sets of parameters of the EOSs of SEI corresponding to γ=1/2 and γ=1, which cover a range of nuclear matter incompressibility between 240 MeV and 280 MeV, are given in Table 1, along with their nuclear matter saturation properties. The relevant analytic expressions of the various quantities required in the study for the Yukawa form of SEI are given in the Appendix. As it can seen from Table 1, the SEI values of the saturation density and binding energy per nucleon in SNM lie within the empirical ranges 0.17±0.03 fm −3 and 16±0.2 MeV [51]. The value of the symmetry energy E s (ρ 0 ) also agrees well with the range 29-33 MeV provided by recent analyses [63,64]. The slope parameter of the symmetry energy L(ρ 0 ) is also compatible with the commonly accepted range between 40 and 70 MeV [65]. The analysis of the excitation energy of the giant monopole resonance provides a range of allowed values of the incompressibility K between 200 and 260 MeV [60]. The SEI used in this work have K values around 240 MeV (γ=1/2) and 280 MeV (γ=1). Although this second value is slightly high, we use it to simulate a stiffer EOS. However, we can see in the middle panel of Figure 1 that the two considered EOS computed with the SEI used in this work lie within the boundaries of allowed values extracted from the analysis of the flow data in heavy-ion collisions [58] and kaon production data [66]. In the left panel of Figure 1 we compare the energy per particle in SNM and PNM, e(ρ) and e N (ρ), respectively, as a function of the density computed using the SEI corresponding to γ=1/2 with the microscopic Dirac-Brueckner-Hartree-Fock calculation using the Bonn B potential [67] and with the variational calculation with the A18+δv+UIX * realistic interaction [68]. It can be seen that, for both SNM and PNM, there is a good agreement between the SEI predictions and the microscopic results up to about ρ = 0.3f m −3 . Above this density, the SNM and PNM EOSs obtained with SEI follow rather well the trend of the variational calculation and differ more from the DBHF results, which grows faster in SNM and PNM compared with the corresponding Table 1. Values of the nine parameters of ANM for the two sets of SEI corresponding to γ=1/2 and γ=1 together with their nuclear matter saturation properties (see text for details). Nuclear matter properties at saturation density where, u n(p) is the neutron (proton) mean field. This quantity can be taken as a measure of the momentum dependence of the neutron (proton) single particle potentials. The neutron-proton effective mass splitting at saturation density and momentum equal to the Fermi momentum at saturation density for the SEI having γ=1/2 is displayed as a function of the asymmetry β in the rightmost panel of Figure 1 compared with the microscopic DHFB [67] and BHF [69] predictions. It can be seen that the SEI results compare very well with the DBHF results over the whole range of asymmetry, while they differ more from the BHF prediction [69]. This fact shows that, at least in normal nuclear matter, the momentum dependence of the SEI mean fields are similar to those of the ab initio DBHF formulation.

Results and Discussion
In order to examine the convergence of the Taylor series expansion of the energy density, we calculate the 2 nd and 4 th order contributions, given in equations (12) and (13), respectively, for our SEI in equation (23) as a function of density. The results are shown in the Figure 2 for the EOS of stiffness γ = 1/2. The 4 th order contribution can be seen to be very small, at most represents 2% of the 2 nd order one over the whole range of density, justifying the validity of the Taylor series expansion of the energy density. It may be mentioned here that only the kinetic and finite range exchange parts of the interactions contribute to the 4 th order term, and its functional dependence in the case of SEI makes it negative beyond a density ρ ≈1.55 f m −3 . For the sake of comparison we have also plotted the symmetry energy density obtained under the empirical parabolic approximation (PA). The results of the 2 nd order Taylor series expansion and the PA compare well differing at most within 2% over the whole range of density. The equilibrium proton fraction has been calculated for the two EOSs of table 1 in the exact analytic case and in the 2 nd and 4 th order approximations of the Taylor series expansion of the energy density by solving the corresponding charge neutrality and β−equilibrium conditions. The β−equilibrium conditions for the 2 nd and 4 th order approximations are given in equations (19) and (20), respectively. In the PA, the proton fraction is calculated from equation (19) with E sym,2 replaced by equation (15). The proton fraction for the EOS of γ =1/2 is shown as a function of density in the left panel of Figure 3. In the right panel, the asymmetric contribution to the nucleonic part of the energy density S N SM (ρ) = [H(ρ, Y p ) − H(ρ)] in the β−equilibrated n + p + e + µ matter, for the four cases mentioned before, is shown as a function of density. The curves of the proton fraction obtained under the different approximations lie very close to each other over the whole density range, with the 2 nd order result at the bottom. The small differences between the 2 nd order Taylor series and the exact curves at different densities, is the measure of the cumulative contributions of all higher order terms of the Taylor series to the equilibrium proton fraction in β−stable matter. The curve for the empirical parabolic approximation of equation (15), where the symmetry energy is defined as the difference between the energy per particle in PNM and SNM, however, remains above the three curves in the density range beyond three times the normal nuclear matter density. The small differences in the composition obtained under the various approximations to the exact energy density do not have any noticeable influence on the nucleonic part of the energy of the β−stable matter. This is evident from the right panel of Figure 3, where the curves of the different approximations practically overlap with the curve of the exact result over the whole range of density. The observations are similar for the EOS of table 1 corresponding to γ =1. The mass-radius relationship in neutron stars is obtained by solving the TOV equations associated to the two EOS given in Table 1. For each EOS, the calculations are performed exactly and using the three approximations discussed in the text. The corresponding results are displayed in Figure 4. It can be seen that the 2 nd and 4 th order Taylor series approximations nicely reproduce the exact result in the whole range of considered densities. This is also the case of the parabolic approach, although some small differences respect to the exact calculation appear for masses of the neutron stars below 1.5 M ⊙ . The reason of the success of the different approaches used in this work to reproduce the exact mass-radius relationship lies on the fact that using these approximations the total (nucleonic+leptonic) energy density and pressure reproduce and e N (ρ), shown as a function of density ρ for the EOS of SEI having γ = 1/2. The corresponding results of microscopic DBHF [67] and variational calculations using a realistic interaction [68] have been given for comparison. (middle panel) Pressuredensity relation in SNM of the two EOSs of SEI having γ=1/2 and 1 shown along with DBHF predictions [67] and the allowed regions extracted from the heavy-ion collision studies (area within the red boundary) [58] and the kaon production data (area within the green boundary) [66]. (right panel) The neutron and proton effective mass difference in normal nuclear matter is shown as a function of asymmetry β for the SEI and compared with the predictions of the microscopic DBHF [67], BHF+3BF and EBHF+3BF [69] calculations. very accurately the exact values, as it can be seen in the two panels of Figure 5. The largest differences at ρ = 2 fm −3 , where the exact energy density and pressure are about 3209 and 1781 MeVfm −3 , respectively, are less than 2 MeVfm −3 for the energy density and 4 MeVfm −3 for the pressure. In order to obtain the mass-radius relationship, our EOS, defined in the core, has been supplemented from a density 0.0582 fm −3 down by the EOS in the crust provided by the Baym-Bethe-Pethick calculation [3,8]. However, it should be pointed out that to join the EOS in the core computed with a given model with the EOS of the crust calculated with a different model, may produce a sizeable effect on the radius of the lightest neutron stars [70], although this aspect does not change the conclusion of our present study. The maximum masses and radii of a NS obtained in the exact as well as in the approximated calculations are 1.88 M ⊙ and 10.12   We shall be using these two EOSs for the crust-core transition study in order to examine the influence of the nuclear matter incompressibility on the transition density. The neutron star crust-core transition is another area where the contributions of higher order terms in the Taylor series expansion cannot be ignored, since the Taylor series approximation of the energy density of ANM can be misleading in the prediction of important properties of NSs. One of such properties is the crustal fraction of the moment as a function of the radius R in km obtained with the exact analytic expression of the energy density, the 4 th and 2 nd order Taylor expansion and empirical parabolic approximation of the energy density, for the two sets of EOSs corresponding to γ =1/2 and γ =1 given in Table 1. of inertia, which critically depends on the pressure at the transition density, as can be seen from equation (21).
The density at which the crust-core transition takes place is calculated from the respective stability conditions in equations (16) and (17) for the 2 nd and 4 th order Taylor series approximations. The corresponding result for the exact treatment of the energy density is also calculated using the stability condition expression in equation (8). The pressure at the transition density P (ρ t ) is a crucial quantity for the calculation of the crustal fraction of the moment of inertia used in the possible explanation of the observed glitches in pulsars (see equation (21)). The results for the transition density ρ t , pressure at the transition density P (ρ t ) and corresponding proton fraction Y P (ρ t ) obtained from the three calculations in each of the two EOSs of table 1 are listed in   [40,30]. Though the results of the different parameter sets differ in magnitude, the same trend of predicting a lower transition density (and lower pressure at the transition density) upon inclusion of the 4 th order term in comparison to the 2 nd order result is found in all the sets of both the relativistic and non-relativistic calculations. Moreover, although the inclusion of the 4 th order term moves the transition density and pressure at the transition density in the right direction, the values are still far from the results obtained in the exact calculation. This clearly shows that in order to have the correct prediction of the transition density of Table 2. Transition density ρ t , pressure at transition density P (ρ t ) and equilibrium proton fraction at transition density Y p (ρ t ) for the exact analytic case, 4 th and 2 nd order Taylor series approximations and the empirical PA of the two EOSs corresponding to γ =1/2 and γ =1 of SEI given in table 1. The results of the calculations with two RMF lagrangians [40] and two non-relativistic interactions (MDI and Skyrme R σ ) [30] are given for comparison. The transition density, ρ t , increases with an increase in the stiffness of nuclear matter.
In SEI, an increase from 0.0788 fm −3 to 0.0845 fm −3 is obtained for the exact case as K(ρ 0 ) increases from 237 MeV (γ = 1/2) to 282 MeV (γ = 1). Similarly, the increase in the values of the transition density in the case of the 4 th order approximation is from 0.0919 fm −3 to 0.0962 fm −3 , and from 0.0954 fm −3 to 0.0994 fm −3 in the 2 nd order case, as γ increases from 1/2 to 1. This is illustrated in Figure 6, where V thermal is plotted as a function of density in the close vicinity of the transition density for these three approximations, in the two EOSs having γ =1/2 and γ =1. The results of ρ t for the exact and the 2 nd and 4 th order Taylor expansion of the energy density given in Table 2 suggest the possibility that the transition density obtained using the Taylor series expansion might not have a convergent behaviour. In order to verify this possibility we examine the behaviour of various properties, namely, energy per particle e(ρ, Y p ), pressure P N (ρ, Y p ) and V thermal computed exactly and using the Taylor expansions and the empirical PA of the energy density. For our analysis we have chosen two densities, ρ = 0.08 fm −3 and 1.2 fm −3 , the former lies in the region of the crust-core transition while the latter corresponds to the NS central density region. The behaviour of e(ρ, Y p ) and P N (ρ, Y p ) as a function of Y p is shown in Figures 7 and  8, respectively. From these figures it can be seen that the exact results and the ones obtained with the different approximations considered in this work are quite similar over a wide range of asymmetry in both the low and high density regions. The small discrepancy between the results obtained exactly and with the Taylor expansion at very large asymmetry, particularly in P N (ρ, Y p ), may require some additional higher order terms in the Taylor expansion series for having a perfect convergence with the exact results. However, a very different behaviour is found if the same study is performed with V thermal , in particular in the crust-core transition density domain. V thermal as a function of Y p is computed using equations(8), (16) and (17) for the exact and 2 nd and 4 th order Taylor expanded cases, respectively. To obtain V thermal in the PA, we use equation (16) together with the definition for symmetry energy in equation (15). The corresponding results, computed at our reference densities ρ = 0.08 fm −3 and ρ = 1.2 fm −3 , are displayed in the two panels of Figure 9. From the right panel of figure 9 it can be seen that in the high density region the agreement between V thermal computed exactly and using the Taylor expansions and the PA of the energy density is similar to that found with e(ρ, Y p ) and P N (ρ, Y p ). The agreement between 4 th order and exact results is very good for proton fractions Y p larger than 0.05. The small differences for low proton fractions with Y p smaller than 0.05 could be, in principle, accounted for by considering additional higher order terms in Taylor series. However, the situation dramatically changes in the low density regime with ρ = 0.08 fm −3 , as it can be observed in the left panel of figure 9. For proton fractions Y p larger than 0.2 the agreement between the V thermal computed exactly and using the 4 th order Taylor expansion is good, but the exact and 4 th order predictions of V thermal exhibit a completely different behaviour when the proton fraction decreases below 0.2. For these low values of the proton fraction the exact V thermal shows a stiff rising behaviour for small values of Y p . The curves of V thermal computed using the 2 nd and 4 th order approximations also show an increasing trend when the proton fraction Y p decreases, however their slopes are clearly smaller than in the exact calculation, which magnifies the differences between the exact and approximated calculations of V thermal for small proton fractions. Although the curve of V thermal for the 4 th order approximation shows some improvements over its 2 nd order counterpart, it cannot reproduce the stiff increasing shown by the exact V thermal curve. The results for the V thermal obtained under the empirical PA coincide with the curves for the 2 nd order Taylor expansion of the energy density. Since the equilibrium proton fractions Y p (ρ t ) at the transition density for the different sets of SEI given in table 2 lie within the range 0.02-0.03, where the behaviour of V thermal for the exact case is sharply increasing, it is unlikely that V thermal obtained by including terms of higher order of the Taylor series expansion of the energy density could reproduce the results for ρ t of the exact calculation in the case of our SEI. In order to understand the stiff rising behaviour of V thermal at high asymmetry in the low density region, the contributions of various parts of the EOS, namely, the kinetic, the finite range exchange (ε l ex + ε ul ex ) part and the zero range (ε 0 + ε γ ) part to V thermal have been calculated at the two reference densities, ρ = 0.08 fm −3 and 1.2 fm −3 , as a function of Y p . This is shown in the two panels of Figure 10, where the Y p -dependence of the contributions of the kinetic and finite range exchange (kin + exch) and of the ε 0 -and ε γ -parts (ε 0 +ε γ ) for the exact case and the 4 th order Taylor approximation calculated from equations (8) and (17), respectively, is displayed at these two reference densities. The contributions of the zero range (ε 0 +ε γ )part are the same for both the exact case and the 4 th order Taylor approximation. But the (kin + exch) contributions, which are negative in both cases, are different as Y p decreases, giving a sharp rising behaviour in the exact case. This behaviour of the (kin + exch)-part of the exact case at low values of Y p is a characteristic feature at all the densities. But the (ε 0 +ε γ )-part that gives a high positive contribution at higher density, as can be seen from the right panel of figure 10, dominates and overshadows this rising feature of the (kin + exch)-part at high density. Therefore, the behaviour of V thermal for the exact case and 4 th order Taylor expansion is nearly similar at high density as has been obtained in the right panel of figure 10. However, the (kin + exch) and the (ε 0 +ε γ ) contributions are opposite and comparable in magnitude at low values of the density as can be seen from the left panel of figure 10, thereby, the characteristic stiff rising behaviour of the (kin + exch)-part at low Y p values is manifested in V thermal for the exact case shown in the left panel of figure 9.
The influence of the symmetry energy parameter E s (ρ 0 ) on ρ t is examined by calculating the transition density for two additional values E s (ρ 0 ) = 30 MeV and 36 MeV in the EOS for γ =1/2, besides the value E s (ρ 0 ) = 33 MeV given in table 1. The splitting of the two SNM parameters ε 0 and ε γ (recall equation (24)) in ANM for each value of E s (ρ 0 ), is made by searching for the value of E γ =1), where E ′ s (ρ 0 ) is determined from the constraint outlined in subsection 2.2, are given in table 1. The procedure of determination of the parameters in ANM adopted for SEI does not allow to make an independent study of the influence of E s (ρ 0 ) and E ′ s (ρ 0 ) (i.e. L(ρ 0 )) on the transition density. The influence of L(ρ 0 ) on the transition density was examined in Refs. [30,35] for the MDI interaction (see figure 6 of Ref. [30] and figure  4a and b of Ref. [35]). In order to compare with the trend of MDI, in the present work we have varied, in each of the two EOSs of table 1, the stiffness of the density dependence of the symmetry energy by starting from a low value of E ′ s (ρ 0 ) and assigning increasing values to E ′ s (ρ 0 ) (now, without imposing the constraint used in subsection 2.2 to fix E ′ s (ρ 0 ) uniquely for a given value of E s (ρ 0 )). The change of ρ t with the variation of E ′ s (ρ 0 ) is depicted in Figure 11. The trend obtained with SEI for the 2 nd and 4 th order Taylor series expansion and the exact case is the same overall as obtained with MDI in Ref. [30]. It points out to the fact that the Taylor series approximated results cannot reproduce the trend of the exact calculation. With the exact analytic expression of the energy density, ρ t shows a decreasing trend with an increase in the value of the slope parameter of the symmetry energy, whereas in the Taylor series approximations ρ t attains a minimum value and thereafter follows a slow increasing trend with increasing slope parameter. In order to examine the influence of the functional form of the finite range part of the interaction on the transition density, calculations have been done with the Gaussian form of f (r) in equation (23)). Here we have constructed the EOS with the Gaussian form of f (r) equivalent to the γ=1/2 EOS of the Yukawa form given in table 1. For the same ρ 0 and E s (ρ 0 ) as of the Yukawa γ=1/2 EOS, the γ value required in the Gaussian form of SEI to predict the same nuclear matter incompressibility is γ=0.42. Following the same parameter determination procedure, the parameters of this equivalent EOS were determined and the crust-core transition density was calculated using the Gaussian EOS for the exact case and for the Taylor series 4 th and 2 nd order cases. The transition densities obtained in these cases are ρ t =0.0788 f m −3 , 0.0916 f m −3 and 0.0952 f m −3 , respectively. These results are very similar to the ones obtained with the Yukawa form of the EOS with γ=1/2. The pressure and proton fraction at the transition density also are found to take similar values as in the Yukawa form of the EOS with γ=1/2. This shows that the use of Yukawa or Gaussian finite range form factors, at least in the SEI case, bears little influence on the prediction of the crust-core transition density in NSs.
The mass and radius of the neutron star as functions of the central density ρ c of the star are calculated under the Taylor series approximations and the exact analytic case in each of the two EOSs of table 1 by solving the corresponding TOV equations. The total moment of inertia of the star is also calculated as a function of the central density from I = J/Ω for the considered cases. Using the mass, radius, transition density and pressure at the transition density of the NS, the crustal fraction of the moment of inertia, ∆I I , is calculated from equation (21). Now using the calculated I, the contribution of the crust can be obtained. The results for the mass, radius, crustal fraction of the moment of inertia and crustal thickness δR (distance from the point of the transition density to the surface of the neutron star) are given in Table 3 at different values of the NS central density ρ c for the EOS having γ =1/2. Comparison of the results shows that although the use of the Taylor expansion of the energy density does not have any influence on the prediction of NS bulk properties such as mass and radius, its influence prominently manifests in the predictions of the crustal fraction of the moment of inertia and the crustal thickness through the dependence on ρ t and P (ρ t ). The two important quantities ∆I I and δR are clearly overestimated if one uses the Taylor expansion of the energy density in the calculation (see table 3). The comparison of the values of these two quantities computed exactly and with the 2 nd and 4 th order Taylor approximations and reported in table 3, shows a similar trend to the one exhibited by the transition density.
We have further examined the consequences of the Taylor series expansion of the energy density on the prediction of the mass-radius constraint of the Vela pulsar. From the analysis of the observed data on glitches of the Vela pulsar, the lower limit for ∆I I is obtained to be 0.014 [10]. Using this condition in equation (21) the radius of the Vela pulsar is expressed as a function of its mass, in the form of a straight line equation for each considered EOS and approximation used. The results are given in Table 4. The slope of the straight line relations for the 2 nd and 4 th order Taylor series approximations decreases as compared to the exact calculation, predicting relatively smaller radius for the NS. This is shown in the left panel of Figure 12, where the straight line curves depicting the predictions of the radius of the Vela pulsar as a function of the mass are plotted for the EOS having γ =1/2 and E s (ρ 0 ) =33 MeV of table 1, under the different considered approximations. In the right panel of figure 12, similar results for the EOSs having different values of γ and E s (ρ 0 ) are given in the exact treatment of the energy density. It can be seen from figure 12 that the predictions for the radius in the Vela pulsar for the EOSs having γ =1/2 and γ =1 with the same E s (ρ 0 ) are very close to each other. But the results show a marked sensitiveness on the value of the symmetry energy E s (ρ 0 ). With increase (decrease) in E s (ρ 0 ) the predictions for the radius in the Vela pulsar decrease (increase), as can be seen in the right panel of figure 12 from the results of the three EOSs having the same stiffness γ =1/2 but different symmetry energy parameter E s (ρ 0 ) =36, 33 and 30 MeV.

Summary and conclusions
The influence of the higher order terms of the Taylor series expansion of the energy density in the study of the crust-core transition in neutron stars is investigated. The calculations are performed with the 2 nd and 4 th order Taylor expansions and with the exact treatment of the energy density for different equations of state based on finite range nuclear forces. We find that the parabolic approximation of the energy in ANM, often used in the nuclear calculations, may be misleading as regards the predictions for the crust-core transition. Even the inclusion of the 4 th order contribution cannot reproduce the results of the exact calculation in the case of SEI. This is due to the sharp rise of the slope of the V thermal computed exactly in the region of low density and low proton fraction that cannot be matched by the Taylor expanded calculations. The transition density and pressure are overestimated under the Taylor series approximation. Therefore, a Taylor expansion of the energy in ANM shall lead to predictions for the properties of the neutron star sensitive to the physical conditions at the crust-core transition region that do not correspond to the actual predictions of the EOS of the model. In this context, the crustal thickness and the crustal fraction of the moment of inertia in neutron stars of different central density are found to take higher values as compared to the exact result when the Taylor series approximation is used. The analytic evaluation of further higher order terms in the Taylor expansion with a finite range interaction becomes a difficult task and prevents one to examine the convergence of the results to the exact prediction. Hence, where possible, the exact analytic expressions should be used in performing studies sensitive to the physical properties in the crust-core transition region. The stiffness parameter γ in nuclear matter has the effect of increasing the transition density with an increase in the incompressibility K(ρ 0 ). Similarly, the transition density also increases with an increase in the value of the symmetry energy parameter E s (ρ 0 ). However, the stiffness of the symmetry energy has an opposite impact and the transition density takes up smaller values when the value of the slope parameter L is larger. All these effects have been examined in the case of the Vela pulsar, for which the lower limit of the crustal fraction of the moment of inertia has been ascertained from the study of the observed glitches. The nuclear matter incompressibility is found to have little influence in the prediction of the Vela pulsar radius. On the other hand, the symmetry energy parameter has a significant influence on the Vela pulsar radius which is predicted to take lower values for a larger symmetry energy parameter. Our present predictions for the crust-core transition density and pressure are based on the thermodynamical method. A natural continuation of this work, to be done in the future, would be to extend the study reported in this paper using the more involved dynamical method.
where the functions J(k i ) and I(k, k i ) with k i = (3π 2 ρ i ) 1/3 (i = n, p) are given by and where x i =k i /Λ (i=n,p), x=k/Λ and Λ = 1/α. The expressions of the energy density for symmetric nuclear matter and pure