A&A 428, 191-197 (2004)
DOI: 10.1051/0004-6361:20041722

Analytical representations of unified equations of state of neutron-star matter

P. Haensel1 - A. Y. Potekhin2,3


1 - N. Copernicus Astronomical Center, Polish Academy of Sciences, Bartycka 18, 00-716 Warszawa, Poland
2 - Ioffe Physico-Technical Institute, Politekhnicheskaya 26, 194021 St. Petersburg, Russia
3 - Isaac Newton Institute of Chile, St. Petersburg Branch, Russia

Received 23 July 2004 / Accepted 7 August 2004

Abstract
Analytical representations are derived for two equations of state (EOSs) of neutron-star matter: FPS and SLy. Each of these EOSs is unified, that is, it describes the crust and the core of a neutron star using the same physical model. Two versions of the EOS parametrization are considered. In the first one, pressure and mass density are given as functions of the baryon density. In the second version, pressure, mass density, and baryon density are given as functions of the pseudo-enthalpy, which makes this representation particularly useful for 2-D calculations of stationary rotating configurations of neutron stars.

Key words: equation of state - dense matter - stars: neutron

   
1 Introduction

The equation of state (EOS) of dense matter is a crucial input for the neutron-star structure calculations. Under standard conditions neutron-star matter is strongly degenerate, and therefore the matter pressure is temperature independent; exceptions are the outermost (a few meters thick) envelopes, newly-born neutron stars, and the envelopes of exploding X-ray bursters. At  $\rho\ga 10^8~\mbox{g~cm$^{-3}$ }$ the EOS is not affected by the magnetic field even as strong as 1014 G and by the temperature $T\la 10^9$ K. Therefore, except for a thin outer envelope (which for the bottom density of  $10^8~{\rm g~cm^{-3}}$ is $\sim$20 m thick and contains only $\sim$ $10^{-9}~{M}_\odot$ of matter), the EOS of neutron-star matter has a one-parameter character. In order to determine the neutron-star structure up to the maximum allowable mass,  $M_{\rm max}$, one has to know the EOS up to a few times  $10^{15}~{\rm g~cm^{-3}}$. While the EOS of the neutron-star crust ( $\rho \la (1{-}2)\times 10^{14}~{\rm g~cm^{-3}}$) is rather well known (e.g., Haensel 2001), the EOS of the liquid core at $\rho\ga 5\times 10^{14}~{\rm g~cm^{-3}}$, which is crucial for determining  $M_{\rm max}$, remains uncertain (e.g., Haensel 2003; Heiselberg & Pandharipande 2000).

The EOS of the crust depends on the crust formation scenario. Two limiting cases are: cold catalyzed matter being at the ground state at a fixed baryon density n, and the accreted crust formed via compression at T<109 K from the thermonuclear ashes of the X-ray bursts in the outer envelopes of accreting neutron stars (Haensel & Zdunik 1990,2003). The EOS of the liquid core does not depend on the formation scenario, and can slightly deviate from the ground state one only when there are deviations from the weak interaction equilibrium - e.g., when rapid matter flow (like in stellar pulsations) is involved. In the present paper we consider the standard case of the ground-state matter.

A "unified EOS'' is obtained in the many-body calculations based on a single effective nuclear Hamiltonian, and is valid in all regions of the neutron star interior. For unified EOSs the transitions between the outer crust and the inner crust, and between the inner crust and the core are obtained as a result of many-body calculation. Alas, up to now only a few models of unified EOSs have been constructed. All other EOSs consist of crust and core segments obtained using different physical models. The crust-core interface has there no physical meaning and both segments are joined using an ad hoc matching procedure. Therefore, neutron-star models based on these EOSs contain a shell with an unphysical EOS. For such "matched EOSs'' it is not possible to study phenomena which are sensitive to the position of the crust-core interface.

In the present paper we consider two unified EOSs, the FPS EOS of Pandharipande & Ravenhall (1989) and the SLy EOS of Douchin & Haensel (2001).

