# Non-linear Evolution of Rotating Relativistic Stars

## 1 Introduction

The numerical evolution of neutron stars in full General Relativity has been the focus of many research groups in recent years . So far, these studies have been limited to initially non-rotating stars. However, the numerical investigation of many interesting astrophysical applications, such as the rotational evolution of proto-neutron stars and merged neutron stars or the simulation of gravitational radiation from unstable pulsation modes, requires the ability of accurate long-term evolutions of rapidly rotating stars. We thus present here the first study of hydrodynamical evolutions of rotating neutron stars in the approximation of a static spacetime. This approximation allows us to evolve relativistic matter for a much longer time than present coupled spacetime plus hydrodynamical evolution codes. Since the pulsations of neutron stars are mainly a hydrodynamical process, the exclusion of the spacetime dynamics has only a limited effect and allows for qualitative conclusions to be drawn.

The rotational evolution of neutron stars can be affected by several instabilities (see  for a recent review). If hot protoneutron stars are rapidly rotating, they can undergo a dynamical bar-mode instability . When the neutron star has cooled to about $`10^{10}`$K after its formation, it can be subject to the Chandrasekhar-Friedman-Schutz instability  and it becomes an important source of gravitational waves. It was recently found that the $`l=m`$ $`r`$-mode has the shortest growth time of the instability  and it can transform a rapidly rotating newly-born neutron star to a Crab-like slowly-rotating pulsar within about a year after its formation . In this model, there are two important questions still to be answered : What is the maximum amplitude that an unstable $`r`$-mode can reach (limited by nonlinear saturation) and is there any transfer of energy to other stable or unstable modes via non-linear couplings? Such questions cannot be answered by computations of normal modes of the linearized pulsation equations, but require non-linear effects to be taken into account. We therefore need to develop the capability of full non-linear numerical evolutions of rotating stars in General Relativity.

Our present 2-D (axisymmetric) code uses high-resolution shock-capturing (HRSC) finite-difference schemes for the numerical integration of the general relativistic hydrodynamic equations  (see  for a recent review of applications of HRSC schemes in relativistic hydrodynamics). In a similar context to the one presented here let us note that such schemes have been succesfully used before in the study of the numerical evolution and gravitational collapse of non-rotating neutron stars in 1-D . An alternative approach, based on pseudospectral methods, has been presented in .

Using our code in 1-D time-evolutions we can accurately identify specific normal modes of pulsation. In 2-D the code is suitable for the evolution of rotating stars, with the additional complication of having to pay special attention to an angular momentum-loss at the (non-spherical) surface of the star, as we will show below.

## 2 Initial Configurations

Our initial models are fully relativistic, stationary and axisymmetric configurations, rotating with uniform angular velocity $`\mathrm{\Omega }`$. The metric in quasi-isotropic coordinates  is

$`ds^2=e^{2\nu }dt^2+B^2e^{2\nu }r^2\mathrm{sin}^2\theta (d\varphi \omega dt)^2+e^{2\alpha }(dr^2+r^2d\theta ^2),`$ (1)   
where $`\nu `$, $`B`$, $`\alpha `$ and $`\omega `$ are metric functions (gravitational units are implied). In the non-rotating limit the above metric reduces to the metric of spherical relativistic stars in isotropic coordinates.

We assume a perfect fluid, zero-temperature equation of state (EOS), for which the energy density is a function of pressure only. The following relativistic generalization of the Newtonian polytropic EOS is chosen:

$`p`$ $`=`$ $`K\rho _0^{1+1/N}`$ (2)   
$`ϵ`$ $`=`$ $`\rho _0+Np,`$ (3)   
where $`p`$ is the pressure, $`ϵ`$ is the energy density, $`\rho _0`$ is the rest-mass density, $`K`$ is the polytropic constant and $`N`$ is the polytropic exponent.

The initial equilibrium models are computed using a numerical code by Stergioulas & Friedman  which follows the Komatsu, Eriguchi & Hatchisu  method (as modified in ) with some changes for improved accuracy (see  for a comparison with other codes). The code is freely available and can be downloaded from the following URL address: http://www.gravity.phys.uwm.edu/Code/rns.

## 3 Relativistic Hydrodynamic Equations

The equations of (ideal) relativistic hydrodynamics are obtained from the local conservation laws of density current, $`J^\mu `$ and stress-energy, $`T^{\mu \nu }`$

$`_\mu J^\mu `$ $`=`$ $`0`$ (4)   
$`_\mu T^{\mu \nu }`$ $`=`$ $`0`$ (5)   
with

