text:
# Inertial modes in slowly rotating stars : an evolutionary description

## I Introduction

In 1998, Andersson anders98  discovered that in any relativistic rotating perfect fluid, r-modes111In this paper, we shall use the terminology of the hydro community : inertial modes are the oscillatory modes of a fluid whose restoring force is the Coriolis force. R-modes \[or purely toroidal (axial) modes\] then form a sub-class of inertial modes. The later were discovered by Thompson in 1880 thom80 . The reader can find a short point of history in Rieutord 2001 rieut01 . are unstable via the coupling with gravitational radiation. This is just a particular case of the Chandrasekhar-Friedman-Schutz (CFS) mechanism discovered in 1970 (see chandra70  and friedschu78 ). However, a characteristic of r-modes is the fact that this instability is expected whatever the speed of rotation of the star (see also Friedman et al. 1998 friedmor98 ). This discovery triggered an important activity in the neutron stars (NS) community since such inertial instabilities might contribute to explanation of the observed relative slow rotation rates of young and recycled NS (see Bildsten bil98 ). Moreover, such an instability could in principle make of NS efficient sources of gravitational waves (GW) for observation by the ground interferometric detectors under construction or their direct descendants. Nevertheless, as various not always well-known physical processes take place which might inhibit the growth of r-modes in real NS, their relevance is still questionable at this time. More details can be found in the reviews by Andersson and Kokkotas 2001 andko01  and Friedman and Lockitch 2001 friedlock01 .

The present study was originally motivated by the results announced by Lindblom et al. in 2001 lintv01 . These authors computed, in the highly non-linear regime, the evolutionary track of the CFS instability of r-modes in a Newtonian NS with a magnified approximate post-Newtonian radiation reaction (RR) force. In their calculation, they found a rapid growth of the mode, until strong shocks appeared and quickly damped it. A recent article by Arras et al. arr02  asserts that this phenomenon was just an artifact linked with the huge RR force. Yet, we found those results so interesting and so surprising that we decided to try to reproduce and better understand them with a different approach. The code of Lindblom et al. lintv01  was written in cylindrical coordinates with 3D finite difference scheme and they studied fast rotating stars. We chose to create a spectral code in spherical coordinates and to begin with a slowly rotating NS.

The use of spectral methods was motivated by the fact that they can provide a deep insight in describing the turbulence generated by quadratic terms. This is analogous to what is done with Fourier analysis in the study of homogeneous turbulence (e.g. Lesieur 1987 lesieur87 ). Furthermore, we chose a slowly rotating NS for several reasons that will be explained in the following.

Since NS are known for their quite rapid rotation and since inertial modes have for restoring force the Coriolis force, an a priori reasoning would probably lead to the conclusion that the only - or most - interesting case of inertial modes to study is the case of modes in a fast rotating NS. However, the final answer is not so easy to decide, as among observational data and among models of cooling of pulsars, many things in fact support slowly rotating newborn single NS. For instance, the estimations of the rotation periods at birth of the two historical X-rays and radio pulsars whose ages are known exactly (the 820 year old 65.8 ms period pulsar in SNR 3C58 MSPS02 , CS2002  and the 948 year old 33 ms period Crab pulsar) are respectively 60 ms MSPS02  and 14 ms SBJM94 . Moreover, compared with assessment of the angular momentum of isolated supergiant stars (whose core collapses will give type II supernovæ and then potentially NS), these numbers show that important losses of angular momentum are expected to happen during the stellar evolution. Furthermore, this is the same concerning the giant phase of less massive stars that give white dwarfs. Indeed, the observed rotation periods of 19 white dwarfs range between 12.1 min and 12 days with a median value of about 1 hour Scmi01  but, if the sun shrank without losing any angular momentum into a white dwarf of the same mass with a radius of 5000 km, its rotation period would be of a few minutes only. In the stellar evolution community, several ways to explain these weak angular momenta can be found and the final answer is still not clear. But a quite admitted idea is that magnetic field plays probably the key-role via magnetic braking type mechanisms (Schatzman 1962 scha62 ). And whatever the actual mechanism, the current conclusion in this community is that a huge amount of angular momentum is supposed to be lost during the stellar evolution for main sequence stars of any mass. Then, the most difficult thing to explain in stellar evolution physics seems to be the reason why these rotation periods at birth (of the orders of 1 hour for white dwarfs and of 1 ms for NS) are so small (e.g. Spruit 1998 spr98  and Spruit and Phinney 1998 spph98 ): stellar evolution scenarios do not expect of baby NS to be fast rotating.

However, it can be mentioned that some isolated fast rotating NS are found. But, for the most part of them, the current idea is that these stars have been recycled in binary systems. In these systems, where fast rotating NS exist, the NS is thought to have been spun up by the accreted matter. The recent discoveries of accreting millisecond pulsars (XTEJ0929-314, SAXJ1808.4-3685, XTEJ1751-305 and 4U 1636-53) whose rotation frequencies are respectively 185 Hz, 401 Hz, 435 Hz and 581 Hz STRO02  (corresponding to respective periods of 5.4 ms, 2.5 ms, 2.3 ms and 1.7 ms) support the above scenario. Yet, there is also an example of single ms pulsar associated to a supernova remnant: PSR J0537-6910 whose period is 16 ms, whose age is about 5000 yr and which is supposed to have had an initial period of a few ms (see for instance mar98 ). Hence, in spite all the above arguments, the existence of rapidly rotating baby NS should not be completely excluded. Moreover, a fast rotating baby NS with very weak magnetic field and consequently very weak losses of angular momentum due to magnetic braking like mechanisms during the supergiant phase would escape to the observations as pulsars but would be very interesting as GW sources.

Nevertheless, such a discussion does not make clear the important point. Indeed, whatever the value of the frequency, it does not directly tell if a pulsar is fast rotating or if it is not. What says, in a particular work, if a NS can be regarded as (quite) slowly rotating is the relative importance of the deformation for this study. In a more general framework, this question is settled by the ratio between the angular velocity of the pulsar and its Keplerian angular velocity. Even for the fastest rotating known pulsars (that are in binary systems) this ratio is less than a third, since the Kepler frequency is around 1 ms (see haen89 ). Furthermore, what is implied in the hydrodynamics of a NS is not this ratio, but its square. Hence, in appropriate units, this factor is less than ten per cent of the Coriolis force for most of the NS. Anyway, there is a last crucial issue. In a Newtonian star, even if this factor is a few per cent correction in the equation of motion, it is fundamental since it creates a coupling between polar and axial part of the velocity. Yet, in a relativistic star, the situation is completely different as noticed Kojima koj98 . Indeed, in a relativistic rotating star, the frame dragging term has the same qualitative result: it makes a coupling between the polar and the axial parts of the velocity. But, in appropriate units, this coupling is scaled by the ratio between the angular velocity and the Keplerian angular velocity and not by the square of this ratio. Hence, even for the fastest rotating known NS, the deformation introduces a kind of second order correction that can be neglected in a first approach222As an example, we verified that, for a rotation frequency of 300 Hz in a star of $`1.7M_{}`$ with a fully relativistic code for stationary configuration of rotating stars, the coupling between the spherical harmonic terms due to the drag effect is one order of magnitude larger than the coupling due to the deformation of the star..

Thus, our choice of using the slow rotation limit as a first step was motivated by all the astrophysical reasons mentioned above: even for ms pulsar, the slow rotation approximation is still quite good. But, secondly, we wanted to look for a possible saturation due to non-linear coupling that may occur before a highly non-linear regime is reached. To better understand such a phenomenon, we thought it was probably wiser to begin with an easier situation in which there are not several effects with consequences of the same orders of magnitude. Thus, we began to build a non-linear hydrodynamics code using the Newtonian theory of gravity and the slow rotation approximation. But, once a linear Newtonian code had been written (first step to a non-linear version), upgrading it to a general relativity (GR) linear code with strong Cowling approximation333In December 2001, during a workshop on r-modes which took place at the Meudon site of Paris Observatory, Carter suggested to call strong Cowling approximation the approximation in which all the coefficients of the metric are frozen. This name was chosen to contrast with what should be called the weak Cowling approximations where some perturbations of the metric are allowed. See for instance the work of Ruoff et al. ruoff01 , ruoff01b  and Lockitch et al. lockandf01  using the Kojima equations koj92 . was quite obvious and we chose it for a second approach to inertial instabilities. Indeed, as it was already mentioned above, Kojima first noticed in 1998 koj98  that the frame dragging phenomenon makes the relativistic r-modes quite different from the Newtonian r-modes. But concerning this linear relativistic study, we also decided to begin with the slow rotation limit to try to get a better understanding of the spherical relativistic case, before putting what can be seen as a second order correction.

This article is organized as follows: in Section II, we start by describing the physical conditions we chose in the Newtonian case. Then we write the Navier Stokes equations (NSE) in dimensionless form, and the RR force we introduced in the linear study. This enables us to calculate the associated characteristic numbers and to observe that the corresponding numerical problem is a stiff one with typical time scales being orders of magnitude different. We discuss some approximations that can be made in Section III. The most important of them is the anelastic approximation. Section IV explains the basic tests of the code, with for instance the linear $`l=m`$ r-modes in a rigidly rotating spherical fluid. This admits an easy and analytical solution of Euler equation (EE) in the divergence free approximation. We demonstrate that, thanks to spectral methods, the velocity profile is exactly determined and preserved within the round-off errors. We also show that the error in the conservation of the energy is only due to the time discretization and that it vanishes as the cube of the time-step $`\mathrm{\Delta }t`$. Finally, this section ends with tests of the anelastic approximation and of the RR force we adopted. Section V deals with results obtained in a differentially rotating NS. As expected, pure r-modes are replaced with inertial modes which are still CFS unstable, even if they are no longer purely axial. Moreover, we show that if the amplitude of the difference of angular velocity between the interior and the surface of the star is sufficiently large, these modes quite rapidly develop a kind of “singularity” in the first radial derivative of the velocity. There is then a kind of concentration of the motion near the surface and mainly near the equatorial plane. Our time evolution approach thus provides results of the same kind of those obtained by Karino et al. in 2001 karino01  with a modes calculation. A viscous term has to be added in order to regularize the solution, but it does not change the qualitative result. We shall discuss this result in the conclusion. In Section VI, we give some first results for the GR case in a slowly rotating NS with the strong Cowling approximation. In this framework, the above conclusions still hold. More results in GR will be given in a following article in preparation. The Section VII contains the conclusion and discussion. A detailed description of the numerical algorithm is given in the appendix.

## II Equations and numbers

### II.1 The barotropic case in Newtonian gravity

Cooling calculations (e.g. Nomoto and Tsuruta 1987 nom87 , Yakovlev et al. 2001 yak01 ) showed that several minutes are enough to enable the NS matter to fall far below its Fermi temperature, that is roughly $`10^{10}`$ K. Furthermore, in the Newtonian non-linear hydrodynamics of a not too young slowly rotating NS, it is not worth trying to use a very sophisticated EOS. We shall therefore adopt a barotropic EOS. With those assumptions, we have

$$P=P[n]\text{,}$$
(1)
where $`P`$ is the pressure and $`n`$ the mass density, and the Poisson equation for the gravitational potential $`U`$ is

$$\mathrm{\Delta }U=4\pi Gn.$$
(2)
The Newtonian Navier-Stokes equations (written in the inertial frame for a rigidly rotating NS) are

$`\left(_t+\mathrm{\Omega }_\phi \right)\stackrel{}{W}+\left(\stackrel{}{W}\stackrel{}{}\right)\stackrel{}{W}+\mathrm{\hspace{0.17em}2}\stackrel{}{\mathrm{\Omega }}\stackrel{}{W}+\stackrel{}{}H=`$ (3)
$`{\displaystyle \frac{1}{n}}\left(\stackrel{}{}\left(\stackrel{}{}\eta \stackrel{}{W}\right)+\stackrel{}{}\left(\stackrel{}{W}\stackrel{}{}\eta \right)+\eta \mathrm{\Delta }\stackrel{}{W}\stackrel{}{W}\mathrm{\Delta }\eta +\stackrel{}{}\left(\left(\zeta +{\displaystyle \frac{\eta }{3}}\right)\stackrel{}{}\stackrel{}{W}\right)\right)+\stackrel{}{A}.`$
Here we take $`(r,\vartheta ,\phi )`$ for a spherical system of coordinates in the inertial frame. In this system, $`\vartheta =0`$ coincides with the direction of the rotation axis of the NS, which is parallel to the angular velocity: $`\stackrel{}{\mathrm{\Omega }}=\mathrm{\Omega }\stackrel{}{e}_z`$. $`\stackrel{}{W}`$ is the velocity in the corotating frame, i.e. the part that is added to the velocity of the rigidly rotating background when a mode is present. Note that both velocities can be of the same order in the non-linear case. This is the reason why we shall not use the term Eulerian perturbation for $`\stackrel{}{W}`$. Otherwise, $`H=\frac{1}{n}\frac{\mathrm{d}P}{\mathrm{d}n}dn`$ is the enthalpy, $`\eta `$ and $`\zeta `$ are respectively the dynamical shear and bulk viscosity coefficients, and finally $`\stackrel{}{A}`$ contains any external accelerations (or force per unit of mass). The modifications that are needed in the case of the linear study with differential rotation or in the relativistic linear case are respectively discussed in sections V and VI.

In what follows, $`\stackrel{}{A}`$ is mainly the effective gravitational acceleration, i.e. the gradient of the difference between the centrifugal potential $`\frac{1}{2}\mathrm{\Omega }^2\rho ^2`$ and the gravitational potential $`U`$, where $`\rho `$ is the distance from the rotation axis. In the linear regime, we will sometimes introduce a RR acceleration (cf. Blanchet 1993 and 1997 blan937 ) of the form