EOSs are usually given in the form of tables. Therefore, in order to use them, one has to employ interpolation between the tabulated points. However, the interpolation procedure is not unique. This introduces ambiguities in the calculated parameters of the neutron star models. Moreover, interpolation should be done respecting exact thermodynamic relations. This turned out to be a particularly serious problem in the high-precision 2-D calculations of models of rapidly rotating neutron stars (Nozawa et al. 1998). In the 3-D calculations of the stationary configurations of a close binary neutron-star system one needs derivatives of pressure with respect to the enthalpy, and tabulated non-polytropic EOSs are not useful in this respect (see, e.g., Gourgoulhon et al. 2001).

In view of the deficiencies and ambiguities inherent in the tabulated EOSs, which are particularly serious in the case of matched EOSs, it is of great interest to derive analytical representations of the EOSs. They have two important advantages over the tabulated ones. First, there is no ambiguity of interpolation and the derivatives can be precisely calculated. Second, these representations can be constructed fulfilling exactly the thermodynamic relations. In this way, analytical EOSs can allow for a very high precision of neutron star structure calculation in the 2-D and 3-D cases.

In the present paper we derive analytical representations of the unified FPS and SLy EOSs. In Sect. 2 we discuss general properties of the EOS used for the calculations of non-rotating and rotating neutron-star models. Analytical representations of the EOSs for non-rotating stars are described in Sect. 3, while Sect. 4 is devoted to the case of rotating stars. Discussion and conclusion are presented in Sect. 5.

   
2 General properties of the one-parameter EOS of neutron-star matter

A given EOS is usually presented in the form of a table containing a grid of calculated values of matter density (full energy density ${\cal E}$ including the rest energies of the matter constituents divided by c2, i.e., $\rho = {\cal E}/c^2$), baryon (number) density n, and pressure P. The EOS  $\lbrace{\rho_i,n_i,P_i\rbrace}$ ( $i=1,\ldots ,N$) is then interpolated between the tabulated points so as to get the one-parameter EOS in the form P=P(n),  $\rho=\rho(n)$. The interpolation involves some degree of indeterminacy (there are many ways of interpolating) and this itself implies some ambiguity as far as the calculated values of the neutron-star parameters (for example, the value of  $M_{\rm max}$) are concerned.

   
2.1 Non-rotating configurations

High precision determination of the baryon number A and gravitational mass M of an equilibrium configuration requires condition of thermodynamic consistency of functions $\rho(n)$ and P(n) to be strictly fulfilled. The first law of thermodynamics in the T=0 limit implies relation


 \begin{displaymath}P(n)=n^2 c^2 \frac{{\rm d}}{{\rm d} n}\left(\frac{\rho}{n}\right), \end{displaymath} (1)

which puts constraints on the interpolation procedure (see, e.g., Haensel & Proszynski 1982); the above relation can be also used in the integral forms,
  
                  $\displaystyle \frac{\rho(n)}{n}$ = $\displaystyle \frac{\rho_{\rm s}}{n_{\rm s}} + \int_{n_{\rm s}}^n {P(n^\prime)\over {n^\prime}^2 c^2} {\rm d}n^\prime,$ (2)
$\displaystyle \ln\left(\frac{n}{n_{\rm s}}\right)$ = $\displaystyle c^2 \int_{\rho_{\rm s}}^\rho \frac{{\rm d}\rho'}{P(\rho') + \rho' c^2},$ (3)

where  $\rho_{\rm s}$ and $n_{\rm s}$ are the values of $\rho $ and n at the neutron-star surface. In the present paper we put  $\rho_{\rm s}$ equal to density of  $^{56}{\rm Fe}$at zero pressure and zero temperature, $\rho_{\rm s}=7.86~{\rm g~cm^{-3}}$. In the outermost neutron-star layers, we fix the value of mass per nucleon as  $m_0=1.66\times 10^{-24}~$g, so that  $n_{\rm s}= \rho_{\rm s}/m_0=4.73494\times 10^{24}~{\rm cm^{-3}}$.