$`J^\mu `$ $`=`$ $`\rho _0u^\mu `$ (6)   
$`T^{\mu \nu }`$ $`=`$ $`\rho _0hu^\mu u^\nu +pg^{\mu \nu },`$ (7)   
for a general EOS $`p=p(\rho ,\epsilon )`$. This choice of the stress-energy tensor limits our study to perfect fluids.

In the previous expressions $`_\mu `$ is the covariant derivative, $`u^\mu `$ is the fluid 4-velocity and $`h`$ is the specific enthalpy

$`h=1+\epsilon +{\displaystyle \frac{p}{\rho _0}}`$ (8)   
with $`\epsilon `$ being the specific internal energy, related to the energy density $`ϵ`$ by

$`\epsilon ={\displaystyle \frac{ϵ}{\rho _0}}1.`$ (9)   
With an appropriate choice of matter fields the equations of relativistic hydrodynamics constitute a (non-strictly) hyperbolic system and can be written in a flux conservative form, as was first shown in  for the one-dimensional case. The knowledge of the characteristic fields of the system allows the numerical integration to be performed by means of advanced high-resolution shock-capturing (HRSC) schemes,using approximate Riemann solvers (Godunov-type methods). The multidimensional case was studied in , within the framework of the 3+1 formulation. Further extensions of this work to account for dynamical spacetimes, described by the full set of Einstein’s non-vacuum equations, can be found in . Fully covariant formulations of the hydrodynamic equations (i.e., not restricted to spacelike approaches) and also adapted to Godunov-type methods, are presented in .

In the present work we use the hydrodynamic equations as formulated in . Specializing for the metric given by Eq. (1), the 3+1 quantities read

$`\stackrel{~}{\alpha }`$ $`=`$ $`e^\nu `$ (10)   
$`\beta _\varphi `$ $`=`$ $`\omega B^2e^{2\nu }r^2\mathrm{sin}^2\theta `$ (11)   
$`\gamma _{rr}`$ $`=`$ $`e^{2\alpha }`$ (12)   
$`\gamma _{\theta \theta }`$ $`=`$ $`r^2e^{2\alpha }`$ (13)   
$`\gamma _{\varphi \varphi }`$ $`=`$ $`B^2e^{2\nu }r^2\mathrm{sin}^2\theta `$ (14)   
where $`\stackrel{~}{\alpha }`$ is the lapse function (the tilde is used to avoid confussion with the metric potential $`\alpha `$) and $`\beta _\varphi `$ is the azimuthal shift.

The hydrodynamic equations are written as a first-order flux conservative system of the form

$`{\displaystyle \frac{𝐮}{t}}+{\displaystyle \frac{\stackrel{~}{\alpha }𝐟^r}{r}}+{\displaystyle \frac{\stackrel{~}{\alpha }𝐟^\theta }{\theta }}=𝐬`$ (15)   
where $`𝐮,𝐟^r,𝐟^\theta `$ and $`𝐬`$ are, respectively, the state vector of evolved quantities, the radial and polar fluxes and the source terms. More precisely, they take the form

$`𝐮`$ $`=`$ $`(D,S_r,S_\theta ,S_\varphi ,\tau )`$ (16)   
$`𝐟^r`$ $`=`$ $`(Dv^r,S_rv^r+p,S_\theta v^r,S_\varphi v^r,(\tau +p)v^r)`$ (17)   
$`𝐟^\theta `$ $`=`$ $`(Dv^\theta ,S_rv^\theta ,S_\theta v^\theta +p,S_\varphi v^\theta ,(\tau +p)v^\theta ).`$ (18)   
The source terms can be decomposed in the following way

$`𝐬=\stackrel{~}{\alpha }𝐬^{}\stackrel{~}{\alpha }𝐟^r{\displaystyle \frac{\mathrm{log}\sqrt{\gamma }}{r}}\stackrel{~}{\alpha }𝐟^\theta {\displaystyle \frac{\mathrm{log}\sqrt{\gamma }}{\theta }}`$ (19)   
with $`\gamma =det\gamma _{ij}`$ and

$`𝐬^{}=(0,T^{\mu \nu }\left[{\displaystyle \frac{g_{\nu j}}{x^\mu }}\mathrm{\Gamma }_{\mu \nu }^\delta g_{\delta j}\right],\stackrel{~}{\alpha }\left[T^{\mu t}{\displaystyle \frac{\mathrm{log}\stackrel{~}{\alpha }}{x^\mu }}T^{\mu \nu }\mathrm{\Gamma }_{\mu \nu }^t\right])`$ (20)   
with $`j=r,\theta ,\varphi `$. The definitions of the evolved quantities in terms of the “primitive” variables ($`\rho ,v_j,\epsilon `$) are