$$A_{RR}^{}{}_{j}{}^{}\widehat{=}\frac{4}{c^2}v^i_j\mathrm{{\rm Y}}_i,$$

where

$$\mathrm{{\rm Y}}_i\widehat{=}\frac{4}{45}\frac{G}{c^5}ϵ_{ijk}x^jx^mS_{}^{km}{}_{}{}^{\left(5\right)}[t]$$
(4)
and

$$S_{km}[t]\widehat{=}ϵ_{pq(k}\mathrm{d}^3xnx_{m)}x^pv^q,$$

with $`\left(km\right)\widehat{=}\frac{1}{2}(km+mk)`$, $`v^i`$ being the full velocity and the superscript (5) the fifth time derivative, which can not be easily calculated in a numerical work (see Rezzolla et al. 1999 rezz99 ). Note that this formula is valid only if written with Cartesian components of the tensors and this is the reason why we did not really make distinction between contra and covariant components. Finally, we insist on the fact that such an acceleration does not include any mass multipole but only the current quadrupole that is the most important coefficient for the emission of GW by axial modes in a slowly rotating NS.

To fulfill the system of equations, we need to add the mass continuity equation:

$$\left(_t+\mathrm{\Omega }_\phi \right)n+\stackrel{}{W}\stackrel{}{}n+n\stackrel{}{}\stackrel{}{W}=0.$$
(5)
Taking into account that the EOS is a barotropic one, the preceding equation can be written in a slightly different way, but we will discuss it in the Section III.

As far as inertial modes are concerned, the typical time scale and length scale are $`2\pi \mathrm{\Omega }^1`$ and $`R`$, characteristic length of the background star. The velocity associated with those values, $`R\mathrm{\Omega }`$, can be very far from the characteristic velocity of the mode. Another velocity, scaling that of the mode, must be introduced. But instead, we introduce $`\alpha `$, a pure number that is defined as the ratio of these two velocities. Otherwise, the typical mass is obviously the star’s, $`M_{}1.4M_{}`$. Bearing all this in mind, we define the following dimensionless variables:

$$\xi \widehat{=}\frac{1}{R}r,\tau \widehat{=}t\mathrm{\Omega },\stackrel{}{V}\widehat{=}\frac{1}{\alpha \mathrm{\Omega }R}\stackrel{}{W},$$

$$d\left[\stackrel{}{\xi }\right]\widehat{=}\frac{1}{n_0}n\left[\stackrel{}{r}\right],s\left[\stackrel{}{\xi }\right]\widehat{=}\frac{1}{\nu n_0\mathrm{\Omega }R^2}\eta \left[\stackrel{}{r}\right],$$
(6)
$$b\left[\stackrel{}{\xi }\right]\widehat{=}\frac{1}{\lambda n_0\mathrm{\Omega }R^2}\zeta \left[\stackrel{}{r}\right],\stackrel{}{e}_z\widehat{=}\frac{1}{\mathrm{\Omega }}\stackrel{}{\mathrm{\Omega }},$$

$$\stackrel{~}{H}\left[\stackrel{}{\xi }\right]\widehat{=}\frac{1}{\alpha \mathrm{\Omega }^2R^2}H\left[\stackrel{}{r}\right]\text{and}\stackrel{}{\mathrm{\Sigma }}\left[\stackrel{}{\xi }\right]\widehat{=}\frac{1}{\alpha \mathrm{\Omega }^2R}\stackrel{}{A}\left[\stackrel{}{r}\right].$$

Thus the motion equations are written \[with same conventions than in Eq. (3)\] as

$`\left(_\tau +_\phi \right)\stackrel{}{V}+\alpha \left(\stackrel{}{V}\stackrel{}{}\right)\stackrel{}{V}+\mathrm{\hspace{0.17em}2}\stackrel{}{e}_z\stackrel{}{V}+\stackrel{}{}\stackrel{~}{H}=`$
$`{\displaystyle \frac{\nu }{d}}\left(\stackrel{}{}\left(\stackrel{}{}s\stackrel{}{V}\right)+\stackrel{}{}\left(\stackrel{}{V}\stackrel{}{}s\right)+s\mathrm{\Delta }\stackrel{}{V}\stackrel{}{V}\mathrm{\Delta }s+\stackrel{}{}\left(\left({\displaystyle \frac{\lambda }{\nu }}b+{\displaystyle \frac{s}{3}}\right)\stackrel{}{}\stackrel{}{V}\right)\right)+\stackrel{}{\mathrm{\Sigma }},`$
where the $`\stackrel{}{}`$ and $`\mathrm{\Delta }`$ operators are now performed with dimensionless variables ($`\xi ,\vartheta ,\phi `$) instead of ($`r,\vartheta ,\phi `$).

Here, we have introduced the following notations:

* $`\alpha `$ a pure number, characteristic of the initial amplitude of the mode, but also of the ratio of the non-linear term and the Coriolis force. With our conventions, $`\alpha `$ is twice the usual Rossby number and $`\stackrel{}{V}`$ is a vector whose norm at $`t=0`$ is equal to 1 in a given point on the equator;
* $`n_0=M_{}/(\frac{4}{3}\pi R^3)`$. For the full slow rotation limit approximation (only terms linear in $`\mathrm{\Omega }`$ and spherical shape), $`R`$ would be the radius of the star and then $`n_0`$ the mean density.
* $`s`$(hear) and $`b`$(ulk) dimensionless functions, and $`\nu `$ and $`\lambda `$ pure numbers, all chosen in such a way that if there is any viscosity, $`s`$ and $`b`$ are equal to $`1`$ in a specific position. Typically, for a single fluid model, this would be the center of the star. Nevertheless, it is worth pointing out that to build a more realistic model of NS, several different layers could be made to coexist. In this case, there would be different $`\nu `$ and $`\lambda `$, hugely depending on the shell (see for instance Haensel et al. 2001 haen01 ). With our definitions, those numbers, $`\nu `$ and $`\lambda `$, are twice the usual Ekman numbers, and they then quantify the ratios between the viscosities and the Coriolis force.
* $`\stackrel{~}{H}`$ and $`\stackrel{}{\mathrm{\Sigma }}`$ respectively the dimensionless enthalpy and the dimensionless external accelerations, both scaled by the inverse of $`\alpha `$.

Concerning the latter quantities, we cut them in two parts and use for variables the difference between their present values and their values in the steady state. For the background parts, we then have

$$\stackrel{}{\mathrm{\Sigma }}_0=\stackrel{}{}\stackrel{~}{U}+\stackrel{}{}\frac{\left(\xi \mathrm{sin}\vartheta \right)^2}{2\alpha }\text{with}\stackrel{}{}\stackrel{~}{H}_0=\stackrel{}{\mathrm{\Sigma }}_0.$$
(7)
This equality enables the background to be solution of the NSE and can be solved separately. From now until the end, we consider that this condition is realized, and we forget about $`\stackrel{}{\mathrm{\Sigma }}_0`$ and $`\stackrel{~}{H}_0`$ (except in the mass conservation equation where the latter still appears as an external parameter).

Finally, some words about the Newtonian gravitational field. Inertial modes are current oscillations associated with small density variations. In this context, it is natural to use the Cowling approximation (Cowling 1941 cow41 ) which consists in forgetting fluctuations of the gravitational field (see the relevance of this approximation for Newtonian r-modes in Saio 1982 saio82 ). We will apply it for the non-linear study. Nevertheless, for the Newtonian linear study, we shall see later that depending on the way mass conservation is treated, it may be not necessary (see Section III for more details). In this case, we have

$$\stackrel{}{\mathrm{\Sigma }}=\stackrel{}{\mathrm{\Sigma }}_0+\stackrel{}{\sigma }\text{with}\stackrel{}{\sigma }=\stackrel{}{}\delta \stackrel{~}{U}$$
(8)
where $`\delta \stackrel{~}{U}`$ is the Eulerian perturbation of the dimensionless gravitational potential.

### II.2 Dimensionless equations of motion for the spherical components of the velocity

The equations of motion for the spherical components of the velocity in the orthonormal basis associated with $`(\xi ,\vartheta ,\phi )`$ are given by

$`_\tau V_r+_\phi V_r+\left(\stackrel{}{V}\stackrel{}{}\right)V_r{\displaystyle \frac{1}{\xi }}\left(V_\vartheta ^2+V_\phi ^2\right)2\mathrm{sin}\vartheta V_\phi +_\xi \stackrel{~}{h}=`$
$`{\displaystyle \frac{\nu }{d}}(\left(_\xi s\right)\left(2_\xi V_r\right)+s(\mathrm{\Delta }V_r{\displaystyle \frac{2}{\xi ^2}}({\displaystyle \frac{1}{\mathrm{sin}\vartheta }}_\vartheta \left(\mathrm{sin}\vartheta V_\vartheta \right)+{\displaystyle \frac{1}{\mathrm{sin}\vartheta }}_\phi V_\phi +V_r))`$
$`_\xi s(\stackrel{}{}\stackrel{}{V})+_\xi (({\displaystyle \frac{\lambda }{\nu }}b+{\displaystyle \frac{s}{3}})\stackrel{}{}\stackrel{}{V}))+\sigma _r,`$
$`_\tau V_\vartheta +_\phi V_\vartheta +\left(\stackrel{}{V}\stackrel{}{}\right)V_\vartheta +{\displaystyle \frac{1}{\xi }}\left(V_\vartheta V_rV_\phi ^2\mathrm{cot}\vartheta \right)2\mathrm{cos}\vartheta V_\phi +{\displaystyle \frac{1}{\xi }}_\vartheta \stackrel{~}{h}=`$
$`{\displaystyle \frac{\nu }{d}}(\left(_\xi s\right)(_\xi V_\vartheta +{\displaystyle \frac{1}{\xi }}(_\vartheta V_rV_\vartheta ))+s(\mathrm{\Delta }V_\vartheta +{\displaystyle \frac{1}{\xi ^2}}(2_\vartheta V_r{\displaystyle \frac{V_\vartheta }{\mathrm{sin}^2[\vartheta ]}}{\displaystyle \frac{2\mathrm{cos}\vartheta }{\mathrm{sin}^2[\vartheta ]}}_\phi V_\phi ))`$ (9)
$`+{\displaystyle \frac{1}{\xi }}_\vartheta (({\displaystyle \frac{\lambda b}{\nu }}+{\displaystyle \frac{s}{3}})\stackrel{}{}\stackrel{}{V}))+\sigma _\vartheta ,`$
$`_\tau V_\phi +_\phi V_\phi +\left(\stackrel{}{V}\stackrel{}{}\right)V_\phi +{\displaystyle \frac{1}{\xi }}\left(V_\phi V_r+V_\phi V_\vartheta \mathrm{cot}[\vartheta ]\right)+2\left(\mathrm{sin}\vartheta V_r+\mathrm{cos}\vartheta V_\vartheta \right)+{\displaystyle \frac{1}{\xi \mathrm{sin}\vartheta }}_\phi \stackrel{~}{h}=`$
$`{\displaystyle \frac{\nu }{d}}(\left(_\xi s\right)(_\xi V_\phi +{\displaystyle \frac{1}{\xi \mathrm{sin}\vartheta }}(_\phi V_r\mathrm{sin}\vartheta V_\phi ))+s(\mathrm{\Delta }V_\phi +{\displaystyle \frac{1}{\xi ^2}}({\displaystyle \frac{2}{\mathrm{sin}\vartheta }}_\phi V_r+{\displaystyle \frac{2\mathrm{cos}\vartheta }{\mathrm{sin}^2[\vartheta ]}}_\phi V_\vartheta {\displaystyle \frac{V_\phi }{\mathrm{sin}^2[\vartheta ]}}))`$
$`+{\displaystyle \frac{1}{\xi \mathrm{sin}\vartheta }}_\phi (({\displaystyle \frac{\lambda b}{\nu }}+{\displaystyle \frac{s}{3}})\stackrel{}{}\stackrel{}{V}))+\sigma _\phi ,`$
where

$$\mathrm{\Delta }f\widehat{=}\frac{1}{\xi }_\xi ^2\left(\xi f\right)+\frac{1}{\xi ^2\mathrm{sin}\vartheta }_\vartheta \left(\mathrm{sin}\vartheta _\vartheta f\right)+\frac{1}{\xi ^2\mathrm{sin}^2\vartheta }_\phi ^2f$$

is the scalar Laplacian, while

$$\left(\stackrel{}{V}\stackrel{}{}\right)f\widehat{=}V_r_\xi f+V_\vartheta \frac{1}{\xi }_\vartheta f+\frac{1}{\xi \mathrm{sin}\vartheta }V_\phi _\phi f$$

and

$$\stackrel{}{}\stackrel{}{V}\widehat{=}\frac{1}{\xi ^2}_\xi \left(\xi ^2V_r\right)+\frac{1}{\xi \mathrm{sin}\vartheta }_\vartheta \left(\mathrm{sin}\vartheta V_\vartheta \right)+\frac{1}{\xi \mathrm{sin}\vartheta }_\phi V_\phi .$$