An example of an interpolation recipe which respects Eqs. (1)-(3) and at the same time yields highly smooth functions P(n) and $\rho(n)$ was presented by Swesty (1996). Only a thermodynamically consistent interpolation yields neutron-star models which strictly satisfy the relation connecting the baryon chemical potential $\mu_{\rm b}=c^2 {\rm d}\rho/{\rm d}n_{\rm b}= (\rho c^2 +P)/n$ and the metric function $\Phi(r)$ at a given circumferential radius r. This strict relation stems from the equation for the metric function

 \begin{displaymath}{{\rm d}\Phi\over {\rm d}r}=-{\rho c^2\over \rho c^2 +P}{{\rm d}P\over {\rm d}r}, \end{displaymath} (4)

and can be written as

 \begin{displaymath}\mu_{\rm b}(r) {\rm e}^{\Phi(r)} = \mu_{\rm b}(R) {\rm e}^{\Phi(R)}, \end{displaymath} (5)

where  $\Phi(R)=\sqrt{1-2GM/Rc^2}$. Here, M and R are the total gravitational mass and circumferential stellar radius, respectively. Strictly speaking, if Eq. (5) does not hold, the calculated configuration is not in hydrostatic equilibrium. Alas, this often happens for the neutron-star models calculated using tabulated EOSs with logarithmic interpolation between the tabulated points (see, e.g., Arnett & Bowers 1977; Harrison et al. 1965; Baym et al. 1971). Resulting inconsistencies may seem minor, but they may lead to serious problems if high precision of a simultaneous determination of M and A is concerned; this is the case of the energy release due to a phase transition in the neutron-star core (see, e.g., Haensel & Proszynski 1982). On the contrary, if Eq. (5) is strictly fulfilled, then the constancy of  $\mu_{\rm b}{\rm e}^\Phi$ and the accuracy of calculating stellar parameters are limited only by the numerical precision of the computer code.

For a thermodynamically consistent EOS, Eq. (5) implies the baryon density profile within a static neutron star,

 \begin{displaymath}n(r)=n_{\rm s}\left[{\rho(r)\over \rho_{\rm s}} + {P(r)\over \rho_{\rm s}c^2}\right]{\rm e}^{\Phi(r)-\Phi(R)}, \end{displaymath} (6)

where we have neglected  $P_{\rm s}=P(R,T)$ because the pressure within the atmosphere is small compared with  $\rho_{\rm s} c^2$. Therefore, what one needs to get n(r) in a static neutron star is just the surface baryon density, $n_{\rm s}$; the  $\lbrace{ n_i \rbrace}$ column of the tabulated EOS turns out to be redundant.

   
2.2 Stationary rotating configurations

Rotation breaks the spherical symmetry of the equilibrium configuration. Stationary configurations of a rigidly rotating star of ideal liquid are solutions of 2-D axially symmetric partial differential equations of hydrostatic equilibrium in coordinates r and $\theta$ (rotation of relativistic stars is reviewed by Stergioulas 2003). As far as the EOS is concerned, it is suitable to parametrize it in terms of a dimensionless pseudo-enthalpy

 \begin{displaymath}H(P)\equiv \int_0^P{{\rm d}P^\prime\over \rho (P^\prime)c^2 +P^\prime}\cdot \end{displaymath} (7)

It can be rewritten in terms of the enthalpy per baryon

 \begin{displaymath}h(P)={\rho c^2 +P\over n} \end{displaymath} (8)

as

 \begin{displaymath}H(P)={\rm ln}\left[{h(P)\over h_{\rm s}}\right], \end{displaymath} (9)