$`D`$ $`=`$ $`\rho _0W`$ (21)   
$`S_j`$ $`=`$ $`\rho _0hW^2v_j`$ (22)   
$`\tau `$ $`=`$ $`\rho _0hW^2pD`$ (23)   
where $`W`$ is the relativistic Lorentz factor

$`W\stackrel{~}{\alpha }u^t={\displaystyle \frac{1}{\sqrt{1v^2}}}`$ (24)   
with $`v^2=\gamma _{ij}v^iv^j`$. The 3-velocity components are obtained from the spatial components of the 4-velocity in the following way

$`v^i={\displaystyle \frac{u^i}{W}}+{\displaystyle \frac{\beta ^i}{\stackrel{~}{\alpha }}}.`$ (25)   
Explicit expressions for the non-vanishing Christoffel symbols for metric (1), appearing in the source terms of the hydrodynamic equations, are presented in .

## 4 Numerical Methods

As stated before, our numerical integration of system (15) is based on Godunov-type methods (also known as HRSC schemes). In a HRSC scheme, the knowledge of the characteristic fields (eigenvalues) of the equations, together with the corresponding eigenvectors, allows for accurate integrations, by means of either exact or approximate Riemann solvers, along the fluid characteristics. These solvers, which constitute the kernel of our numerical algorithm, compute, at every interface of the numerical grid, the solution of local Riemann problems (i.e., the simplest initial value problem with discontinuous initial data). Hence, HRSC schemes automatically guarantee that physical discontinuities appearing in the solution, e.g., shock waves, are treated consistently (the shock-capturing property). HRSC schemes are also known for giving stable and sharp discrete shock profiles. They have also a high order of accuracy, typically second order or more, in smooth regions of the solution.

We perform the time update of system (15) according to the following conservative algorithm:

$`𝐮_{i,j}^{n+1}=𝐮_{i,j}^n`$ $``$ $`{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{\Delta }r}}(\widehat{𝐟}_{i+1/2,j}\widehat{𝐟}_{i1/2,j})`$ (26)   
$``$ $`{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{\Delta }\theta }}(\widehat{𝐠}_{i,j+1/2}\widehat{𝐠}_{i,j1/2})+\mathrm{\Delta }t𝐬_{i,j}.`$   
Index $`n`$ represents the time level and the time (space) discretization interval is indicated by $`\mathrm{\Delta }t`$ ($`\mathrm{\Delta }r,\mathrm{\Delta }\theta `$). The “hat” in the fluxes is used to denote the so-called numerical fluxes which, in a HRSC scheme, are computed according to some generic flux-formula, of the following functional form (suppressing index $`j`$):

$`\widehat{𝐟}_{i\pm \frac{1}{2}}={\displaystyle \frac{1}{2}}\left(𝐟(𝐮_{i\pm \frac{1}{2}}^L)+𝐟(𝐮_{i\pm \frac{1}{2}}^R){\displaystyle \underset{\alpha =1}{\overset{5}{}}}\stackrel{~}{\lambda }_\alpha \mathrm{\Delta }\stackrel{~}{\omega }_\alpha \stackrel{~}{r}_\alpha \right).`$ (27)   
Notice that the numerical flux is computed at cell interfaces ($`i\pm 1/2`$). Indices $`L`$ and $`R`$ indicate the left and right sides of a given interface. Quantities $`\lambda `$, $`\mathrm{\Delta }\omega `$ and $`r`$ denote the eigenvalues, the jump of the characteristic variables and the eigenvectors, respectively, computed at the cell interfaces according to some suitable average of the state vector variables. Generic expressions for the characteristic speeds and eigenfields can be found in . Our code has the ability of using different approximate Riemann solvers: the Roe solver , widely employed in fluid dynamics simulations, with arithmetically averaged states and the Marquina solver , which has been extended to Relativity in . The computations presented here were obtained using Marquina’s scheme.

A technical remark: the equilibrium star is supplemented by a low-density uniform atmosphere, which is necessary for computing non-singular solutions of the hydrodynamic equations everywhere in the computational domain. After each time-step we reset the atmosphere’s density and pressure to their initial values, avoiding unwanted accretion of matter onto the star. The influence of the atmosphere is thus restricted to the surface grid-cells.

## 5 Pulsations of Non-rotating Stars

Since our code uses spherical polar coordinates, it can also be employed to study the evolution of non-rotating stars in 1-D. In the evolution of initially static non-rotating stars, we observe the following properties (note that our numerical grid is Eulerian):