In the above equations, we assume $`s`$ to only depend on $`\xi `$ in the $`(\xi ,\vartheta ,\phi )`$ system of coordinates. This assumption is linked with the slow rotation limit we should apply from now. This approximation supposes that the star is slowly rotating compared with its Kepler frequency. As explained in the introduction, observational data make this assumption credible, while known pulsars seem to rotate with a velocity smaller than a third of their Kepler velocity. In the lowest framework, only terms linear in $`\mathrm{\Omega }`$ are kept, and the background star is supposed to have preserved its spherical shape. Yet, to improve the slow rotation approximation, one can go a step further and take into account the deformation of the star, assuming it is a small order effect (or small quantum number $`l`$ effect). To make this improvement, it is sufficient to introduce a new variable $`\chi [\xi ,\vartheta ]`$ that is equal to the unity at the surface, that coincides with isosurfaces of pressure, enthalpy and density, and that decomposes on $`m=0`$ Legendre functions444Restricting to the case $`m=0`$ is sufficient for an isolated NS, but in a binary system, $`m=2`$ should also be included.. With our algorithm, the easiest way to proceed would be to keep the spherical basis for the vectors, and to express all spatial operators depending on ($`\xi ,\vartheta ,\phi `$) as the sum of the new operators depending on ($`\chi ,\vartheta ,\phi `$) and of terms to interpret as applied accelerations. More details can be found in Bonazzola et al. 1997, 98, 99 bonagm97 , bonagm98  and bonagm99 . This will not be done in this article, for reasons that were also explained in the introduction. Finally, note that we do not really solve these equations. Indeed, we use the Helmholtz theorem (see for instance Morse and Feshbach 1953 mors53 ) that says that any vector of $`^3`$ can be in an unique way written as the sum of a divergence free vector and of the gradient of a potential, given its normal component over the boundary:

$$\stackrel{}{V}=\stackrel{}{}\psi +\stackrel{}{}\stackrel{}{𝒱}.$$
(10)
Instead of working with the above equations \[Eq.(II.2)\], we use the equation on the scalar potential $`\psi `$ and the equations on the $`\vartheta `$ and $`\phi `$ components of the divergence free vector $`\stackrel{}{}\stackrel{}{𝒱}`$. Most of the time, these equations cannot be reached analytically. The way to proceed is then to separate the potential and the divergence free parts of the initial equations by numerically solving Poisson like equations obtained by taking their divergence. The reader can find in the appendix more details about these algorithms.

### II.3 Characteristic numbers

To better understand what are the dominant processes in the dynamics of NS, it is worth estimating characteristic time scales and associated numbers. The acoustic time, i.e. the duration of the acoustic waves travel across the NS is by far the shortest. For a typical NS, it is about $`2.10^4`$ seconds 555This time is of the same order of magnitude as the inverse of the frequency of the fundamental pressure mode, the so-called p-mode.. As we are dealing with slowly rotating NS, the period of rotation should be more than $`2.10^3`$ seconds ($`500Hz`$). It follows that this time is at least one order of magnitude greater than the acoustic time. It is then the same for the typical period of inertial modes, their frequency being (at linear order) proportional to $`\mathrm{\Omega }`$.

Another time scale is the viscous damping time associated with viscosity:

$$T_v\frac{n_0R^2}{\mathrm{max}[\eta ,\zeta ]}$$

or

$$T_v\frac{1}{2\mathrm{max}[E_s,E_b]\mathrm{\Omega }}$$

where $`E_s`$ and $`E_b`$ are respectively the Ekman number for shear and bulk viscosity. In other words, the Ekman number must be interpreted as characteristic of the ratio between the period and the viscous time. This number is typically less than $`10^7`$ and the viscous time is more than 7 orders of magnitude greater than the period (see Cutler et al. 1987 cut87  for more precise values). A flow will be said “rotation dominated” if the above Ekman and Rossby numbers are small compared to the unity. Note that another usual hydrodynamical number appears with them: the Reynolds number, prodrome of turbulence. In a rotating fluid, it is defined by the ratio between Rossby and Ekman numbers, or in any fluid by the ratio between the non-linear and the viscous terms.

The last typical time to evaluate here is the instability rising time $`T_g`$ associated with the RR force. To get an idea of its value, we come back to the dimensionless RR acceleration, formulae (4) and definitions (6). Analytical calculation with $`V^i`$ being the $`l=m`$ linear r-mode (with time dependent amplitude)

$$\stackrel{}{V}=𝒜[\tau ]\frac{1}{m}\mathrm{}\left[\xi ^m\stackrel{}{\xi }\stackrel{}{}Y_{mm}\right]$$
(11)
or in the spherical orthonormal basis

$$\stackrel{}{V}=𝒜[\tau ]|\begin{array}{c}0\hfill \\ \xi ^m(\mathrm{sin}\vartheta )^{m1}\mathrm{sin}[m\phi ]\hfill \\ \xi ^m(\mathrm{sin}\vartheta )^{m1}\mathrm{cos}[\vartheta ]\mathrm{cos}[m\phi ]\hfill \end{array}$$
(12)
(where $`\mathrm{}`$ means the real part of the complex function), gives in dimensioned variables at the lowest order for the $`m=2`$ r-modes

$$_t𝒜[t]=\frac{GR^7\mathrm{\Omega }^6n_0}{c^7}\frac{2^{16}\pi }{3^8\mathrm{\hspace{0.17em}5}^2\mathrm{\hspace{0.17em}7}}𝒜[t]\widehat{=}\frac{1}{T_g}𝒜[t].$$
(13)
For a typical spherical NS with $`R=10`$ km, $`M=1.4M_{}`$ and $`\mathrm{\Omega }(2\pi )200`$ Hz, $`T_g`$ is something like $`10^8`$ periods of the NS. Thus, depending on the viscosity given by the EOS, it will or will not be larger than $`T_v`$, and then, the inertial instability will or will not be relevant. At this step, a new typical number seems natural to introduce. We shall propose to call “Chandra number”666Chandrasekhar was the first who studied the gravitational radiation driven instability for the $`l=m=2`$ fundamental modes of uniform density MacLaurin spheroid in 1970. See chandra70 . the ratio between the viscous time $`T_v`$ and the rising time $`T_g`$. In the same spirit as the Rossby and Ekman numbers, it should also quantify the ratio between the viscous and RR forces:

$$Ch=\frac{2^{16}\pi }{3^8\mathrm{\hspace{0.17em}5}^2\mathrm{\hspace{0.17em}7}}\frac{Gn_{0}^{}{}_{}{}^{2}R^9\mathrm{\Omega }^6}{c^7\mathrm{max}[\eta ,\zeta ]}.$$
(14)
Note finally that with this factor ahead of the physical parameters, we ensure the bifurcation value of $`Ch`$ to be of the order of the unity, at least for the linear $`l=m=2`$ r-mode. Indeed, this point is easily illustrated by looking at NSE for the linear $`l=m=2`$ r-mode with time dependent amplitude and a shear viscosity of the form $`s=1\xi ^2`$. This viscosity that vanishes at the surface implies no need to add more boundary conditions (BC) and gives the exact (at the lowest order for the RR force) differential equation for the amplitude:

$$_t𝒜[t]=\frac{1}{T_g}𝒜[t]\left(1\frac{2}{Ch}\right)$$

where one easily sees that with this shape, the viscosity wins the battle against RR force for $`Ch<2`$. All the previous numbers are gathered in Table 1.

To end with this short discussion, we should insist on the most important conclusion that is summarized in tab. 2: the above figures show how stiff is the numerical problem of finding a dynamical solution of the NSE in this framework. This is a physical situation in which several very different time scales appear. In order to get an efficient and accurate code, some approximations have to be made.

## III Mass conservation, boundary conditions and approximations

To be consistent with the variables we chose in NSE, we have to write the mass conservation equation with the dimensionless enthalpy. With the decomposition explained at the end of the Section II.1, we have

$$\stackrel{~}{H}[\tau ,\stackrel{}{\xi }]=\stackrel{~}{H}_0\left[\stackrel{}{\xi }\right]+\alpha \stackrel{~}{h}[\tau ,\stackrel{}{\xi }].$$
(15)
The first term $`\stackrel{~}{H}_0`$ is the background, the second $`\stackrel{~}{h}`$ corresponds to the mode itself and $`\alpha `$ still quantifies the non-linearity. This gives

$`\left(_\tau +_\phi \right)\stackrel{~}{h}+\left(\stackrel{}{V}\stackrel{}{}\stackrel{~}{H}_0+\mathrm{\Gamma }[n]\stackrel{~}{H}_0\stackrel{}{}\stackrel{}{V}\right)+\alpha \left(\stackrel{}{V}\stackrel{}{}\stackrel{~}{h}+\mathrm{\Gamma }[n]\stackrel{~}{h}\stackrel{}{}\stackrel{}{V}\right)=0`$ (16)
where $`\mathrm{\Gamma }[n]=\frac{\mathrm{d}\mathrm{ln}H}{\mathrm{d}\mathrm{ln}n}`$. In the polytropic case, $`P=kn^\gamma `$, $`\mathrm{\Gamma }`$ is constant and reduces to $`\mathrm{\Gamma }=\gamma 1`$. Exception done of the case $`\gamma =\mathrm{\hspace{0.17em}1}`$ where the enthalpy is a logarithm. In the linear case, we have to neglect the last part of this equation, i.e. to do $`\alpha =0`$ in Eq.(16). We shall now describe some different choices that can be made to deal with this equation. The reader can find in appendix B the algorithms to implement all the following schemes in the framework of spectral methods.

### III.1 Solving the exact system of equation

If one tried to solve numerically the exact Navier-Stokes and mass conservation equations, one would have two main possibilities. The first would be to employ an explicit scheme. But by this way, the Courant conditions for the following of the acoustic waves would impose a time step very small compared with the period of the star. This would almost forbid any hope to make evolutions during durations long from the inertial modes point of view. The second way to proceed would be to use an implicit scheme for the divergence free part (see the appendix). Indeed, this would make it possible to take a time step not too small compared with the period of the NS. But here the problem would be to estimate the errors done. Hence, in a first study, some approximations can be introduced to solve the problem in a easier way.

### III.2 The divergence free approximation

To avoid the problem of solving acoustic waves for better studying inertial modes, the easiest solution is to use the divergence free approximation. It consists in replacing the usual mass conservation equation by $`\stackrel{}{}\stackrel{}{V}=0`$. By this way the time step can be chosen larger than in the solving of the exact system and then allows the following of the inertial modes themselves. From a numerical point of view, this is a fast and robust approximation quite useful in the exploring phase of the numerical work. Yet, it should be noticed that this drastic approximation may be not too bad for inertial modes in NS. Indeed, the $`l=m=2`$ r-modes which are the most interesting for GW have a divergence free limit in the linear order. Furthermore, the latest figures about bulk viscosity (Haensel et al. 2001haen01 ) that are quite huge, could mean a fast damping of all modes, except of those that are divergence free. Finally, note that in the linear limit of EE, this approximation gives an evolution of the mode $`\stackrel{}{V}`$ that is independent of the background star if it is rigidly rotating.

### III.3 The anelastic approximation

Yet, instead of being quite “savage” with the equation and imposing the divergence free condition, one could look for a cleverer way of doing physics and think about inertial modes. A main feature of these modes is that their frequency is of the same order as the angular velocity of the star in which they occur. As their damping and growing times are larger than their period (cf. Section II.3), one would like to take it for a characteristic time scale. From practical and numerical points of view, it means a time unit (or time step) choice not too small compared with the period, in order not to waste memory and computational time calculating non relevant physics. From a physical angle, the philosophy is to neglect acoustic waves by assuming that time derivatives of the pressure and density perturbations do not play a key role in the phenomenon. It can be done in a consistent way with the anelastic approximation.

The anelastic approximation was first introduced in atmospheric physics by Batchelor in 1953 bat53  and then derived from a rigorous scale analysis by Ogura and Phillips in 1962 oguph62  who gave it its name. In astrophysics, it appeared in 1976 in a paper by Latour et al. lat76  concerning convection and was then widely used in that field and others (such as stellar oscillations) where one can neglect temporal variations of the perturbation in density but not necessarily spatial ones. For a recent critical approach of this approximation within astrophysics, one can read Dintrans and Rieutord 2001 dinrieu01  and Rieutord and Dintrans 2002 rieudin02 .

In the Eq.(16), anelastic approximation consists in neglecting $`\left(_\tau +_\phi \right)\stackrel{~}{h}`$ which is the time derivative of the enthalpy in the rotating frame. By this way, we cut acoustic waves and then have

$`\stackrel{}{V}\stackrel{}{}\stackrel{~}{H}_0+\mathrm{\Gamma }[n]\stackrel{~}{H}_0\stackrel{}{}\stackrel{}{V}+\alpha \left(\stackrel{}{V}\stackrel{}{}\stackrel{~}{h}+\mathrm{\Gamma }[n]\stackrel{~}{h}\stackrel{}{}\stackrel{}{V}\right)=0.`$ (17)
As a conclusion on this approximation, we would like to insist on a particularity of the linear case. With anelastic approximation and linearization of all equations, we obtain an equation for the mass conservation that does not depend on the Eulerian perturbation of the enthalpy. But, taking the curl of the linear EE (or NSE) equation gives an equation with exactly the same feature (remember that the curl of a gradient is zero) with an interesting additional fact: it neither depends on the Eulerian perturbation of the gravitational potential. Furthermore, it is easy to verify that for a background enthalpy depending only on $`\xi `$, the boundary condition is that the radial velocity should vanish at the surface. Thus, the situation is that finding the Eulerian perturbation of velocity should give exactly the same result, whatever the hypothesis on the Eulerian perturbation of the gravitational potential. It means, we do not need Cowling approximation with anelastic approximation in the Newtonian linear case. Cowling approximation would play a role only if we wanted to find what is the enthalpy and what is the gravitational potential in the source term of the gradient part of EE.