where  $h_{\rm s}=\rho_{\rm s}c^2/n_{\rm s}$. The pseudo-enthalpy H vanishes at the stellar surface and increases monotonically towards the stellar center. It is a very useful variable in numerical calculations of the stationary configurations of rotating stars, because the definition of H allows one to easily write down the first integral of motion

 \begin{displaymath}H(r,\theta)+\Phi(r,\theta)-\ln\gamma(r,\theta)= H(0,0)+\Phi(0,0), \end{displaymath} (10)

where $\Phi$ is the metric function and $\gamma$ is the Lorentz factor connecting the two measured values of a physical quantity: that measured in a local observer's reference frame at a point  $(r,\theta,\phi)$and that measured by the distant observer in an inertial reference frame at infinity,

 \begin{displaymath}\gamma = \left( 1- U^2\right)^{-1/2}, \end{displaymath} (11)

where U is the fluid velocity in the azimuthal ($\phi$) direction, as measured by a local observer. Equation (10) is the general relativistic version of the Bernoulli theorem for a stationary ideal fluid flow in gravitational field. Consequently, the most useful parametrization of the EOS for rotating stars is  $\rho=\rho(H)$, P=P(H). Then, the baryon number density is given by an exact analytical formula

 \begin{displaymath}n(H)= n_{\rm s}\left[{\rho(H)\over \rho_{\rm s}}+ {P(H)\over \rho_{\rm s} c^2}\right] {{\rm e}^{-H}}, \end{displaymath} (12)

which is the rotating-star analogue of Eq. (6).

   
3 Analytical representations of the EOS

There are three qualitatively different domains of the interior of a neutron star, separated by phase transition points: the outer crust (consisting of the electrons and atomic nuclei), the inner crust (consisting of the electrons, nuclei, and dripped neutrons), and the core which contains the electrons, neutrons, protons, $\mu^-$-mesons, and possibly $\pi$- and K-mesons, some hyperons, or quark matter. The latter species are contained in the innermost stellar domain called the inner core. In addition, there can be density discontinuities at the interfaces between layers containing different nuclei in the crust. In the fitting, we neglect these small discontinuities and approximate the EOS by fully analytical functions. However, the different character of the EOS in the different domains is reflected by the complexity of the fit, which consists of several fractional-polynomial parts, matched together by virtue of the function

\begin{displaymath}f_0(x) = \frac{1}{{\rm e}^x+1}\cdot \end{displaymath} (13)

We rely on a tabulated unified EOS (FPS[*] or SLy) at  $\rho>5\times10^{10}$ g cm-3. At lower densities, $10^8~\mbox{g~cm$^{-3}$ }\la\rho<5\times10^{10}$ g cm-3, the crustal matter is described by the EOS of Haensel & Pichon (1994, HP94), based on experimental nuclear data, supplemented by the EOS for cold catalyzed matter due to Baym et al. (1971, BPS) at still lower density  $\rho\la10^8$ g cm-3. The lowest-density parts of the tables at  $\rho < 10^5~\mbox{g~cm$^{-3}$ }$ have not been used in the fitting, because at such low density the EOS is no longer one-parametric, but depends also on temperature (see Fig. 1).
  \begin{figure} \par\includegraphics[width=8.3cm,clip]{1722fig1.eps} \end{figure} Figure 1: Neutron-star EOS for non-rotational configurations: BPS (triangles), Haensel & Pichon (1994) (HP94, stars), SLy (dots), OPAL at T=106, 107, and 108 K (dashed lines), the fit (14) (solid line) and the fit modified at low $\rho $ (dotted line).
Open with DEXTER

   
3.1 Non-rotating stars

For non-rotating configurations, we have parametrized the pressure as function of density. Let us denote  $\xi=\log(\rho/{\rm g cm}^{-3})$, $\zeta = \log(P/\mbox{dyn~cm$^{-2}$ })$. Then the parametrization reads

 
                          $\displaystyle \zeta$ = $\displaystyle \frac{a_1+a_2\xi+a_3\xi^3}{1+a_4~\xi}~f_0(a_5(\xi-a_6))$  
    $\displaystyle + (a_7+a_8\xi)~f_0(a_9(a_{10}-\xi))$  
    $\displaystyle + (a_{11}+a_{12}\xi)~f_0(a_{13}(a_{14}-\xi))$  
    $\displaystyle + (a_{15}+a_{16}\xi)~f_0(a_{17}(a_{18}-\xi)).$ (14)