1. Small-amplitude radial pulsations are triggered by the truncation errors of the finite-differencing scheme.
2. The radial pulsations are dominated by a set of discrete frequencies, which correspond to the normal modes of pulsation of the star.
3. The numerical viscosity of the finite-difference scheme damps the pulsations and the damping is stronger for the higher frequency modes.
4. The presence of a constant density atmosphere affects the finite differencing at the surface grid-cells, which increases the numerical damping of pulsations and also causes a continuous but very small drifting of the density distribution.

The initial amplitude of the radial pulsations and the small drift in density converge to zero at a second order rate with increasing resolution. The value of the density in the atmosphere region has a large effect on the damping of the pulsations. If it is too large, the damping is strong. To minimize this effect, we typically set the density of the atmosphere equal to $`10^6`$ times the density of the last grid point inside the star.

### 5.1 1-D evolutions

We study the numerical evolution of a nonrotating $`N=1.5`$ relativistic polytrope with $`M/R=0.056`$. The star is immediately set into radial pulsation, triggered by the finite difference truncation errors. The time-evolution of the radial velocity $`v_r`$, computed at 25% of the star radius using a grid of 400 zones, is shown in Fig. 1. The vertical axis is dimensionless ($`c=G=M_{}=1`$). The radial velocity is initially a very complex function of time. As we will show next, the pulsation consists mainly of a superposition of normal modes of oscillation of the fluid. The high frequency normal modes are damped quickly and after 20ms the star pulsates mostly in lowest frequency modes. Because these oscillations are caused only by the truncation errors, the magnitude of the radial velocity is extremely small. As shown in Fig. 1 it is only a few times larger than the non-zero residual velocity around which the star is oscillating. This residual velocity converges to zero as second order with increased resolution.

The small-amplitude radial pulsations in the non-linear, fixed spacetime evolutions correspond to linear normal modes of pulsation in the relativistic Cowling approximation, in which perturbations of the spacetime are ignored. A Fourier transform of the density or radial velocity time-evolution can be used to identify the normal mode frequencies. Fig. 2 shows the Fourier transform of the radial velocity evolution shown in Fig. 1. The normal mode frequencies stand out as sharp peaks on a continuous background. The width of the peaks increases with frequency. The frequencies of radial pulsations identified from Fig. 2 are shown in Table 1.

To compare the obtained frequencies to linear normal mode frequencies, we use a different code that solves the linearized relativistic pulsation equations for the stellar fluid, in the Cowling approximation , as an eigenvalue problem. In Table 1 we present the results of this comparison. The typical agreement between frequencies computed by the two methods is better than 0.5% for the fundamental $`F`$-mode and the lowest frequency harmonics $`H_1H_4`$ and better than 0.8% for the the higher harmonics $`H_5H_9`$. This is a strong test for the accuracy of the evolution code and our results can be used as a testbed computation for other relativistic multi-dimensional evolution codes.

### 5.2 2-D evolutions

In a similar way, small-amplitude non-radial pulsations can be studied with the present evolution code and the obtained frequencies can be compared to perturbation results. We find that the truncation errors of the finite difference scheme do not excite non-radial pulsations to a sufficiently large amplitude compared to the amplitude of radial pulsations, so that one cannot identify them accurately in a Fourier transform. Instead, one has to perturb the initial configuration, using an appropriate eigenfunction for each nonradial angular index $`l`$. Such a perturbation can be constructed using the eigenfunctions of linear pulsation modes, computed with the perturbation code in the Cowling approximation. The frequencies of the non-radial modes are then found from a Fourier transform of the time-evolution of the velocity component $`v_\theta `$.

Table 2 shows a similar comparison as in Table 1 for the quadrupole ($`l=2`$) pulsations of the same $`N=1.5`$ relativistic polytrope. Since the non-radial modes have to be computed on a 2-D grid, we cannot use resolutions as high as in the 1-D computations. For a small grid-size of $`80\times 80`$ zones and a total evolution time of 6ms, the agreement between frequencies computed by the two methods is better than 1.4% for the fundamental $`f`$-mode and the $`p`$-modes $`p_1p_6`$. For this grid-size, frequencies higher than the $`p_6`$ mode could not be computed accurately, because the grid is to coarse to resolve their eigenfunctions (higher harmonic eigenfunctions have a larger number of nodes in the radial direction).

## 6 Rotating Stars