### III.4 Boundary conditions

The only missing information is now the boundary conditions we chose. In an actual NS, the fluid core is supposed to be surrounded by a more or less rigid crust, an ocean and an atmosphere. Their physics is by itself quite a complex subject (see for example Haensel 2001b haen01b ). But even for a simple toy-model, for instance a crust made of only one type of nuclei, depending on the temperature and on the amplitude of the motion of the inner fluid, the physical state of the crust can be quite complex to describe, something like icepack on the sea (see for instance Lindblom et al. 2000 linowus00 ). There is no need to explain how difficult would be to translate this BC in a mathematical language… This is the reason why, for the EE, we then chose to begin with two different and extreme BC to get an idea of the limit cases. The first is the free surface BC, i.e. the absence of any crust. The second is the presence of a rigid crust at $`\xi =1`$, if we take the inner radius of the crust for the typical length $`R`$. In the first case, one has $`\stackrel{~}{H}_{\xi =1}=0`$ and in the second $`V_r_{\xi =1}=0`$. The numerical ways to take into account those BC can be found in the appendix B.3. For the NSE, as it was already said in Section II.3, in order to avoid the need of more BC, we chose a degenerate shear viscosity of the form $`s=1\xi ^2`$.

## IV Test and calibration of the code

Now that we have presented our physical framework, we will focus on the code. As it was already mentioned in the previous sections, it uses spherical coordinates and spectral methods to solve NSE. More precisely, it solves the equations coming from the decomposition of NSE into a potential and a divergence free parts. We will not give more details here and send the reader to the appendix for informations about the algorithms and spectral methods. The rest of this article, and the discussion about numerical stability in the appendix, are devoted to the linear study. Non-linear work is still in progress and will be described later. Furthermore, in this section, we only deal with rigid rotation in order to have some analytical solutions to make tests with.

### IV.1 Conservation of the energy

The first test of the code was obviously the free evolution of the linear $`l=m`$ r-mode in a rigidly rotating inviscid and incompressible fluid. The vector defined in Eq.(12) is indeed an eigenvector of EE whose frequency in the inertial frame is $`w_i=\frac{(m1)(m+2)}{m+1}`$. Moreover, note that the absence of radial component implies that both free surface and spherical rigid crust BC are automatically satisfied. Yet, for all the preliminary calculations using the divergence free approximation, the code was built to work with the rigid crust BC.

The great advantage of spectral methods is to change any linear spatial operation in linear algebra calculations, which can be done exactly. As the velocity of Eq.(12) has a very small number of coefficients in the reciprocal space, it is easy to verify (for example by directly looking at the time evolution of those coefficients) that there are only errors due to round-off and to time discretization. In figure (1) is illustrated the time evolution of $`V_\vartheta `$ at the equator for the $`m=2`$ linear r-mode. The spatial lattice was of the shape $`(8\times 6\times 4)`$ for $`(N_r,N_\vartheta ,N_\phi )`$777In fact, symmetries of the spherical harmonics are used and the effective value of $`l_{max}`$ is roughly twice the number of points in $`\vartheta `$. with 100 time steps per period of the mode. The duration of this run was chosen for the evolution to last exactly 100 times the expected period. This calculation was done on a DEC Alpha Station with 500 MHz processor and took 217 seconds of CPU time. Comparison of the amplitude at the first step and at the last steps shows the growth of the amplitude due to errors that come from the time discretization. This is easier to see in figure (2) where the time evolution of the error in energy appears. The different curves illustrate results obtained with the same grid and duration, but in runs with different numbers of time steps per oscillation. Here are pictured results for 100, 150 and 200 \[see also associated power spectra in figure (3)\]. The last two calculations on the same computer took respectively 269 and 342 seconds of CPU time.

In figure (4), we drew for several runs the relative error in the energy per oscillation of the mode, versus the number of time steps per oscillation. Both are in logarithmic scales. This error appears to exactly varies as $`\mathrm{\Delta }t^3`$ and then as the inverse of the cube of the number of steps per period (regression slope of $`3.01`$). This is due to the second order scheme we use to solve the NSE or EE. Indeed, when the energy is calculated, the error of order $`\mathrm{\Delta }t^2`$ done on the velocity is multiplied by a source term linear in $`\mathrm{\Delta }t`$.

As we are now only dealing with the linear code, we will not discuss the stability of the non-linear code, nor give test pictures corresponding to modes with azimuthal numbers different from $`2`$. Indeed, there is a coupling between different values of $`m`$ only in the non-linear versions of NSE and EE. The modes that are the most interesting for GW are $`m=2`$ modes and in the following, we will always choose this kind of initial data. Nevertheless, the same calibration tests for the other linear $`l=m`$ r-modes were done and gave the same power law relations between $`\frac{\delta E}{E}`$ and the number of steps per oscillation. Finally, we also verified that the error in the phase of the modes is proportional to the square of the number of steps per period. The relative errors in the phase for the $`m=2`$ mode with 100, 150 and 200 steps per oscillation were respectively $`\mathrm{1.64\hspace{0.17em}10}^3,\mathrm{\hspace{0.17em}7.33\hspace{0.17em}10}^4\text{and}\mathrm{\hspace{0.17em}4.12\hspace{0.17em}10}^4`$. It is easy to verify that the logarithms of these numbers are on a straight line with a slope very close of $`2`$.

We shall see now how we tried to improve the conservation of energy. For the free r-mode, the energy is expected to be exactly conserved. Furthermore, we know our numerical error in the velocity is of the order of $`\mathrm{\Delta }t^2`$. Remember that it comes from the second order scheme used to calculate source terms:

$$S^{j+1/2}=\frac{3S^jS^{j1}}{2}.$$
(18)
The basic idea is thus to modify this quantity, in the coefficients space, with some additional term of the order of $`\mathrm{\Delta }t^2`$, in such a way to retrieve the conservation of the energy without changing the error in the velocity. We then modified the Eq.(18) and took

$$S^{j+1/2}=\frac{(3+\epsilon )S^j(1+\epsilon )S^{j1}}{2}$$
(19)
where $`\epsilon `$ is of the order of $`\mathrm{\Delta }t^2`$. It is calculated using values of the energy at the last two instants and in such a way to “impose” on the energy to be conserved. Obviously, when there are forces such as viscosity or RR, their power has to be taken into account in the energy balance equation. In fig.(5) are the same curves as in fig.(2) but in logarithmic scales. We added the result of this “improved conservation of energy” for the run with 100 steps per oscillation. As the correction is local in time (done at almost each time step), it leads to a time independent error that should be able to remain the same as long as one could wish. This calculation took 320 seconds of CPU time.

### IV.2 Modes driven to instability in the rigidly rotating case

After trying to improve the conservation of energy for the free case, we did tests with a RR force included. For instance, we switched on the RR force in order to drive toward instability the linear $`m=2`$ r-mode. As it has already been implied in Section II.3, we did it by using a modified version of the formulae (4) (see this section). For these calculations, we compared numerical results with analytical calculations, both reached with the same approximation. We took for the background star a homogeneous NS with $`M=1.4M_{}`$, $`R=10`$ km and $`\mathrm{\Omega }\left(2\pi \right)500Hz`$. In fig.(6) appear two time evolutions of the ratio between the energy and the inertial energy for 100 oscillations of the mode. The first calculation was done without the enhanced conservation of energy and the second with this improvement. The analytical calculation gives for the ratio a final value of the order of $`1.0575`$ whose difference with the improved numerical value, $`1.0582`$, is less than $`10^3`$.

The next step was to switch on the RR force, but without the linear r-mode for initial data, even if we still assumed that the fluid was divergence free. Instead, we chose for the initial velocity a random Gaussian noise with a first moment equal to $`0`$ and a second moment equal to $`1`$. It should be noticed that in this case and others where the initial energy is not only distributed within quite low frequency modes, the above method to improve the conservation of energy is no longer appropriate and we shall not use it. Moreover, even if this method does not change the angular momentum when applied to the $`m=2`$ part of the velocity, it can if one naively applied it to the $`m=0`$ part. But, as we already mentioned, here we only deal with the linear code and mainly with $`m=2`$ velocities. And in this case, we exactly control the error done, without changing the angular momentum. We drew in figures (7) and (8) the time evolutions of the $`r`$ and $`\vartheta `$ components of the velocity and their power spectra. In fig.(9) is illustrated the time evolution of one of the two independent components of the $`S_{ij}`$ tensor, with the associated power spectrum. These calculations were done with a lattice of the shape $`(12\times 8\times 4)`$, with $`200`$ steps per oscillation and lasted $`200`$ times the period of the $`m=2`$ r-mode. The included RR force corresponds to what should exist in a NS with $`M=1.4M_{}`$, $`R=10`$ km and $`\mathrm{\Omega }(2\pi )\mathrm{\hspace{0.17em}1050}Hz`$. Thus, this calculation is not really physical (due to the use of the slow rotation approximation) but just aims at giving a qualitative idea of what should happen with physical values. As expected, there is an axial mode (no radial velocity) driven to instability with exactly the frequency ($`w=\frac{4}{3}\mathrm{\Omega }`$) and coefficients ($`l=2`$) of the r-mode. Moreover, a single look at the figure (9) shows that even when the velocity is mainly noisy, the tensor that plays the key-role in the instability is quite smooth. In this figure, we also drew a zoom corresponding to $`t\mathrm{\Omega }<50`$ and the associated power spectrum. It shows two frequencies in this component of the tensor. One of them is the unstable r-mode with $`3w/\mathrm{\Omega }=4`$ and the other (with $`3w/\mathrm{\Omega }2.70`$) disappears with a longer run, as the scale is adapted to the growing mode. Yet, even in the spectrum of the full evolution, a trace of it can still be seen. We achieved exactly the same features with others noisy initial conditions, for instance a Dirac kick into the NS888We mean an initial velocity equal to $`0`$ anywhere except in an arbitrary point.. We also did the same calculations with other spatial lattices (up to $`64\times 48\times 4`$), and this result did not change at all.

### IV.3 Test of the anelastic approximation

The last calculations we did in a rigid and Newtonian background were to test the effect of the anelastic approximation. The idea was to compare the previous results obtained with the divergence free approximation and a rigid crust BC with results coming from the same initial conditions but with the anelastic approximation and the free surface BC. For the rest of this section, the background star is a polytrope with $`\gamma =2`$.

First, we took for initial data the linear $`l=m=2`$ r-mode and let it freely evolve. Here again, the BC do not play any role. The spatial lattice was $`(8\times 6\times 4)`$ with 100 time steps per period of the mode. The time evolution of $`V_\vartheta `$ at the equator \[same as in fig.(1)\] did not show any difference with the divergence free case, even for the power spectrum. But this is quite obvious to understand when looking at both the Eq.(17) and the EE. Indeed, we see that the linear r-mode, which is divergence free and has no radial component, is still an eigenvector of this system of equations. Then taking it for initial data in the divergence free case or in the anelastic approximation should give exactly the same evolution. We verified \[see fig.(10)\] that the radial component of the velocity does not grow and neither does its divergence. This calculation took 344.6 of CPU time without any optimization and with a very basic solving of the anelastic equation.

The next step was the driving toward instability of a mode with noise for initial data. We took exactly the same conditions (duration, lattice, initial data) as in the divergence free case, exception done of the boundary conditions. For these calculations and all that follows with the anelastic approximation, we will use the free surface condition that is automatically satisfied (see the appendix for more details). The figures (11), (12) and (13) that show respectively the time evolutions of $`V_r`$, of $`V_\vartheta `$ and of a component of $`S_{ij}`$ give a kind of summary of the results:

* the radial velocity still remains noisy and does not grow;
* the $`\vartheta `$ component grows in the same way as in the divergence free case. There is a difference between the final amplitudes that comes from the fact that the coupling to the RR force depends on an integral on the whole star of a function that is proportional to the density (and then depends on the EOS). Note that the power spectra have the same secondary peaks as in the divergence free case.
* the spectrum of the $`S_{ij}`$ tensor shows that there is only one unstable mode which has again the frequency of the linear r-mode.

The conclusion is then that the anelastic approximation gives results very close of those obtained in the divergence free case. This is due to the fact that, in spite of the presence of all modes in the initial conditions, the only growing mode is the r-mode that has no radial velocity and is divergence free.

## V Differential rotation

As we already mentioned it in the introduction, Kojima first noticed (see Kojima 1998 koj98 ) that relativistic r-modes may be quite different from what they are in a Newtonian background due to frame dragging. Yet, since the main reason for this difference is the modification of the NSE and the possible appearance of what is called a continuous spectrum (see Ruoff et al. 2001 ruoff01 , ruoff01b  and Beyer et al. 1999 beyer99 ), the same may append even in the Newtonian case if the star is not rigidly but differentially rotating. And there are several reasons why a NS may not be in rigid rotation. First, the birth conditions of the NS themselves. Secondly, the non-linear coupling of modes. Then, a possible drift induced by the existence of a magnetic field. As we are here only dealing with linear hydrodynamics and tests of the code, we will not give more details about those processes (and send the reader to the following articles for more details: Spruit 1999 spr99 , Rezzolla et al. 2000 rezls00 , 2001 rezlms01  and rezlms01b , and Schenk et al. 2002 sch02 ). What we have done is just to take as given that the background star is differentially rotating with quite an arbitrary law and to look at the influence of this law on the existence of the modes. Nevertheless, instead of asking the question “Is there any r-mode left in a differentially rotating NS?” that is not well defined, we decided to try to answer to two different and more precise questions:

* is there anything growing when a RR force is applied on noise in a differentially rotating background?
* what does happen to the linear r-modes if they are chosen to be the initial data in such kind of background?

Some lights on these questions are in the following subsections, but we shall begin with some words about the modifications implied on the equations by differential rotation.

### V.1 Modifications

Assuming differential rotation of the background star involves slight modifications of what we said in the Section II. First, terms coming from the spatial derivatives of the non constant $`\stackrel{}{\mathrm{\Omega }}[\xi ,\vartheta ,\phi ]`$ must be added in the NSE. Secondly, we shall now specify what exactly means $`\mathrm{\Omega }`$ in the definition of dimensionless variables \[cf. Eq.(6)\].

Concerning $`\mathrm{\Omega }`$, we chose to take for a time unit the inverse of its value at the equator. This is uniquely determinate due to the fact we have always assumed that the rotation law is of the form $`\mathrm{\Omega }[r,\vartheta ,\phi ]=\mathrm{\Omega }[r\mathrm{sin}[\vartheta ]]`$ or of the form $`\mathrm{\Omega }[r,\vartheta ,\phi ]=\mathrm{\Omega }[r]`$. The first case corresponds to what must be this law for the background to be stable with respect to the Newtonian EE, and the second case is in a way inspired by GR even if it is not a solution of the full Newtonian EE. For more details see Section VI.

In the dimensionless NSE, the modifications induced by differential rotation are quite simple. First, $`_\phi `$ is replaced with $`\stackrel{~}{\mathrm{\Omega }}[r,\vartheta ]_\phi `$ where $`\stackrel{~}{\mathrm{\Omega }}[r,\vartheta ]`$ is the dimensionless profile of rotation. Then, we have to add new terms coming from $`\stackrel{}{\mathrm{e}}_\phi r\mathrm{sin}[\vartheta ]\stackrel{}{V}\stackrel{}{}\left(\stackrel{~}{\mathrm{\Omega }}[r,\vartheta ]\right)`$ that are in the spherical orthonormal basis

$$|\begin{array}{cc}& 0\hfill \\ & 0\hfill \\ V_rr\mathrm{sin}[\vartheta ]_r\stackrel{~}{\mathrm{\Omega }}\hfill & +V_\vartheta \mathrm{sin}[\vartheta ]_\vartheta \stackrel{~}{\mathrm{\Omega }}.\hfill \end{array}$$
(20)
### V.2 Noise with huge RR force

Once again, our goal was to stay as close as possible of the basic and well understood situation to minimize the number of unknowns. Thus, we took noisy initial data, put it in a differentially rotating background, switched on the RR force and looked to what was to happen. By this way, the question was not to look for the existence of r-modes, but simply to look for the existence of modes driven to instability by the RR force in a non rigidly rotating background.

The only difference with the basic study done in the case of rigid rotation was that we had to choose a law for $`\mathrm{\Omega }`$. The first choice was very simple and corresponded to a background stable with respect to the Newtonian EE. It was of the form

$$\stackrel{~}{\mathrm{\Omega }}[r,\vartheta ]=𝒲(1+\beta _nr^2\mathrm{sin}[\vartheta ]^2)$$
(21)
where $`\beta _n`$ was a constant depending on the run and $`𝒲`$ calculated to have $`\stackrel{~}{\mathrm{\Omega }}[1,\frac{\pi }{2}]=1`$. As explained above, we also tried a law inspired by GR:

$$\stackrel{~}{\mathrm{\Omega }}[r,\vartheta ]=𝒲\left(1+\beta _{gr}r^2\right)$$
(22)
or laws coming from the already quoted article by Karino et al. karino01 :

$$\stackrel{~}{\mathrm{\Omega }}[r,\vartheta ]=\frac{𝒲\beta _{L}^{}{}_{}{}^{2}}{r^2\mathrm{sin}[\vartheta ]^2+\beta _{L}^{}{}_{}{}^{2}}$$
(23)
or

$$\stackrel{~}{\mathrm{\Omega }}[r,\vartheta ]=\frac{𝒲\beta _v}{(r^2\mathrm{sin}[\vartheta ]^2+\beta _{v}^{}{}_{}{}^{2})^{1/2}}.$$
(24)
Since in all these laws, none of the free parameters $`\beta `$ corresponds to a physical variable, we will not give quantitative results. It will be done in the relativistic study (and then in another article) where there is a parameter that physically quantifies the way the equations are far from the Newtonian EE: the compactness of the star. Here we will only discuss in a qualitative way the results obtained in all the previous cases and give a very representative example: what happens in a $`\gamma =2`$ polytrope with anelastic approximation (free surface) and the rotation law given by the Eq.(21) with $`\beta _n=0.4`$. Exception done of these choices, the calculation was done with exactly the same conditions as in the last study with noisy initial data in the previous section.

The main difference with the case of rigid rotation is that the mode that grows is no longer purely axial. Indeed, we can see on fig.(14) that the radial velocity is now also driven to instability. Moreover, comparing power spectra in this figure and in fig.(15) and (16) shows that this is really a single mode with frequency very close of the frequency of the r-mode. By a single mode, we mean that this is exactly the same frequency for several physical quantities and several positions in the star. This is not a sufficient condition to talk about “discrete spectrum”, but it is a necessary condition. This unavoidable existence of a polar part of the velocity in a differentially rotating Newtonian star with barotropic EOS should be compared with the results achieved by Lockitch et al. lockandf01  in the relativistic framework. Indeed, in GR, the main reason for the coupling between axial and polar parts of the velocity is the frame dragging that is imitated in a Newtonian framework by differential rotation. Finally, note that the value of the frequency in the Newtonian case is depending on the way we chose to normalize $`\mathrm{\Omega }`$. Then, to summarize all our calculation, we can say that we observed that the polar counter part of the mode appears as soon as there is differential rotation, and whatever the chosen law. Yet, the more the law for $`\mathrm{\Omega }`$ is far from the rigid case, or in more pragmatic way the greater the free parameter is, the more the unstable mode has a polar counter part.

### V.3 Free evolution

The other question we asked was: “What does happen to linear r-modes if they are put into a differentially rotating background?”. The idea in taking such kind of initial data was to have something quite close of an eigenvector, assuming if there is one it should be quite similar to the linear r-mode.

Once again, we did not make quantitative calculations and postponed it to the GR case. As in the case of noisy initial data driven toward instability by a RR force, the free evolution of the linear r-mode showed the growth of a polar part of the velocity. It can be seen on figures (17) and (18) that respectively show the time evolution of the ratio between the polar energy and the total energy and the time evolution of the radial velocity at the point of coordinates $`\left(\xi =\frac{1}{2},\vartheta =\frac{\pi }{2},\phi =0\right)`$. This results comes from a calculation done with a spatial lattice of the shape $`(64,48,4)`$, the anelastic approximation (with a free surface) and a rotation law given by the Eq.(21) with $`\beta _n=0.4`$. For reasons that will be explained later, we included degenerate viscosity (cf. Section III.4) with an Ekman number $`E_s=\mathrm{\hspace{0.17em}5.\hspace{0.17em}10}^5`$ in order to regularize the solution. The evolution was done to last $`30`$ periods of the linear r-mode with $`100`$ steps per oscillation. In the figure (17), we see that, after a while, the ratio between the energy in the polar part of the mode and the total energy reaches a kind of stationary state with a coupling between different modes. The existence of this “hybrid final state” was verified during other runs with other physical conditions and is once again to compare with results achieved in GR by Lockitch et al. lockandf01 .

Concerning the existence of modes, looking simultaneously at fig.(18), (19) and (20) shows that apparently one single mode mainly appears in the reaction force (or in the component of the corresponding tensor) even if both axial and polar part of the velocity contain several modes. Yet, the spectrum of the $`S_{ij}`$ tensor is quite noisy due to the fact that the RR force is not switched on and that the run is short. We verified this feature during other calculations. The conclusion is that for small values of the $`\beta `$ parameter (these values are depending of the chosen rotation law), the main effect of differential rotation on a “free linear r-mode” is to give it a polar counter part and to widen its spectrum. But for larger values of the parameter, even if the velocity’s spectrum becomes very noisy, the $`S_{ij}`$ tensor is always less noisy. Yet, here we will not give more details about this points, the linear Newtonian study being not really interesting from the NS point of view.

To end with this section, we will mention a quite amazing feature found in all our calculations and illustrate it with an example coming from the previous free run of a linear r-mode put in the differentially rotating background. We always found that quite fast (in something like $`15`$ periods of the linear r-mode) the main part of the velocity is most of the time concentrated in a region close of the surface of the star and of the equator. This is illustrated in the figure (21) where we drew the shape of the $`\phi `$ component of the velocity versus the radius of the star for several values of the $`\vartheta `$ angle. Note that we got the same results with different BC and even when we add viscosity to regularize the solution (and this is the case in this calculation) and to be sure this is not a numerical artifact. In fig.(22) appears the shape of $`\vartheta `$ component for an evolution done with the same grid, viscosity and initial data but with $`\beta _n=0.8`$. With such a huge value of the parameter that governs the rotation law, there is a lot of different modes that appear in the velocity spectrum. But what seems to be the most important result is the very high concentration of the velocity near the surface. This is analogous to the results achieved by Karino et al. in 2001 karino01 . In the conclusion, we will shortly discuss what is, in our point of view, a possible important repercussion of this result.

## VI Strong Cowling approximation in GR

NS are among the most relativistic macroscopic objects made of usual matter in the Universe. Moreover, the instability of r-mode is due to GR. It is then quite natural to try to study r-modes in the GR framework. The problem is that it is quite difficult to evolve dynamically the full GR and hydrodynamics equations taking into account GW. This is the reason why we decided to use the strong Cowling approximation at the beginning.

Here we will only give some preliminary results obtained in the GR case. A more complete and physical study will be in a following article in preparation. Our goal is just to show that the above conclusions concerning the appearance of a polar part of the velocity and the “concentration of motion” for Newtonian differentially rotating NS still hold in the framework of relativistic rigidly rotating NS.

As we are studying slowly rotating NS, which remain spherical, a very convenient and good approximation is to assume that the $`3`$-space is conformally flat (see Cook et al. 1996 cook96 ). Furthermore, as we are using the strong Cowling approximation, there is no GW even without this approximation. Then, in unit $`c=1`$, the infinitesimal space-time interval can be written

$`ds^2=(N[r]^2a[r]^2r^2\mathrm{sin}[\vartheta ]^2N^\phi [r]^2)dt^22a[r]^2r^2\mathrm{sin}[\vartheta ]^2N^\phi [r]dtd\phi +a[r]^2d\stackrel{}{r}^2`$
where, following the $`3+1`$ formulation, $`N[r]`$ is the lapse function and $`N^\phi [r]`$ the third contravariant component on the spherical coordinate basis of the shift vector, $`a[r]`$ is the conformal factor and finally $`d\stackrel{}{r}^2`$ is the infinitesimal interval in the $`3`$-space. In this equation, we use the isotropic gauge and the slow rotation limit to assume that all these unknown functions are only depending on $`r`$ (and this is the reason why, in the Newtonian study, we called “inspired by GR” the law of differential rotation in which the functions are only depending on $`r`$). In Hartle 1967 hartl67 , the system of coordinates is not the same as the one we use, but it is easy to show that the conclusions are identical since the mapping between the two systems is quite trivial. Furthermore, from the same article we know that $`N^\phi [r]`$ is of the first order in $`\mathrm{\Omega }`$.

To get the equations of motion, the relativistic equivalent of EE, we write the energy tensor of a perfect fluid

$$T^{\mu \nu }=(\rho +p)U^\mu U^\nu +pg^{\mu \nu }$$
(25)
where $`\rho `$ is the energy density, $`p`$ the pressure and $`U^\mu `$ the $`4`$-velocity of the perfect fluid. Then, we apply the conservation of energy:

$$_\mu T^{\mu \nu }=0.$$
(26)
If we define the generalized enthalpy:

$$\mathrm{d}H\widehat{=}\frac{\mathrm{d}P}{\rho +p},$$
(27)
it gives

$$U^\mu U^\nu _\nu H+^\mu H+U^\nu _\nu U^\mu =\mathrm{\hspace{0.17em}0}.$$
(28)
It is well-known that due to thermodynamics, the system of equations obtained by adding the baryonic number conservation $`_\mu \left(nU^\mu \right)`$ (where $`n`$ is the baryonic density) to those four equations is a degenerate system. Then, we chose to work with the baryonic number conservation and the projections on the $`3`$-space of Eq.(28).

Following the Newtonian case, we just have to add Eulerian perturbations to the rigid rotation and to linearize the equations with respect to both the amplitude of these perturbation and $`\mathrm{\Omega }`$. For the generalized enthalpy, the definition of the perturbation is obvious, but for the full $`4`$-velocity, we use the well-known results about rigid rotation of relativistic stars and write it as

$$U^\mu [t,r,\vartheta ,\varphi ]=\frac{1}{N[r]}|\begin{array}{cc}1+\delta U^0[t,r,\vartheta ,\varphi ]\hfill & \\ \left(\frac{N[r]}{a[r]}\right)^2\delta U^r[t,r,\vartheta ,\varphi ]\hfill & \\ \left(\frac{N[r]}{a[r]}\right)^2\frac{\delta U^\vartheta [t,r,\vartheta ,\varphi ]}{r}\hfill & \\ \mathrm{\Omega }+\left(\frac{N[r]}{a[r]}\right)^2\frac{\delta U^\phi [t,r,\vartheta ,\varphi ]}{r\mathrm{sin}[\vartheta ]}\hfill & \end{array}$$
(29)
Note that in this equation, $`\delta U^0`$ is not a dynamical variable but is determined according to the constraint that $`U^\mu `$ is a $`4`$-velocity ($`U^\mu U_\mu =1`$). Furthermore, $`\delta U^r,\delta U^\vartheta `$ and $`\delta U^\phi `$ are not contravariant components of a $`4`$-vector but convenient variables that are the components of a $`3`$-vector on the orthonormal basis associated with the spherical system of coordinates for the flat $`3`$-space. It enables us to write the motion equations in a way very similar to the Newtonian EE. Indeed, writing the $`3`$-velocity on the orthonormal basis associated with the spherical coordinates as