The parameters ai for FPS and SLy EOSs are given in Table 1. The typical fit error of P is 1-2% (for $\xi\ga5$). The maximum error is determined by the jumps near the phase transitions in the tabulated EOSs, which are smoothed by the fit (14). For FPS, the maximum error is 3.6% at $\xi=14.22$(crust-core interface). For SLy, the maximum error is 2.9% at $\xi=8.42$(62Ni-64Ni interface in the HP94 part of the table).

Figure 1 shows $\log P$ against $\log\rho$for a tabulated EOS (symbols) and the corresponding fit (solid line). Triangles correspond to BPS, stars to HP94, and dots to SLy data. By construction, the fit is accurate at $\rho\ga10^5$ g cm-3. As stated above, at lower density the EOS becomes temperature-dependent (for T values typical for neutron star envelopes). This is illustrated by the dashed lines, that show the OPAL EOS (Rogers et al. 1996) for T=106, 107, and 108 K.[*] However, a reasonable continuation of the fit to lower densities can be constructed by a simple interpolation. For instance, the dotted line in Fig. 1 corresponds to $P=P_{\rm fit}+P_0$, where  $P_{\rm fit}=10^\zeta$ is given by Eq. (14) (where we always assume $\xi>0$), and  $P_0=3.5\times10^{14}~\rho$ approximates the OPAL EOS near  $\rho\sim\rho_{\rm s}$ at T=107 K (here P is in dyn cm-2 and $\rho $ in g cm-3).

Table 1: Parameters of the fit (14).


  \begin{figure} \par\includegraphics[width=8.3cm,clip]{1722fig2.eps} \end{figure} Figure 2: Comparison of the data and fits for SLy and FPS EOSs for non-rotational configurations. Upper panel: rarefied tabular data (symbols) and the fit (14) (lines); lower panel: relative difference between the data and fit. Filled dots and solid line: SLy; open circles and dot-dashed line: FPS (triangles and stars on the upper panel are BPS and HP94 data at  $\rho <5\times 10^{10}$ g cm-3).
Open with DEXTER

Table 2: Parameters of the fits (15) and (16).

In Fig. 2 we compare the FPS and SLy EOSs. Symbols on the upper panel show the data (triangles, stars, dots, and open circles for BPS, HP94, SLy, and FPS, respectively) and lines show the fit (solid for SLy and dot-dashed for FPS). In order to make the differences between the data and fits and between SLy and FPS EOSs visible, we plot the difference  $\log P - 1.4 \log \rho$, where P is in dyn cm-2 and $\rho $ in g cm-3. The lower panel shows the relative difference between the tabulated and fitted EOSs (solid and dot-dashed lines for SLy and FPS, respectively). It illustrates the accuracy of the fit (14).