We now turn to the evolution of initially stationary, uniformly rotating neutron stars. In these evolutions, we observe the same qualitative properties as for non-rotating stars (section 5) and an additional important property: the angular momentum of the star is not conserved at the surface layer. This is due to the fact that the velocity component $`v_\varphi `$ of the fluid has a maximum at the surface, while the numerical scheme (although second-order accurate in smooth regions of the solution) is only first-order accurate at local extrema. Moreover, the code evolves the relativistic momenta, $`S_i`$, and the velocity components (as well as the rest of “primitive” variables) must be recovered through a root finding procedure which involves dividing by the density. At the surface of the star (where the density is very small) this contributes to obtaining less than second-order accuracy.

A representative example of the evolution of a rotating star is presented in Fig. 3, which shows the evolution, at different times, of the velocity component $`v_\varphi `$. The star is again a $`N=1.5`$ polytrope with the same central density as the non-rotating star presented in section 5 and rotating at 74% of the mass-shedding limit at same central density. The evolution was for one rotation period on a $`96\times 60`$ grid. The vertical axis is dimensionless ($`c=G=M_{}=1`$). The figure shows that the $`\varphi `$-velocity in the interior of the star remains close to its initial value, while it decreases as a function of time in the outer layers. We find that this is a generic property of the present numerical scheme for any rotation rate and for any grid-size. By comparing evolutions with different grid sizes, we verified that the loss of angular momentum at the surface improves as first-order with resolution, while the evolution of the $`\varphi `$-velocity in the interior is second-order accurate. However, as the evolution proceeds in time, the first-order surface effect gradually affects the interior of the star.

## 7 Quasi-radial Modes of Rotating Stars

As a first application of our code, we compute quasi-radial modes (i.e. modes that in the non-rotating limit reduce to radial modes) of rapidly rotating relativistic stars in the Cowling approximation. Previously, these modes have been computed for fully relativistic stars only in the slow-rotation limit (but without the assumption of a fixed spacetime) by Hartle & Friedman  (see also ). We compute the three lowest-frequency quasi-radial modes for a sequence of rotating stars of same central density. The non-rotating member of the sequence is the non-rotating star of section 5. Table 3 and Fig. 4 show our results for a low resolution grid of $`100\times 80`$ zones (note that our computational grid assumes equatorial plane symmetry). For this resolution we estimate the accuracy of the frequencies to be of the order of $`12`$%.

For the sequence of stars considered here, the frequencies of the quasi-radial modes decrease with increasing rotation rate. This agrees with previous slow-rotation computations which predict a decrease as $`\mathrm{\Omega }^2`$, where $`\mathrm{\Omega }`$ is the angular velocity of the star. For fast rotation, the change in the frequencies of quasi-radial modes is affected by higher order terms in $`\mathrm{\Omega }`$, because of the large deformation of the equilibrium star. Also, for rapidly rotating stars the quasi-radial mode frequencies are more “closely packed” than in non-rotating stars.

## 8 Discussion

Our axisymmetric relativistic hydrodynamical code is capable to evolve rapidly rotating stars in a fixed spacetime. We find that, for non-rotating stars, small amplitude oscillations have frequencies that agree with linear normal mode frequencies in the Cowling approximation and we compute the quasi-radial modes of rapidly rotating stars.

Modern HRSC numerical schemes (as the ones used in our code), satisfying the “total variation diminishing” (TVD) property , are second-order accurate in smooth regions of the flow, but only first-order accurate at local extrema. In our rotating stars runs we find that this results in a loss of angular momentum of the surface layers of the star, which gradually also affects the interior of the star. This angular momentum loss only vanishes as first-order with incresing resolution and we thus conclude that for accurate long-term evolutions of rotating neutron stars it is essential to use rather fine grids. Furthermore, to reduce the computational cost, one could use surface-adapted coordinates or fixed-mesh refinement. It would also be interesting to see whether the loss of angular momentum per rotation period will be significantly smaller in a frame co-rotating with the star. An alternative solution to this problem, which we plan to investigate, could be the use of “essentially non-oscillatory” (ENO) schemes, which maintain high-order of accuracy even at local extrema .

All previous considerations are important for the study of the non-linear dynamics of unstable toroidal oscillations ($`r`$-modes) in 3-D, which have a long growth time and thus require highly accurate long-term evolutions.

## Acknowledgements

We thank John L. Friedman, Curt Cutler, Philippos Papadopoulos and Tom Goodale for helpful discussions. We also thank S. Yoshida for sending us for comparison unpublished results on quasi-radial modes of rotating stars in the Cowling approximation, computed with a linear perturbation code. J.A.F acknowledges financial support from a TMR grant from the European Union (contract nr. ERBFMBICT971902). K.D.K. is grateful to the Max-Planck-Institut für Gravitationsphysik (Albert-Einstein-Institut), Potsdam, for generous hospitality.

## References