$$\stackrel{}{V}\widehat{=}|\begin{array}{c}\delta U^r\hfill \\ \delta U^\vartheta \hfill \\ \delta U^\phi \hfill \end{array}$$
(30)
and defining on the same basis

$$\stackrel{}{\varpi }[t,r,\vartheta ,\varphi ]\widehat{=}\{\begin{array}{cc}\left(\mathrm{\Omega }N^\phi [r]\right)\mathrm{cos}[\vartheta ]\hfill & \\ \mathrm{sin}[\vartheta ]\left(\mathrm{\Omega }N^\phi [r]+A[r]\right)\hfill & \\ 0\hfill & \end{array}$$
(31)
with $`A[r]\widehat{=}\left(\mathrm{\Omega }N^\phi [r]\right)\left(\frac{\mathrm{d}\mathrm{ln}[\frac{a[r]}{N[r]}]}{\mathrm{d}\mathrm{ln}[r]}\right)\frac{rN^\phi ^{}[r]}{2}`$ where  is the derivative versus the radial coordinate, we have

$$\left(_t+\mathrm{\Omega }_\phi \right)\stackrel{}{V}+\mathrm{\hspace{0.17em}2}\stackrel{}{\varpi }\stackrel{}{V}+\stackrel{}{}h=0.$$
(32)
This equation is very close of the Newtonian EE, the main difference being that the $`3`$-vector that appears instead of $`\stackrel{}{\mathrm{\Omega }}`$ is now depending on the coordinates (as in the case of differential rotation) but no longer parallel to the rotation axis.

Concerning the baryonic number conservation, writing it in a Newtonian like way, we have in the slow rotation limit

$$\left(_t+\mathrm{\Omega }_\phi \right)\stackrel{~}{n}+\left(\frac{N[r]}{a[r]}\right)^2\mathrm{div}\left(\stackrel{~}{\mathrm{n}}\stackrel{}{\mathrm{V}}\right)=0$$
(33)
where $`\stackrel{~}{n}`$ is defined as $`\stackrel{~}{n}\widehat{=}nN^2a`$, $`n`$ being the baryonic number density. The natural generalization of the anelastic approximation is then

$$\mathrm{div}\left(\stackrel{~}{\mathrm{n}}\stackrel{}{\mathrm{V}}\right)=0.$$
(34)
As we are using the strong Cowling approximation, the background star appears in the equations of motion only as “external” (from the point of view of the mode) data. For any relativistic calculation, what we do is to calculate the background configuration using the already existing code illustrated in Bonazzola et al. 1993 bona93  and then to use the resulting lapse, shift, conformal factor and their derivatives with respect to the radial coordinate in our equations.

Here we will only give one single example, showing we have in this relativistic case with the strong Cowling approximation results similar to those obtained in the case of the Newtonian differential rotation. We took a $`\gamma =2`$ polytrope with $`1.74`$ solar masses and a radius equal to $`12.37`$ km. The star was very slowly rotating (the ratio between its kinetic energy and its mass energy was about $`10^{16}`$) and we took for a time unit the inverse of the angular velocity. The anelastic approximation with the free surface BC was used and to regularize the solution, we added a degenerate viscosity with $`E_s=\mathrm{5.\hspace{0.17em}10}^5`$ once the final equations were written in a Newtonian way. We insist on the fact that this viscous term does not come from relativistic calculations and just aims at regularizing the solution. In figures (23), (24), (25), (26) and (27) are plotted exactly the same quantities as in the Section V. The conclusions for this preliminary relativistic calculation are the same as in the Newtonian case: a polar counter part of the velocity appears from the beginning of the evolution and most of the time, the velocity is concentrate near the equator and the surface.

## VII Conclusion

This article deals with a new hydrodynamical code based on spectral methods in spherical coordinates. This first version of the code was written to study inertial modes in slowly rotating NS, but it can easily be modified to include fast rotation. As it was explained in the introduction, the slow rotation limit is a very good approximation for a first step. The present version of the code can be used to study both the Newtonian case and the relativistic case with the so-called strong Cowling approximation. Here, we have given the main algorithms it uses and shown some tests of the linear version. The way we overcame the numerical instability is also explained.

Furthermore, in order to work with very different time scales, we have introduced an approximation to deal with mass (or baryonic number) conservation that proved itself quite robust and useful to study inertial modes: the anelastic approximation. Even if this approximation is not a necessary one and can be abandoned in further studies, it is very useful to have an idea of the main properties of the inertial modes in slowly rotating stars. The results presented in this article were obtained in the linear case using very basic models of NS, such as the divergence free case with a rigid crust ($`V_r_{r=1}=0`$) or a $`\gamma =2`$ polytrope with anelastic approximation and free surface. Deeper studies in the GR framework with more quantitative results on the effects of EOS, of stratification of the star and of BC will follow. Despite the naive aspect of our studies, some common features appeared. Indeed, whatever the approximation (divergence free or anelastic) with different BC (rigid crust for the divergence free case and free surface for the anelastic approximation), we saw the appearance of a polar part in the mode, as soon as the background is Newtonian and differentially rotating, or relativistic and rigidly rotating (in the latter case, this phenomenon is due to the frame dragging). Whether in the case of noisy initial data with a radiation reaction force acting on them, or in the case of the free evolution of an initial linear r-mode, from the very beginning the energy of the polar part is at least one per cent of the energy of the axial part. This result is compatible with the analytical work by Lockitch et al. lockandf01  who proved, in the relativistic framework, the existence of a non-vanishing polar part of the velocity in a sequence of NS with a barotropic equation of state and decreasing rotation rate. Furthermore, another interesting result that we obtained is the “concentration of the motion” near the surface that quickly appears: after less than $`15`$ periods of the linear r-mode in the calculations done up to now. Adding viscosity has shown that this is a robust and physical feature. Moreover, it seems to neither depend on the physical conditions or on the BC. This phenomenon is very similar to the results found by Karino et al. in 2001 karino01  with a eigenvectors calculation. But, if this is really a singularity of the derivative of the velocity near the surface, as it seems to be in our time evolutions, it would mean that a mode calculation with finite difference schemes can only give quantitative results completely depending on the number of points. Indeed, with this scheme, there is an intrinsic numerical viscosity that depends on the resolution. Moreover, if this phenomenon is verified in following studies in GR, the stratification of the NS or the physics of the crust could be very important to know if inertial modes are relevant for GW emission.

###### Acknowledgements.

We would like to thank E. Gourgoulhon and B. Carter for carefully reading this paper. The numerous remarks of the anonymous referee were also very useful for us in trying to improve this article. Then, L. Blanchet and J. Ruoff gave us the opportunity to have friendly and fruitful discussions. Finally, we would like to thank the computer department of the observatory for the technical assistance.
The numerical algorithm adopted to solve the NSE (or the EE) is based on the pseudo-spectral methods (PSM) widely used in hydro or MHD problems. Before explaining this algorithm in more details, we will begin with a short summary of the PSM in order to make more evident the peculiarity of the solving of vectorial equations like NSE. Then, follows a second appendix that deals explicitly with EE and aims at explaining the way we implemented the approximation done on the mass conservation equation.

## Appendix A Spirit of the pseudo-spectral methods

Our group developed algorithms and routines library (see Bonazzola et al. 1990, 1999 bonam90 ,bonagm99 )999In 1980, one of us (S.B.) started to build a library of routines based on spectral methods to solve PDE in different geometries. Today, this library contains more than 700 routines written in FORTRAN 70 and 90 languages. These routines are highly hierarchised and allow us to assemble codes in modular way. We call this library “Spectra”. A part of this library (the highest in the hierarchy) was written in $`C^{++}`$ language by J.A. Marck and E. Gourgoulhon in order to allow the use of an object oriented language. This library is called “Lorene”. The code described in this section uses the “Spectra” library and is written in FORTRAN 90 language. allowing us to solve partial differential equations (PDE) in different geometries, mainly in domains diffeomorphic to a sphere. First, with an example of scalar PDE, we will look at the singularities contained in operators expressed in spherical like coordinates \[see for instance the different components of the NSE: Eq.(II.2)\] and at our choice of spectral basis. Then, we shall discuss the further difficulties that arise in vectorial PDE and explain the way we overcome them.

### A.1 The scalar heat equation

Consider the heat equation:

$$\frac{F}{t}=\alpha \mathrm{\Delta }F+S$$
(35)
where $`\mathrm{\Delta }F=\frac{^2F}{r^2}+\frac{2}{r}\frac{F}{r}+\frac{1}{r^2}(\frac{^2F}{\vartheta ^2}+\frac{\mathrm{cos}\vartheta }{\mathrm{sin}\vartheta }\frac{F}{\vartheta }+\frac{1}{\mathrm{sin}\vartheta ^2}\frac{^2F}{\phi ^2})`$ is the Laplacian in spherical coordinates, $`\alpha `$ the heat conductivity supposed to be constant, and $`S`$ a source term (that may include non-linear terms if they exist). The idea of spectral methods is to look for the solution of the Eq.(35) on the form

$$F[r,\vartheta ,\phi ,t]=\underset{n,l,m=0}{\overset{\mathrm{}}{}}F_{nlm}[t]H_n^l[r]K_l^m[\vartheta ]L_m[\phi ]$$
(36)
with ($`0r1;0\vartheta \pi ;0\phi 2\pi `$) and where $`H_n^l[r],K_l^m[\vartheta ],L_m[\phi ]`$ are a well chosen complete set of functions. The problem is then to find the time evolution of the coefficients $`F_{nlm}[t]`$.

In a spherical geometry, it is quite natural to choose $`L_m[\phi ]=\mathrm{exp}^{im\phi }`$ and $`K_l^m[\vartheta ]=P_l^m[\vartheta ]`$, where $`P_l^m[\vartheta ]`$ are the Legendre functions. With these choices, the Eq.(35) can be written

$`{\displaystyle \frac{F_{lm}[r,t]}{t}}\alpha \left({\displaystyle \frac{^2F_{lm}[r,t]}{r^2}}+{\displaystyle \frac{2}{r}}{\displaystyle \frac{F_{lm}[r,t]}{r}}{\displaystyle \frac{l(l+1)}{r^2}}F_{lm}[r,t]\right)=S_{lm}[r,t]`$ (37)
where $`S_{lm}[r,t]`$ are the Fourier-Legendre coefficients of the function $`S`$ at the radius $`r`$ and the instant $`t`$.

In order to handle the singularity at $`r=0`$, we shall consider separately the cases $`l=0`$ , $`l=1`$ and $`l>1`$:

* For $`l=0`$, we use the fact that a $`C^1`$ function symmetric with respect to the inversion $`rr`$ has its first derivative that vanishes at least as $`r`$ at $`r=0`$. Therefore, for such a function, the term $`\frac{^2F}{r^2}+\frac{2}{r}\frac{F}{r}`$ is regular. Even Chebyshev polynomials $`T_{2n}[r]`$ have this property, and the choice $`H_n^0[r]=T_{2n}[r]`$ ($`n0`$) then satisfies the regularity conditions.
* For $`l=1`$, it is almost the same, but the final choice is $`H_n^1[r]=T_{2n+1}[r]`$.
* The case $`l>1`$ is more delicate to handle. Indeed, it can be shown (See Bonazzola et al. 1990 bonam90 ) that for $`F[r,\vartheta ,\phi ,t]`$ to be a $`C^{\mathrm{}}`$ function, the coefficients $`F_{lm}`$ must vanish as $`r^l`$ at $`r=0`$. It means that the $`F_{lm}[r,t]`$ are symmetric with respect to the inversion $`rr`$ if $`l`$ is even, and anti-symmetric in the opposite case. We shall then distinguish the two cases $`l`$ even and $`l`$ odd.

Case $`l>1`$ even
The functions $`H_n^l[r]=T_{2n}[r]T_{2n+2}[r]`$ are even and vanish as $`r^2`$ at the origin. Therefore the quantity $`(\frac{\mathrm{d}^2}{\mathrm{dr}^2}+\frac{2}{r}\frac{\mathrm{d}}{\mathrm{dr}}\frac{l(l+1)}{r^2})H_n^l[r]`$ is regular at the origin and is retained.

Case $`l>1`$ odd
The functions $`H_n^l[r]=(2n+1)T_{2n+1}[r]+(2n1)T_{2n+3}[r]`$ are anti-symmetric with respect to the inversion $`rr`$ and vanish as $`r^3`$ at $`r=0`$. Therefore the quantity $`(\frac{\mathrm{d}^2}{\mathrm{d}r^2}+\frac{2}{r}\frac{\mathrm{d}}{\mathrm{d}r}\frac{l(l+1)}{r^2})H_n^l[r]`$ is also regular and completes our basis.

Note that by this way, we expand the solution with a set of functions that satisfy minimal conditions of regularity, excepted for $`l=0`$ and $`l=1`$. It means that they vanish as $`r^2`$ or $`r^3`$, in order to make regular any term in the equation, instead of vanishing as $`r^l`$ as they should to form a $`C^{\mathrm{}}`$ function. Until now, in all the problems that our group treated, these minimal regularity conditions were sufficient (see Bonazzola et al. 1999 bonagm99  for more details). But we shall see later that they make the Euler equations unstable.

To end with the scalar case, we shall show how the Eq.(37) can be written when the time discretization is performed and a $`2^{nd}`$ order implicit scheme used101010Due to the use of Chebyshev polynomials, the Courant conditions for the explicit problem formulation are quite severe. In the present case $`\mathrm{\Delta }t_{max}`$ is proportional to $`1/N_{max}^4`$ (see Gottlieb et al. 1977 gott77 ), therefore implicit or semi-implicit formulation of the problem is required.:

$`F_{nlm}^{j+1}\alpha {\displaystyle \frac{\mathrm{\Delta }t}{2}}{\displaystyle \underset{p=0}{\overset{\mathrm{}}{}}}A_{pn}^lF_{plm}^{j+1}=F_{nlm}^j+\mathrm{\Delta }t\left(\alpha {\displaystyle \frac{1}{2}}{\displaystyle \underset{p=0}{\overset{\mathrm{}}{}}}A_{pn}^lF_{plm}^j+S_{nlm}^{j+1/2}\right).`$ (38)
Here $`A_{jn}^l=H_j^lO^lH_n^l`$ is the matrix of the operator $`O^l=\frac{\mathrm{d}^2}{\mathrm{d}r^2}+\frac{2}{r}\frac{\mathrm{d}}{\mathrm{d}r}\frac{l(l+1)}{r^2}`$. The index $`j+1/2`$ means a quantity $`S`$ computed at the time $`t^j+\mathrm{\Delta }t/2`$, obtained by extrapolating this quantity using its (known) values at time $`t^{j1}`$ and $`t^j`$: $`S^{j+1/2}=\left(3S^jS^{j1}\right)/2`$. The left hand side operator can be easily reduced to a pentadiagonal one, and then inverted with a number of arithmetic operations that is proportional to the maximum value $`N_{max}`$ of $`N`$.

In the case of a space dependent heat conductivity $`\alpha [r,\vartheta ,\phi ]`$, we can use a semi-implicit method (see Gottlieb et al. 1977gott77 ). It consists in solving the equation

$`F_{nlm}^{j+1}(a+br^2){\displaystyle \frac{\mathrm{\Delta }t}{2}}{\displaystyle \underset{p=0}{\overset{\mathrm{}}{}}}A_{pn}^lF_{plm}^{j+1}=`$ (39)
$`F_{nlm}^j+{\displaystyle \frac{\mathrm{\Delta }t}{2}}(a+br^2){\displaystyle \underset{p=0}{\overset{\mathrm{}}{}}}A_{pn}^lF_{plm}^j+\mathrm{\Delta }t\left(S_{nlm}^{j+1/2}+\left(\alpha abr^2\right)\left({\displaystyle \underset{p=0}{\overset{\mathrm{}}{}}}A_{pn}^lF_{plm}\right)^{j+1/2}\right)`$
where $`a`$ and $`b`$ are two coefficients that satisfy the condition $`a+br^2\alpha `$ everywhere and that we choose in such a way that $`a+br^2`$ takes the same values as $`\alpha `$ at the boundaries.

The new matrix in the L.H.S. is again a pentadiagonal matrix. Boundary conditions are imposed by adding a homogeneous solution of the Eq.(38) or of the Eq.(39)111111The case $`\alpha [1,\vartheta ,\phi ]=0`$ is called degenerate. In this case no BC are allowed, and $`a+b`$ must be $`=0`$.. The reader can find in the quoted literature more details on the spectral methods, and on how to implement boundary conditions.

To conclude, note that in the present example, we expand the solution in spherical harmonics. But it turns out that in general it is more convenient to use linear combinations of Chebyshev polynomials to treat the expansion in $`\vartheta `$. The philosophy of the spectral methods then consists in performing simple operations such as computing derivatives, primitives, integrals, multiplications or divisions by $`r`$, $`r^2`$, $`\mathrm{sin}\vartheta `$ or $`\mathrm{cos}\vartheta `$ in the coefficients space, and by making multiplications of functions in the configuration space. If necessary, the passage to a Legendre representation is performed with a matrix multiplication.

### A.2 Analysis in vectorial case

#### A.2.1 Stability of numerical vectorial PDE

We have seen how to handle coordinates singularities appearing in scalar PDE when spherical coordinates are used. But when treating vectorial PDE, as EE or NSE, the singularities are more malicious: a single look at the components of NSE \[Eq.(II.2)\] should be sufficient to convince the reader. The singular terms cannot be handled only by laying down analytical properties, as it was done in the above scalar equation. Indeed, the dependence existing among the different spherical components must be taken into account in order to make singular terms compensate each others. A simple example will illustrate this crucial point.

Consider the following constant divergence free vector $`\stackrel{}{V}`$ of Cartesian components: $`V_x=0,V_y=0,V_z=1`$. Its spherical components are $`V_r=\mathrm{cos}\vartheta ;V_\vartheta =\mathrm{sin}\vartheta ;V_\phi =0`$. Then, we have $`\mathrm{div}\stackrel{}{V}=\frac{V_r}{r}+\frac{2}{r}V_r+\frac{1}{r}(\frac{V_\vartheta }{\vartheta }+\frac{\mathrm{cos}\vartheta }{\mathrm{sin}\vartheta }V_\vartheta )=0`$. A small error (for example a round-off error) in the computation of components will obviously forbid an exact compensation between the two singular terms $`2V_r/r`$ and $`\frac{1}{r}(_\vartheta V_\vartheta +\frac{\mathrm{cos}\vartheta }{\mathrm{sin}\vartheta }V_\vartheta )`$. The consequence is the creation of high order coefficients and possible instabilities appearing in the iterative process.

Indeed, we found that the hydro-code for solving linearized EE was unstable. As expected, high frequency terms were exploding after few hundred time-steps (few periods of the r-mode) either in the case of the anelastic approximation or in the case of the incompressible approximation ($`\mathrm{div}\stackrel{}{W}=0`$). As explained above, the main reason of this instability was the non exact compensation of the singular terms contained in the different source terms in the Poisson equations (56) and (58) (see below).

To overcome this problem, we define two angular potentials $`T_h`$ and $`P_o`$ in a such a way that

$`W_\vartheta =_\vartheta P_o{\displaystyle \frac{1}{\mathrm{sin}\vartheta }}_\phi T_h`$ (40)
$`W_\phi ={\displaystyle \frac{1}{\mathrm{sin}\vartheta }}_\phi P_o+_\vartheta T_h.`$
If $`rP_o`$ and $`T_h`$ are any arbitrary set of regular functions, we are now sure that the corresponding components $`W_\vartheta `$ and $`W_\phi `$ will make the divergent terms appearing in the vectorial PDE compensate each other. Moreover, the system of equations (40) can be easily inverted. First, taking the angular divergence of $`W_\vartheta `$ and $`W\phi `$: $`\mathrm{div}_{\vartheta \phi }\stackrel{}{W}\widehat{=}(_\vartheta +\mathrm{cot}\vartheta )W_\vartheta +\frac{1}{\mathrm{sin}\vartheta }_\phi W_\phi `$ gives $`\mathrm{\Delta }_{\vartheta \phi }P_o=\mathrm{div}_{\vartheta \phi }\stackrel{}{W}`$ where $`\mathrm{\Delta }_{\vartheta \phi }\widehat{=}_\vartheta ^2+\mathrm{cot}\vartheta _\vartheta +\frac{1}{\mathrm{sin}^2\vartheta }_\phi ^2`$. If $`\mathrm{div}_{\vartheta \phi }\stackrel{}{W}`$ is then expanded in spherical harmonics, we immediately have

$$P_{o,lm}=\frac{1}{l(l+1)}(\mathrm{div}_{\vartheta \phi }\stackrel{}{W})_{lm}.$$
(41)
It is the same to compute $`T_h`$. Taking the angular curl of $`\stackrel{}{W}`$, we obtain

$$\mathrm{\Delta }_{\vartheta \phi }T_h=\mathrm{curl}_{\vartheta \phi }\stackrel{}{\mathrm{W}}$$
(42)
where $`\mathrm{curl}_{\vartheta \phi }\stackrel{}{\mathrm{W}}\widehat{=}\frac{1}{\mathrm{sin}\vartheta }(_\phi \mathrm{W}_\vartheta _\vartheta (\mathrm{sin}\vartheta \mathrm{W}_\phi ))`$.

To complete this procedure, we have to guarantee that $`W_r`$ can also compensate the singular terms generated by the components $`W_\vartheta `$ and $`W_\phi `$. A satisfactory way to proceed is to take as a new variable the quantity $`\mathrm{div}\stackrel{}{W}`$. Indeed, once the quantity $`g[r,\vartheta ,\phi ]=\mathrm{div}\stackrel{}{W}`$ is known, it turns to be easy to find $`W_r`$. Bearing in mind that $`\mathrm{div}\stackrel{}{W}=_rW_r+2W_r/r+1/r^2\mathrm{div}_{\vartheta \phi }\stackrel{}{W}`$, we have, after an expansion in spherical harmonics121212It is important to take into account the following properties of the spherical harmonics decomposition of the functions $`P_o`$ and $`T_h`$: for a given $`l`$ the product by $`r`$ of the component $`W_{r,lm}`$ of $`\stackrel{}{W}`$, a regular vectorial function, vanishes as $`r^{l1}`$. The corresponding poloidal part $`P_{o,lm}`$ behaves in the same way. The associated toroidal part is $`T_{o,(l1)m}`$ and is therefore a regular function.

$$W_{r,lm}=\frac{1}{r^2}_0^ru\left(ug[u,\vartheta ,\phi ]_{lm}+l(l+1)P_{o,lm}[u,\vartheta ,\phi ]\right)du.$$
(43)
Thus, at each time step, the strategy consists in calculating the potentials $`P_o`$ and $`T_h`$, and then in calculating back the components $`W_\vartheta `$ and $`W_\phi `$. The component $`W_r`$ is also computed at each time step by using the Eq.(43).

Finally, note that with this procedure, the potentials $`P_o`$, $`T_h`$ and the component $`W_r`$ have the correct analytical behaviour on the axis $`e_z`$, i.e they vanish as $`\mathrm{sin}^m\vartheta `$. Thus, the instability of the EE hydro-code was drastically reduced, but not completely suppressed. It still kept on growing, but about ten times slower than before. The reason of this residual instability was that, at each time step, the round-off errors did not have the correct analytical properties at $`r=0`$. Indeed, they did not vanish as $`r^l`$. Consequently high spatial frequency terms were still generated that did not have the good analytical properties and a runaway toward the instability was once again generated. To avoid this phenomenon, we project, at each time step, all the scalar quantities $`(\mathrm{div}\stackrel{}{W})_{lm}[r]`$, $`P_o,_{lm}[r]`$ and $`T_h,_{lm}[r]`$ on a Legendre space $`P_n^l[r]`$ in such a way to satisfy the analytical conditions. Note that the necessity of this procedure is due to the fact that in the EE, solved with spectral methods, no dissipative term (numerical or physical) is present. Consequently the numerical instability are not dumped.

#### A.2.2 Solving the vectorial heat equation

Now, more details concerning the algorithm to solve vectorial PDE with spectral methods. In solving NSE, we need to solve a vectorial heat equation that can be reduced to an equation of the type

$$\frac{\stackrel{}{\widehat{B}}}{t}+(\mu \stackrel{}{\widehat{B}})=\stackrel{}{\widehat{J}}$$
(44)
where $`\stackrel{}{\widehat{J}}`$ is the divergence free term of the source. We remember that it is obtained by the decomposition of the source term $`\stackrel{}{J}`$ in a potential part $`\stackrel{}{}\varphi `$ and a divergence free part $`\stackrel{}{\widehat{J}}`$: $`\stackrel{}{J}=\stackrel{}{\widehat{J}}+\stackrel{}{}\varphi `$.

The Eq.(44) can be written in the following form:

$$\frac{\stackrel{}{B}}{t}=\mu \mathrm{}\stackrel{}{B}+(\stackrel{}{B})\stackrel{}{}\mu +\stackrel{}{\widehat{J}}.$$
(45)
From a numerical point of view, the term $`\stackrel{}{S}=(\stackrel{}{\widehat{B}})\stackrel{}{}\mu `$ is considered as a source and is computed by a second order scheme using its values at the times $`t_{j1}`$ and $`t_{j2}`$.

Therefore, we have to solve the equation

$$\frac{\stackrel{}{B}}{t}=\mu \mathrm{\Delta }\stackrel{}{B}+\stackrel{}{\widehat{J}}+\stackrel{}{S}$$
(46)
with the condition $`\stackrel{}{}\stackrel{}{\widehat{B}}=0`$.

The technique generally used was already described in Bonazzola et al. 1999 bonagm99 . Nevertheless, we adopted here a slightly different approach which gives more accurate results.

To solve the Eq. (46) it is convenient to introduce the poloidal and toroidal potentials $`P_o`$ and $`T_h`$ as defined in appendix A.2.1. Here we will only consider the case of $`\mu =1`$, the semi-implicit scheme which must be used in the more general case being straightforward. With a second order scheme, we have in the harmonic representation

$`P_{o,lm}^{j+1}{\displaystyle \frac{1}{2}}\mathrm{\Delta }t\left({\displaystyle \frac{\mathrm{d}^2P_{o,lm}^{j+1}}{\mathrm{d}r^2}}+{\displaystyle \frac{2}{r}}{\displaystyle \frac{\mathrm{d}P_{o,lm}^j}{\mathrm{d}r}}{\displaystyle \frac{l(l+1)}{r^2}}P_{o,lm}^{j+1}{\displaystyle \frac{2}{r^2}}\widehat{B_r}_{lm}^{j+1}\right)=`$ (47)
$`P_{0,lm}^j+\mathrm{\Delta }t\left({\displaystyle \frac{1}{2}}\left({\displaystyle \frac{\mathrm{d}^2P_{o,lm}^j}{\mathrm{d}r^2}}+{\displaystyle \frac{2}{r}}{\displaystyle \frac{\mathrm{d}P_{o,lm}^j}{\mathrm{d}r}}{\displaystyle \frac{l(l+1)}{r^2}}P_{o,lm}^j{\displaystyle \frac{2}{r^2}}\widehat{B_r}_{lm}^j\right)+\mathrm{\Pi }_{lm}^{j+1/2}\right)`$
where $`\mathrm{\Pi }`$ is the poloidal component of $`\stackrel{}{\widehat{J}}+\stackrel{}{S}`$. Note the presence of a singular term $`\frac{2}{r^2}B_r^{j+1}`$ in this equation. Bearing in mind that $`\mathrm{div}\stackrel{}{\widehat{B}}=0`$ and that $`P_{o,lm}^{j+1}`$ vanishes as $`r^{l1}`$, it is obvious that there should be a compensation between these two terms. Nevertheless, we shall slightly modify the Eq.(47) in order to obtain more easily this compensation:

$`P_{o,lm}^{j+1}{\displaystyle \frac{1}{2}}\mathrm{\Delta }t\left({\displaystyle \frac{\mathrm{d}^2P_{o,lm}^{j+1}}{\mathrm{d}r^2}}+{\displaystyle \frac{2}{r}}{\displaystyle \frac{\mathrm{d}P_{o,lm}^j}{\mathrm{d}r}}{\displaystyle \frac{l(l1)}{r^2}}P_{o,lm}^{j+1}\right)=`$ (48)
$`P_{0,lm}^j+\mathrm{\Delta }t\left({\displaystyle \frac{1}{2}}\left({\displaystyle \frac{\mathrm{d}^2P_{o,lm}^j}{\mathrm{d}r^2}}+{\displaystyle \frac{2}{r}}{\displaystyle \frac{\mathrm{d}P_{o,lm}^j}{\mathrm{d}r}}{\displaystyle \frac{l(l1)}{r^2}}P_{o,lm}^j\right){\displaystyle \frac{1}{r^2}}\left(\widehat{B_r}_{lm}^j+lP_{0,lm}^j+\widehat{B_r}_{lm}^{j+1}+lP_{0,lm}^{j+1}\right)+\mathrm{\Pi }_{lm}^{j+1/2}\right).`$
Once $`P_{o,lm}^{j+1}`$ known, $`\widehat{B}_{r,lm}^{j+1}`$ is obtained as explained above \[see Eq.(43), with $`g=0`$\].

Finally, the equation for the toroidal part of $`\stackrel{}{\widehat{B}}`$ reduces to an ordinary scalar heat equation:

$`T_{h,(l1)m}^{j+1}+{\displaystyle \frac{1}{2}}\mathrm{\Delta }t\left({\displaystyle \frac{\mathrm{d}^2T_{h,(l1)m}^{j+1}}{\mathrm{d}r^2}}+{\displaystyle \frac{2}{r}}{\displaystyle \frac{\mathrm{d}T_{h,(l1)m}^j}{\mathrm{d}r}}{\displaystyle \frac{l(l1)}{r^2}}T_{h,(l1)m}^{j+1}\right)=`$ (49)
$`T_{h,(l1)m}^j+\mathrm{\Delta }t\left({\displaystyle \frac{1}{2}}\left({\displaystyle \frac{\mathrm{d}^2T_{h,lm}^j}{\mathrm{d}r^2}}+{\displaystyle \frac{2}{r}}{\displaystyle \frac{\mathrm{d}T_{h,lm}^j}{\mathrm{d}r}}{\displaystyle \frac{l(l1)}{r^2}}T_{h,lm}^j\right)+\vartheta _{(l1)m}^{j+1/2}\right)`$
where $`\vartheta `$ is the toroidal component of $`\stackrel{}{\widehat{J}}+\stackrel{}{S}`$.

## Appendix B Solving Euler or Navier-Stokes equations

In this section, we shall show how to implement the different schemes proposed in the Section III to handle the different time scales appearing in the EE for a slowly rotating star. In all cases, for simplicity, we will only consider the linearized equations. First, we will begin with the case where the exact system of equations is to solve. Then we shall consider the implementation of physical approximations and finally the problem of BC.

### B.1 Exact Euler equations

As usual, we start by doing a decomposition of the velocity vector $`\stackrel{}{W}`$ in a divergence free component $`\stackrel{}{\widehat{W}}`$ and a potential one:

$$\stackrel{}{W}=\stackrel{}{\widehat{W}}+\stackrel{}{\mathrm{\Psi }},\mathrm{div}\stackrel{}{\widehat{W}}=0.$$
(50)
In a second order time scheme, the divergence free component of the EE is then written as

$$\stackrel{}{\widehat{W}}^{j+1}=\stackrel{}{\widehat{W}}^j+\mathrm{\Delta }t[2\stackrel{}{\mathrm{\Omega }}\stackrel{}{W}+\stackrel{}{F}]_{DF}^{j+1/2}$$
(51)
while the potential part is

$$\mathrm{\Psi }^{j+1}=\mathrm{\Psi }^j\mathrm{\Delta }t[\phi _\mathrm{\Omega }+h+\phi ]^{j+1/2}.$$
(52)
In the Eq.(51), the subscript $`DF`$ means the divergence free component of the external force. In the Eq.(52), the subscript $`\mathrm{\Omega }`$ is for the potential component of the Coriolis force $`2\stackrel{}{\mathrm{\Omega }}\stackrel{}{W}`$, $`\phi `$ is the potential component of some other exterior force and $`h`$ is the variation of the enthalpy. Moreover, we also have the baryonic number (or mass) conservation equation:

$$h^{j+1}=h^j\mathrm{\Delta }t\left(\left(\widehat{W}_r+\frac{\mathrm{\Psi }}{r}\right)\frac{H_0}{r}+\mathrm{\Gamma }H_0\mathrm{\Delta }\mathrm{\Psi }\right)^{j+1/2}$$
(53)
where $`H_0`$ is the enthalpy of the non perturbed configuration. In this explicit second order scheme, the required typical value of the time step $`\mathrm{\Delta }t`$, which guarantees the numerical stability, is determined by the term $`_rH_0`$ that is much larger than the Coriolis term for a slowly rotating NS.

The implicit version of the system of equations (52), (53) is the following:

$`\mathrm{\Psi }^{j+1}=\mathrm{\Psi }^j\mathrm{\Delta }t\left({\displaystyle \frac{1}{2}}\left(h^{j+1}+h^j\right)+[\phi _\mathrm{\Omega }+h+\phi ]^{j+1/2}\right)`$ (54)
$`h^{j+1}=h^j{\displaystyle \frac{1}{2}}\mathrm{\Delta }t\left(\left(_r\mathrm{\Psi }^{j+1}+_r\mathrm{\Psi }^j+\widehat{W}_r^{j+1}+\widehat{W}_r^j\right)_rH_0+\mathrm{\Gamma }H_0\left(\mathrm{\Delta }\mathrm{\Psi }^{j+1}+\mathrm{\Delta }\mathrm{\Psi }^j\right)\right)`$
where $`\widehat{W}_r^{j+1}`$ is obtained from the Eq.(51). Solving the above system of equations reduces to invert, for each value of $`l`$ and $`m`$ an enneadiagonal ($`2N_r\times 9`$) matrix (9 diagonals) where $`N_r`$ is the number of degrees of freedom in $`r`$. The generalization to the NSE is quite obvious and we shall not give more details about it.

### B.2 Implementation of physical approximations

As it was already explained in the Section III, for solving the EE or NSE, we chose to use approximations in order to better control the results.

#### B.2.1 The divergence free approximation

The spirit of this approximation consists in replacing the mass conservation equation (53) by the condition

$$\mathrm{div}\stackrel{}{W}=0$$
(55)
or equivalently $`\mathrm{\Psi }=0`$. The problem is then to find the enthalpy $`h`$ in a such a way that the condition (55) is satisfied. This can be done easily by solving the Poisson equation

$$\mathrm{\Delta }h\mathrm{div}F^{j+1/2}=0$$
(56)
reached by taking the divergence of the EE. As the gradient of an harmonic function is a divergence free vector, such a function can then be added to a particular solution of the Eq.(56) in order to satisfy the boundary conditions. Once $`h`$ is obtained, the EE has to be solved with $`\stackrel{}{\widehat{F}}=\stackrel{}{F}\stackrel{}{}h`$. If viscous terms are present, NSE case, a vectorial type equation must be solved.

#### B.2.2 Anelastic approximation

As it was already explained, the anelastic approximation consists in neglecting the term $`_th`$ in the Eq.(53). Once again, the game consists in finding the enthalpy $`h`$ in a such a way that the following equation is satisfied:

$$\mathrm{\Gamma }H_0\mathrm{div}\stackrel{}{W}+W_r_rH_0=0.$$
(57)
To find $`h`$, it is sufficient to derive with respect to time the Eq.(57), and to replace $`_t\mathrm{div}\stackrel{}{\mathrm{W}}`$ by its value reached from the EE. It gives

$$\mathrm{\Gamma }H_0\mathrm{\Delta }h+_rH_0_rh=\mathrm{\Gamma }H_0\mathrm{div}\stackrel{}{F}^{j+1/2}+F_r^{j+1/2}_rH_0$$
(58)
where $`\stackrel{}{F}`$ contains all the force terms (Coriolis force included).

The solution of the above equation is achieved with a semi-implicit scheme very similar to the one used to solve the Eq.(35). The problem is then reduced to solve the following Equation (for simplicity, we shall consider only the case of a polytropic EOS):

$$(\gamma 1)(a+br^2)\mathrm{\Delta }h+2br_rh=\overline{H_0}\mathrm{\Delta }h+_r\overline{H_0}_rh+\mathrm{div}\stackrel{}{\mathrm{F}}^{\mathrm{j}+1/2}+_\mathrm{r}\mathrm{H}_0\mathrm{F}_\mathrm{r}^{\mathrm{j}+1/2}$$
(59)
where $`a`$ and $`b`$ are two constants defined as was done in solving the Eq.(35). The L.H.S. operator can be easily inverted and the Eq.(59) can be solved by iteration (see Gourgoulhon et al. 2001 gougc01 ). The iteration at each time step being quite time consuming, a more convenient strategy (although less accurate) consists in replacing $`h`$ in the R.H.S. of the Eq.(59) by $`h^{j+1/2}`$.

### B.3 Boundary conditions

It was already explained how to impose BC for divergence free case. For the system of equations (54) or the Eq.(58), it is worth distinguishing two cases: either the enthalpy vanishes or does not vanish at the boundary of the integration domain. Here we will only consider the case of EE. Indeed, for NSE we chose to take a viscosity that vanishes at the surface to avoid the need of further BC.

If $`H_0|_{r=1}>0`$, the solution is quite easy : the Eq.(58) or the system of equations (54) admit an homogeneous solution that can be used to satisfy one BC.

If $`H_0(1)`$ vanishes, the system of equations (54) or the Eq (58) are degenerate and therefore no BC can be imposed. But we shall show that, in this case, the correct BC are automatically satisfied.

Indeed, consider first the case of the linearized exact system of equations (54). On the surface of the non perturbed star, we have $`_th+\stackrel{}{W}\stackrel{}{}H_0`$ that is the the correct BC. The surface $`H_0+h=0`$ then defines the profile of the perturbed star. To show that the solution in the anelastic approximation also satisfies the BC, it is more convenient to first examine the non-linear case. Here, $`h`$ is not infinitesimal, and the Eq.(57) must be replaced with the equation

$$\mathrm{\Gamma }(H_0+h)\mathrm{div}\stackrel{}{W}+\stackrel{}{W}\stackrel{}{}(H_0+h)=0$$
(60)
that is satisfied if $`h`$ is a solution of the non linearized Eq (58):

$$\mathrm{\Gamma }(H_0+h)\mathrm{\Delta }h\mathrm{div}\stackrel{}{F}^{j+1/2}+(\stackrel{}{}h\stackrel{}{F}^{j+1/2})\stackrel{}{}(H_0+h))=0.$$
(61)
Once again, the surface of the star is defined by $`H_0+h=0`$. But from the Eq.(60) we have the correct BC: $`\stackrel{}{W}\stackrel{}{}(H_0+h)|_{H_0+h=0}=0`$. The surface of the star is then an unknown quantity that is determined by the relation $`H_0+h=0`$ in solving the Eq.(61) by iteration. This technique was developed in a different context in the already quoted paper gougc01 .

In the linear case, we have $`W_r=0`$ at $`r=1`$ \[cf. Eq.(57)\], and the enthalpy does not vanish on the non-perturbed surface of the star. But once again the free surface of the star can be obtained by looking for the position where the total enthalpy $`H_0+h=0`$ vanishes. This surface coincides within first order quantities with the non perturbed surface $`r=1`$. Note that if $`\gamma >0`$ the pressure $`P|_{r=1}h^{\gamma +1}|_{r=1}`$ vanishes within terms $`o(h)`$.