Now, $n(\rho)$ can be easily obtained from Eq. (3). Doing this and substituting $P(\rho')$ in the integrand of Eq. (3) from Eq. (14), we recover the original tabular values with maximum difference <0.4% for FPS and <0.12% for SLy.

In some applications, it may be convenient to use n as independent variable, and treat $\rho $ and P as functions of n. For this purpose one can use the following fit:

 
$\displaystyle \frac{\rho}{n m_0}$ = $\displaystyle 1+ \frac{p_1 n^{p_2}+p_3 n^{p_4}}{(1+p_5 n)^2}~ f_0(-{p_6}(\log n + p_7))$ $\displaystyle + \frac{n}{8\times10^{-6} + 2.1~n^{0.585}}~ f_0({p_6}(\log n + p_7)),$ (15)

where n is in fm-3. The inverse fit $n(\rho)$ is given by the formula
 
$\displaystyle \frac{x}{n}$ = $\displaystyle 1+ \frac{q_1 x^{q_2}+q_3 x^{q_4}}{(1+q_5 x)^3}~ f_0({q_6}(q_7 - \log \rho))$ $\displaystyle + \frac{x}{8\times10^{-6} + 2.1~x^{0.585}}~ f_0({q_6}(\log \rho - q_7)),$ (16)

where  $x=\rho/m_0$ and $\rho $ is in g cm-3. Coefficients pi and qi of the fits (15) and (16) are given in Table 2. The fractional fit residual in $\rho $ and n varies from $\la$10-9 at $\rho\la10$ g cm-3 to a fraction of percent near the maximum  $\rho\sim10^{16}$ g cm-3, while the difference  $(\rho-n m_0)$is approximated with the typical accuracy of a few percent. In Fig. 3 (upper panel) we plot the ratio  $\rho/(n m_0)$against n for the two EOSs. We see that in the crust, at $n\la0.1$ fm-3, $\rho $ is to a good accuracy proportional to n. The lower panel, which shows the relative difference between fitted and tabulated values of $\rho(n)$, confirms the accuracy of Eq. (15).


  \begin{figure} \par\includegraphics[width=8.3cm,clip]{1722fig3.eps} \end{figure} Figure 3: Dependence of $\rho $ on n. Upper panel: rarefied data (symbols) and the fit (15) (lines); lower panel: relative difference between the data and fit. Filled dots and solid line: SLy; open circles and dot-dashed line: FPS.
Open with DEXTER

It should be stressed that thermodynamics requires Eq. (1) to be satisfied exactly. To achieve this, one should not totally rely on the fits (15) and (16); otherwise thermodynamic consistency will be violated on the scale of these fits' errors (a fraction of percent). So, if $\rho $ is used as an input, then $n(\rho)$ should be calculated from Eq. (3). Alternatively, if the input is n, then, having calculated  $\rho_{\rm fit}(n)$ from Eq. (15) and  $P(n)=P(\rho_{\rm fit}(n))$ from Eq. (14), one should refine $\rho(n)$ using relation (2).

   
3.2 Rotating stars

As explained in Sect. 2.2, for rotating stars it is most useful to parametrize density and pressure as functions of the pseudo-enthalpy H, which can be written in terms of the enthalpy per baryon h. Let us define  $\eta\equiv h/m_0 c^2-1$. In view of relation (8), the function which we intend to parametrize, $\xi(\eta)$, is not independent of the function  $\zeta(\xi)$parametrized by Eq. (14). In order to fulfill Eq. (8) as closely as possible, we first calculate $\eta(\xi)$ using Eqs. (14) and (8), and then find the inverse fit $\xi(\eta)$. The best fit reads:

 
                               $\displaystyle \xi$ = $\displaystyle \left[ b_1+b_2 \log\eta+\frac{b_3\eta^{b_4}}{1+b_5\eta} \right] f_0(b_6(\log\eta-b_7))$  
    $\displaystyle + \frac{b_8+b_9\log\eta+(b_{10}+b_{11}\log\eta)(b_{12}\eta)^7}{ 1+b_{13}\eta+(b_{12}\eta)^7}$  
    $\displaystyle \times f_0(b_6(b_7-\log\eta)) + b_{14}~f_0(b_{15}(b_{16}-\log\eta)).$ (17)

Where the parameters bi are given in Table 3. The comparison of the fit and the data is presented in Fig. 4. The typical fit error of $\rho $, according to Eq. (17), is about 1% at  $\eta\ga10^{-7}$(corresponding to $\xi\ga 3$), and the maximum error <4% occurs near the neutron drip and crust-core phase transitions.

Table 3: Parameters of the fit (17).


  \begin{figure} \par\includegraphics[width=8.3cm,clip]{1722fig4.eps} \end{figure} Figure 4: SLy and FPS EOS for rotational configurations. Upper panel: rarefied data to be fitted, calculated according to Eq. (14) (symbols) and the fit (17) (lines); lower panel: relative difference between the data and fit. Filled dots and solid line: SLy; open circles and dot-dashed line: FPS.
Open with DEXTER

When used in combination, the fits (14) and (17), together with Eq. (3) or Eq. (16) give the parametrizations of $\rho(H)$, P(H), and n(H) needed for calculations of the stationary rotating configurations. In this case, the function  $P(H)=P(\rho(H))$ obtained using Eqs. (17) and (14), reproduces the tabular values with a typical discrepancy of 1-2%, with a maximum within 10% near the crust-core boundary.

The remark on the thermodynamic consistency, made at the end of Sect. 3.1, applies also here: one should refine either n or $\rho $ fitted values, using the exact relations (3) or (2).

   
4 Adiabatic index

An important dimensionless parameter characterizing the stiffness of the EOS at given density is the adiabatic index, defined by

\begin{displaymath}\Gamma=\frac{n}{P}~\frac{{\rm d}P}{{\rm d}n} = \left[ 1+ \fr... ... c^2} \right] \frac{\rho}{P}\frac{{\rm d}P}{{\rm d}\rho}\cdot \end{displaymath} (18)

Using our fit (14), we obtain the analytical expression
 
                                             $\displaystyle \frac{\rho}{P}\frac{{\rm d}P}{{\rm d}\rho} = \frac{{\rm d}\zeta}{{\rm d}\xi} = \bigg[ \frac{ a_2-a_1 a_4 + 3a_3\xi^2+2a_3 a_4 \xi^3}{(1+a_4\xi)^2}$  
    $\displaystyle \quad-~a_5~\frac{1_1+a_2\xi+a_3\xi^3}{1+a_4\xi}~f_0(a_5(a_6-\xi))\bigg] f_0(a_5(\xi-a_6))$  
    $\displaystyle \quad+ \sum_{i=2}^4 f_0(a_{4i+1}(a_{4i+2}-\xi)) \big[ a_{4i}$  
    $\displaystyle \quad+~a_{4i+1} (a_{4i-1}+a_{4i}\xi)f_0(a_{4i+1}(\xi-a_{4i+2})) \big].$ (19)

Different regions of neutron-star interior are characterized by distinct behavior of $\Gamma $, displayed in Fig. 5. Precise values of $\Gamma $ calculated by Douchin & Haensel (2001) are shown by dotted line, and the analytical approximations according to Eqs. (14), (19) by solid line (for SLy) and dot-dashed line (for FPS).


  \begin{figure} \par\includegraphics[width=8.3cm,clip]{1722fig5.eps} \end{figure} Figure 5: Adiabatic index $\Gamma $ for SLy and FPS EOSs. Solid line: analytical approximation (SLy); dotted line: precise values (SLy); dot-dashed line: analytical approximation (FPS).
Open with DEXTER

In the outer crust, the value of $\Gamma $ depends quite weakly on density. At  $10^8~\mbox{g~cm$^{-3}$ }\la\rho\la3\times10^{11}~\mbox{g~cm$^{-3}$ }$, this value would be  $\Gamma\approx4/3$ if A,Z values were fixed, because in this case P is mainly determined by the pressure of ultrarelativistic electron gas which behaves as $\propto$ $(Z\rho/A)^{4/3}$. For example,  $\Gamma\approx4/3$ within each shell with constant A and Z of the EOS derived by Haensel & Pichon (1994). However, the compressible liquid drop model used by Douchin & Haensel (2001) effectively smoothes the discontinuities caused by transitions from one to another (A,Z)-species with increasing density, which leads to an effective continuous increase of the A/Z ratio and corresponding decrease of $\Gamma $, seen in Fig. 5.

A dramatic drop in $\Gamma $ occurs at neutron drip threshold, which corresponds to strong softening of the EOS. The analytical expression (19) somewhat smoothes this drop. The behavior of $\Gamma $ in the inner crust results from an interplay of several factors, with stiffening due to interaction between dripped neutrons, a softening effect of neutron-gas - nuclear-matter coexistence, and the softening Coulomb contribution.

At the crust-core interface, matter strongly stiffens, and $\Gamma $jumps from $\approx$1.7 to $\approx$2.2, which results from disappearance of nuclei. The analytical approximation smoothes this jump also, though reflects the stiffening. This approximation is also smooth across a small discontinuous drop of $\Gamma $at  $\rho\approx2\times10^{14}~\mbox{g~cm$^{-3}$ }$, where muons start to replace a part of ultrarelativistic electrons. However, the electrons and muons give only minor contribution to the pressure (and therefore behavior of $\Gamma $) in the core, because the main contribution comes from interactions between nucleons.

   
5 Summary, discussion, and conclusion

Analytical representations of the EOSs used in the modern 2-D and 3-D simulations of neutron star dynamics have many important advantages over the tabulated EOSs. Analytical EOSs do not require any interpolations, are thermodynamically consistent, and allow for a very high precision of calculations. In the present paper we constructed analytical representations, in terms of the continuous and differentiable functions of a single chosen variable, of the SLy and FPS EOSs. For these analytical representations, the thermodynamic relations are exactly satisfied at any point. Two choices of the independent variable were considered. The first one is $\rho $; function $P(\rho)$ given by Eq. (14) fits the original tables in the density interval  $10^{5}~{\rm g~cm^{-3}}<\rho<10^{16}~{\rm g~cm^{-3}}$within typical error of 1-2%. Function $n(\rho)$ can be calculated either from Eq. (3) to satisfy exactly the first law of thermodynamics, or from the fit (16) with a typical error $\sim$0.1%. A variant which ensures the same accuracy is to choose nas an independent variable and calculate $\rho(n)$ from the fit (15) and $P(\rho)$ from Eq. (14). The other choice of the independent variable is to use the pseudo-enthalpy H. This choice is particularly advantageous for applications to 2-D and 3-D numerical simulations of neutron star dynamics, such as rotation and inspiraling stage of the evolution of relativistic neutron-star - neutron-star binary. We represented both EOSs by the continuous and differentiable functions P(H), $\rho(H)$, and n(H), where $\rho(H)$ is given by Eq. (17) with typical accuracy within a few percent, while P(H) and n(H)are calculated from the functions $P(\rho)$ and $n(\rho)$, respectively. Differentiation of $P(\rho)$ then yields analytical representations of the adiabatic index for the SLy and FPS EOSs; this quantity is important, for example, for numerical simulations of dynamics of the neutron-star - neutron-star system at the inspiral phase.

The quality of our analytical representations of the EOSs was tested by evaluation of the virial identities GRV2 (Bonazzola 1973; Bonazzola & Gourgoulhon 1994) and GRV3 (Gourgoulhon & Bonazzola 1994) in the numerical simulations of the 2-D stationary rotation of neutron stars. GRV2 and GRV3 are integral identities which must be satisfied by a stationary solution of the Einstein equations and which are not imposed in the numerical procedure (see Nozawa et al. 1998 for the details of computation of GRV2 and GRV3). For rotating configurations we get ${\rm GRV2,~GRV3}\sim 10^{-6}{-}10^{-5}$, which is excellent.

The subroutines for the numerical applications of analytical EOSs (Fortran and C++ versions) can be downloaded from the public domain http://www.ioffe.ru/astro/NSG/NSEOS/

Acknowledgements
We are very grateful to M. Bejger for performing the 2-D calculations for rapidly rotating neutron stars with our EOSs, and to M. Shibata for pointing out to significant misprints in Eqs. (15) and (16) of the preliminary version of this paper. The work of P.H. was partially supported by the KBN grant No. 1-P03D-008-27. The work of A.P. was partially supported by the RFBR grants 02-02-17668 and 03-07-90200 and the Russian Leading Scientific Schools grant 1115.2003.2.

References

 
Copyright ESO 2004