A&A 408, 331-335 (2003)
DOI: 10.1051/0004-6361:20030917

Differential rotation and hydrodynamic stability of hot neutron stars

V. Urpin

A. F. Ioffe Institute for Physics and Technology, SU 194021, St. Petersburg, Russia
Isaak Newton Institute of Chili, Branch in St. Petersburg, Russia

Received 3 February 2003 / Accepted 22 May 2003

Young neutron stars formed in core collapse or originating in the merger of a binary neutron star system may possess significant differential rotation. We consider the stability properties of such hot (non-superfluid) differentially rotating neutron stars. Compared to ordinary stars, the criteria of hydrodynamic instability of neutron stars can be different because generally the thermal diffusivity and viscosity are of the same order of magnitude and are relatively large. For example, the well-known Goldreich-Schubert instability can manifest itself only in rapidly rotating stars and is strongly suppressed if rotation is slow. The rotational instabilities in neutron stars can grow on the thermal or dynamical timescale.

Key words: MHD - stars: neutron - stars: rotation - stars: pulsars: general

1 Introduction

It is generally believed that young neutron stars formed in the collapse of a massive star should be rapidly rotating due to the conservation of the angular momentum of the collapsing core. Most likely, this rotation is non-uniform because even if the progenitor rotates uniformly, collapse will generate a significant amount of differential rotation. Numerical simulations of rotating core collapse indicate that the angular velocity in a newly formed neutron star should be appreciably non-uniform (see Mönchmeyer & Müller 1989; Janka & Mönchmeyer 1989; Zwerger & Müller 1997; Rampp et al. 1998). During a convective stage that lasts $\sim $30-40 s (Miralles et al. 2000, 2002) and seems to be inevitable in protoneutron stars, the original distribution of the angular velocity may evolve to even a more complex one. During this stage, the angular momentum is basically transported by convection and neutron finger instability and can be very rapid. Perhaps, turbulent Reynolds stresses result in the angular velocity dependent on both the radial and polar coordinates like that of the Sun where the departure from a uniform rotation reaches about 30%. When convection is exhausted and the neutron star becomes transparent to neutrinos it is still rather hot and, therefore, the matter is likely to be in a normal (non-superfluid) state in the core. The non-superfluid neutrino transparent stage can last from hours to thousands of years depending on the model of cooling and superfluidity, and our consideration will concern this stage. Note that accretion of supernova remnant material may also drive differential rotation during this stage (Watts & Andersson 2002). Differential rotation can also be maintained by stellar oscillations. For instance, recent papers on r-modes argued that the modes can drive differential rotation via non-linear effects (Rezzolla et al. 2000; Levin & Ushomirsky 2001).

There are only a few possibilities that can be contemplated as mechanisms of the angular momentum transport during the neutrino-transparent non-superfluid stage: viscosity, hydrodynamic and hydromagnetic instabilities, and large scale motions caused by the magnetic field if the latter has been generated in the core by convection. Viscosity seems to be not very efficient at large scales because it can redistribute the angular momentum only on a relatively long timescale. The characteristic timescale of viscous dissipation, $\tau_{\nu}$, can be estimated as $L^{2}/\nu$ where L is the lengthscale of the angular velocity and $\nu$ is the kinematic viscosity. If the star is transparent to neutrinos, viscosity is mainly determined by the scattering of neutrons and has been calculated in the classical paper by Flowers & Itoh (1976). A simple fitting formula for their numerical results has been obtained by Cutler et al. (1990). Using this fit and assuming that L is comparable to the neutron star radius, we can estimate $\tau_{\nu}$ as 102- 103 yrs that probably is comparable to (or longer than) the duration of a non-superfluid stage.

If there is a sufficiently strong magnetic field in the neutron star then large-scale motions induced by the field can redistribute the angular momentum by magnetic braking (Shapiro 2000). Differential rotation twists lines of a frozen-in poloidal magnetic field and amplifies the toroidal field. This process generates Alfvén waves, which can transport the angular momentum within the star and even carry out some angular momentum from the star to the surrounding plasma. Note, however, that the efficiency of this mechanism is very sensitive to assumptions regarding the density of external plasma, $\rho_{\rm ex}$. For example, a simple estimate of the magnetic braking timescale proposed by Shapiro (2000) shows that this process is more efficient than viscous damping only at $\rho_{\rm ex}/\rho > 10^{-14}$ even if the magnetic field is as strong as 1012 G, where $\rho$ is the density of the neutron star. For smaller $\rho_{\rm ex}$, Alfvén waves inside the neutron star dissipate on a viscous timescale.

The angular momentum can also be transported by turbulent motions if differential rotation is unstable. Since a spatial dependence of the angular velocity is likely to be complex, the neutron star can be subject to hydrodynamic instabilities. The stability properties of differentially rotating neutron stars and radiative zones of ordinary stars can be essentially different because the viscous and thermal timescales are generally relatively short for small-scale disturbances in neutron stars and are generally comparable. This can change qualitatively, for example, the criterion of the well-known Goldreich-Schubert instability (Goldreich & Schubert 1967) that arises on the thermal timescale.

Hydrodynamic instabilities can be of importance for neutron stars originating in the merger of a binary neutron star system as well. In their numerical simulations, Rasio & Shapiro (1999) and Shibata & Uryu (2000) found that coalescence will form a differentially rotating remnant with the core rotating faster than the envelope. Note that differentially rotating neutron stars can support significantly more rest mass than their uniformly rotating counterparts (Baumgarte et al. 2000). Hydrodynamic instabilities in the remnant will destroy differential rotation and may lead to delayed collapse and a delayed gravitational wave burst.

In this paper, we derive the stability criteria for differentially rotating neutron stars. Our consideration assumes the neutron stars that are relatively cold to be transparent to neutrinos but are still sufficiently hot to be in a normal (non-superfluid) state. This evolutionary stage begins at the age $\sim $40 s and can last from hours to $\sim $103 years depending on the model. The paper is organized as follows. In Sect. 2, we derive the dispersion equation governing the rotational modes of a differentially rotating neutron star. In Sect. 3, the conditions of instability are obtained. The growth rate of instability is calculated in Sect. 4. A discussion of the results is given in Sect. 5.

2 Dispersion equation for rotational modes

Consider the core of a neutron star rotating with the angular velocity $\Omega = \Omega(s, z)$, where (s, $\varphi$, z) are cylindrical coordinates. We consider axisymmetric short wavelength perturbations with the space-time dependence $\exp(\gamma t -
i \vec{k} \cdot \vec{r})$ where $\vec{k}= (k_{\rm s}, 0, k_{z})$ is the wavevector, $\vert\vec{k} \cdot \vec{r}\vert \gg 1$. Small perturbations will be indicated by subscript 1, whilst unperturbed quantities will have no subscript, except for indicating vector components when necessary. In the unperturbed state, the star is in hydrostatic equilibrium,

\begin{displaymath}\frac{\nabla p}{\rho} = \vec{G}, \;\;\;
\vec{G} = \vec{g} + \Omega^{2} \vec{s}
\end{displaymath} (1)

where $\vec{g}$ is the gravity.

We use the Boussinesq approximation since the e-folding time of short wavelength instabilities associated with shear is typically much longer than the period of a sound wave with the same wavelength. Then, the linearized equations governing the behaviour of small perturbations in the lowest order in $\vert\vec{k} \cdot \vec{r}\vert^{-1}$ read

$\displaystyle \gamma \vec{V}_{1} + 2 \vec{\Omega} \times \vec{V}_{1} +
...a =
\frac{i \vec{k} p_{1}}{\rho} - \beta \vec{G} T_{1}
- \nu k^{2} \vec{V}_{1},$     (2)

\begin{displaymath}\vec{k} \cdot \vec{V}_{1} = 0 \;,
\end{displaymath} (3)

\begin{displaymath}\gamma T_{1} + \vec{V}_{1} \cdot (\Delta \nabla T) =
- \chi k^{2} T_{1} - \varepsilon_{1}/ \rho c_{\rm p} \;,
\end{displaymath} (4)

where $\vec{V}_{1}$, p1 and T1 are perturbations of the hydrodynamic velocity, pressure and temperature, respectively; $\varepsilon_{1}$ is the perturbation of the neutrino emissivity; $\beta= - (\partial \ln \rho / \partial T)_{\rm p}$ is the thermal expansion coefficient; $\nu$ and $\chi$ are the kinematic viscosity and thermal diffusivity, respectively; $c_{\rm p}$ is the specific heat at constant pressure; $(\Delta \nabla T) = \nabla T - \nabla_{\rm ad} T$ is the difference between the actual and adiabatic temperature gradients; we denote by $\vec{e}_{\varphi}$ the unit vector in the azimuthal direction. In Eq. (2) it is assumed that the density perturbation in the buoyancy force is mainly determined by the temperature perturbation, thus $\rho_{1}= - \rho \beta T_{1}$, in accordance with the idea of the Boussinesq approximation. We take into account viscosity in Eq. (2) because its timescale is generally comparable to the thermal timescale, but we neglect viscous dissipation in Eq. (4) where its contribution is much smaller than advection. Note the absence of terms proportional to p1 in the thermal balance Eq. (4) since their contribution is negligible in the Boussinesq approximation. We took into account the effect of thermal diffusivity since it can be significant for perturbations with short wavelengths.

The neutrino emissivity depends on details of the internal structure but can generally be expressed in terms of $\rho$ and T, $\varepsilon = \varepsilon (\rho, T)$ (see, e.g., Maxwell 1979). Then,

\begin{displaymath}\varepsilon_{1} = \frac{\partial \varepsilon}{\partial T} T_{...
...\rho \frac{\partial \varepsilon}{\partial \rho} \right) T_{1}.
\end{displaymath} (5)

The dispersion equation can now be obtained from Eqs. (2)-(4) in a standard way. Equating the determinant of Eqs. (2)-(4) to zero, we obtain

\begin{displaymath}\gamma^{3} + b_{2} \gamma^{2} + b_{1} \gamma + b_{0} = 0
\end{displaymath} (6)

$\displaystyle b_{2} = \omega_{\rm T} + 2 \omega_{\nu} , \;\;\; b_{1} = Q^{2} +
\omega_{\rm g}^{2} + \omega_{\nu} (\omega_{\nu} + 2 \omega_{\rm T}),$      
$\displaystyle b_{0} = \omega_{\rm T} (Q^{2} + \omega_{\nu}^{2}) + \omega_{\nu}
\omega_{\rm g}^{2}.$                                              

                                               Q2 = $\displaystyle 4 \Omega^{2} \frac{k_{z}^{2}}{k^{2}} + 2 \Omega s
...artial \Omega}{\partial s}
- k_{s} \frac{\partial \Omega}{\partial z} \right) ,$  
$\displaystyle \omega_{\rm g}^{2}$ = $\displaystyle - \beta (\Delta \nabla T) \cdot \left[ \vec{G} -
\frac{\vec{k}}{k^{2}} (\vec{k} \cdot \vec{G}) \right] ,
\;\; \omega_{\nu} = \nu k^{2} ,$  
$\displaystyle \omega_{\rm T}$ = $\displaystyle \omega_{\chi} + \omega_{\varepsilon} , \;\;
\omega_{\chi} = \chi ...
...}{\partial T}
- \beta \rho \frac{\partial \varepsilon}{\partial \rho}
\right) ,$  

$\omega_{\rm g}$ is the frequency of buoyancy waves; $\omega_{\chi}$ and $\omega_{\nu}$ are the inverse time-scales of dissipation due to the thermal conductivity and viscosity, respectively; $\omega_{\varepsilon}$ characterizes the cooling rate due to neutrino emission. In Eq. (6), the terms proportional to Q2 describe the effects associated with the angular velocity and its gradient, the terms containing $\omega_{\rm g}^{2}$ and $\omega_{\chi}$ are caused by the buoyancy force and the energy exchange between the perturbations and surrounding plasma, respectively.

3 The criteria of instability of neutron star cores

Equation (6) describes three low-frequency modes that can exist in rotating non-magnetic neutron star. The sound waves are excluded from consideration in the Boussinesq approximation. The condition that at least one of the roots of Eq. (6) has a positive real part (unstable mode) is equivalent to one of the following inequalities

\begin{displaymath}b_{2} < 0 , \;\;\; b_{1} b_{2} < b_{0} , \;\;\; b_{0} < 0
\end{displaymath} (7)

being fulfilled (see, e.g., Aleksandrov et al. 1985; DiStefano III et al. 1994). Since $\omega_{\rm T}$ and $\omega_{\nu}$ are positive defined quantities, the first condition, b2 < 0, will never apply, and only the two other conditions determine the instability of neutron stars.

Substituting the values of b0, b1 and b2, we can rewrite the second condition (7) as

\begin{displaymath}(\omega_{\rm T} + \omega_{\nu}) \omega_{\rm g}^{2} + 2 \omega_{\nu}
[Q^{2} + (\omega_{\rm T} + \omega_{\nu})^{2}] < 0.
\end{displaymath} (8)

In the case of an inviscid fluid with $\Omega=0$ and $\vec{g} \parallel
\Delta \nabla T$, Eq. (8) reduces to the well-known Schwarzschild criterion of convection,

\begin{displaymath}\omega_{\rm g}^{2} < 0.
\end{displaymath} (9)

The third criterion (7), b0 < 0, is equivalent to

\begin{displaymath}\omega_{\rm T} (Q^{2}+ \omega_{\nu}^{2}) + \omega_{\nu} \omega_{\rm g}^{2} < 0.
\end{displaymath} (10)

If viscosity is negligible then we recover the well-known criterion of the Goldreich-Schubert instability,

Q2 < 0, (11)


\begin{displaymath}\frac{k_{z}^{2}}{k^{2}} \cdot \Omega_{\rm e}^{2} -
...{2}} \cdot 2 \Omega s \frac{\partial \Omega}
{\partial z} < 0,
\end{displaymath} (12)

where we denote

\begin{displaymath}\Omega_{\rm e}^{2} = \frac{1}{s^{3}} \frac{\partial}{\partial s}(s^{4}
\end{displaymath} (13)

If $\partial \Omega / \partial z \neq 0$ then there exists a region in the plane $(k_{\rm s}, k_{z})$ where the condition Q2 < 0 is fulfilled for any dependence of $\Omega$ on s and z. The components of a wave vector should satisfy the inequality

\begin{displaymath}\left\vert \frac{k_{s}}{k_{z}} \right\vert > \frac{\vert\part...
...rt}{\vert 2 \Omega s^{4} \partial \Omega /
\partial z \vert},
\end{displaymath} (14)

and ks and kz are of the same/opposite sign if the quantity $[\partial (s^{4} \Omega^{2}) / \partial s]/ [\partial
\Omega/ \partial z]$ is positive/negative.

Consider the criteria (8) and (10) in more detail. In neutron stars transparent to neutrinos, we have typically $\omega_{\nu}
\sim \omega_{\rm T}$. Since $\omega_{\nu}$ and $\omega_{\rm T}$ increase rapidly when the wavelength decreases, neither of the conditions (8) or (10) can be satisfied for very short wavelengths. Therefore, we consider Eqs. (8) and (10) at

\begin{displaymath}\omega_{\nu} \sim \omega_{\rm T} \ll {\max}(\vert\omega_{\rm g}\vert, \vert Q\vert).
\end{displaymath} (15)

Substituting the fitting expression for $\nu$ (Cutler et al. 1990) and estimating $\vert Q\vert \sim 10^{3}$ s-1, we obtain from Eq. (15) the following inequality for a wavelength, $\lambda = 2 \pi /k$,

\begin{displaymath}\lambda \gg \lambda_{\rm c} = 2.1 \; \frac{\rho_{14}^{5/8}}{T_{9}} \;\;
{\rm cm},
\end{displaymath} (16)

where $\rho_{14}$ has units 1014 g/cm3, and T9 has units 109 K. Taking into account that usually $\omega_{\chi}
\gg \omega_{\varepsilon}$, Eqs. (8) and (10) can be transformed at $\lambda \gg \lambda_{\rm c}$ into

\begin{displaymath}Q^{2} < - \frac{1}{2} \left(1 + \frac{\chi}{\nu} \right) \omega_{\rm g}^{2},
\end{displaymath} (17)


\begin{displaymath}Q^{2} < - \frac{\nu}{\chi} \omega_{\rm g}^{2},
\end{displaymath} (18)

respectively. If the star is far from the rotational distortion ( $g \gg \Omega^{2}s$) then $\Delta \nabla T$ and $\vec{g}$ are approximately radial, and

\begin{displaymath}\omega_{\rm g}^{2} \approx - \beta \; \frac{k_{\theta}^{2}}{k...
...\nabla T) = \frac{k_{\theta}^{2}}{k^{2}}
\omega_{\rm BV}^{2},
\end{displaymath} (19)

where $\omega_{\rm BV}^{2} = - \beta ( \vec{g} \cdot \Delta
\nabla T)$ is the Brunt-Vaisala frequency, and $k_{\theta}$ is the $\theta$-component of a wavevector, $k_{\theta}^{2}= k_{s}^{2} \cos^{2} \theta - 2 k_{s} k_{z} \cos
\theta \sin \theta + k_{z}^{2} \sin^{2} \theta$, $\theta$ is the polar angle;

The conditions (17) and (18) depend on the direction of $\vec{k}$. To obtain the true necessary conditions of instability, we have to minimize these expressions as functions of $\vec{k}$. For both the inequalities (17) and (18), a dependence on $\vec{k}$ has the same form and can be represented as

F(ks, kz) = A ks2 - C kz ks + D kz2 < 0, (20)

where A, C and D are functions of $\omega_{\rm BV}^{2}$, Q2, and $\theta$. The left hand side reaches its minimum at kz = C ks/ 2 B. The minimum value is

\begin{displaymath}\min F(k_{s}, k_{z}) = k_{s}^{2} \left( A - \frac{C^{2}}{4 D}
\end{displaymath} (21)

This minimum is negative if

C2 > 4 A D. (22)

Substituting the values A, C and D, corresponding to the condition (17), into Eq. (22), we obtain the following criterion

\begin{displaymath}\left( \frac{\partial \Omega}{\partial z} \right)^{2} >
...\partial \Omega}{\partial z} \cos \theta
\sin \theta \right),
\end{displaymath} (23)

where $\mu = (1 + \chi/\nu)^{1/2}$, $q = \omega_{\rm BV}/ \Omega$, and $\alpha = \mu^{2} q^{2}/2$. If we now subsitute into Eq. (22) the coefficients A, C and D, corresponding to the condition (18), then we recover Eq. (23) with $\alpha =
\eta q^{2}$, where $\eta = \nu / \chi$. Therefore, a further consideration of instabilities given by Eqs. (17) and (18) can be done in parallel. We assume that the Rayleigh stability criterion is usually fulfilled in neutron stars and $\Omega_{\rm e}^{2}
> 0$.

Solving the inequality (23) for $\partial \Omega / \partial z$, we obtain that Eq. (23) is fulfilled if

\begin{displaymath}\frac{\partial \Omega}{\partial z} > \frac{\alpha \Omega}{s}
...}^{2}}{\alpha \Omega^{2}} \right)^{1/2} -
\sin \theta \right],
\end{displaymath} (24)


\begin{displaymath}\frac{\partial \Omega}{\partial z} < -\frac{\alpha \Omega}{s}...
...}^{2}}{\alpha \Omega^{2}} \right)^{1/2} +
\sin \theta \right].
\end{displaymath} (25)

Note that the conditions (24) and (25) are qualitativy different to the classical Goldreich-Schubert criterion that allows the instability at any vertical gradient of the angular velocity. In neutron stars, $\partial \Omega / \partial z$ should exceed a some threshold value to start the instability. This threshold value depends on the parameter $\alpha$ and can generally vary over a wide range.

In the case of a rapidly rotating star, $\Omega \gg \omega_{\rm BV}$, we have $q^{2} \ll 1$ and $\alpha \ll 1$, and Eqs. (24) and (25) yield

\begin{displaymath}\left\vert \frac{\partial \Omega}{\partial z} \right\vert >
...rt\cos \theta\vert \; \frac{\sqrt{\Omega_{\rm e}^{2}}}{s}\cdot
\end{displaymath} (26)

Since $\alpha$ is small even a relatively small vertical gradient of $\Omega$ can cause instability in rapidly rotating stars. However, the situation is quite different in slowly rotating stars with $\Omega < \omega_{\rm BV}$ ( $\alpha \gg 1$). For such stars, Eqs. (24) and (25) yield

\begin{displaymath}2 \Omega \vert \cos \theta \sin \theta \vert \frac{\partial \...
...\theta}{s^{4}} \frac{\partial}{\partial s}
(s^{4} \Omega^{2}),
\end{displaymath} (27)


\begin{displaymath}\frac{\partial \Omega}{\partial z} < - \frac{\alpha \Omega}{s}
\vert \sin \theta \cos \theta \vert.
\end{displaymath} (28)

The inequality (27) requires a large positive vertical gradient of $\Omega$ whereas the condition (28) can be satisfied only at very large negative gradients. It is unlikely that these sorts of profiles of the angular momentum can be formed in neutron stars and, most likely, the rotational instabilities cannot manifest themselves in slowly rotating stars.

4 The growth rate of instabilities

Consider the growth rate of rotational instabilities. Since the coefficients of Eq. (6) are real there exist three real roots or one real and two complex conjugate roots. The number of roots with a positive real part is determined by Routh criterium (DiStefano III et al. 1994), which states that the number of unstable modes of a cubic Eq. (6) is given by the number of changes of sign in the sequence

\begin{displaymath}\left\{ 1, \; b_{2}, \; \frac{b_{2} b_{1} - b_{0}}{b_{2}},
\; b_{0} \right\}\cdot
\end{displaymath} (29)

Therefore, if the condition (8) holds but the condition (10) does not then the two complex modes are unstable. On the contrary, if the inequality (10) is fulfilled but Eq. (8) is not then only one real mode is unstable. In the case when both the inequalities (8) and (10) hold, only one real mode is again unstable.

The roots $\gamma_{i}$ (i=1,2,3) of the cubic Eq. (6) can be represented as $\gamma_{i} = x_{i} - b_{2}/3$. The expressions for xi are

\begin{displaymath}x_{1} = u+v, \;\; x_{2}=
\varepsilon_{1} u + \varepsilon_{2} v, \;\;
x_{3} = \varepsilon_{2} u +\varepsilon_{1} v ,
\end{displaymath} (30)


\begin{displaymath}\varepsilon_{1, 2}= - \frac{1}{2} \pm i \frac{\sqrt{3}}{2},
\;\; (u, v) = \Big(-q \pm \sqrt{q^{2}+p^{3}}\Big)^{1/3},
\end{displaymath} (31)

                                      q = $\displaystyle \frac{\omega_{\rm T}}{6} \left( 2 Q^{2} - \omega_{\rm g}^{2} +
...frac{26}{9} \omega_{\nu}^{2} -
\frac{10}{9} \omega_{\rm T} \omega_{\nu} \right)$  
    $\displaystyle +\frac{\omega_{\nu}}{6} \left( \omega_{\rm g}^{2} - 2 Q^{2}
+ \fr...
...frac{20}{9} \omega_{\nu} \omega_{\rm T}
- \frac{2}{9} \omega_{\nu}^{2} \right),$  
p = $\displaystyle \frac{1}{3}
\left[ Q^{2} + \omega_{\rm g}^{2} - \frac{1}{3}( \omega_{\rm T} -
\omega_{\nu} )^{2} \right]$  

(see, e.g., Bronstein & Semendyaev 1957). These expressions are rather cumbersome but they can be simplified in the most interesting case (15).

Under the condition (15) (or (16)), we have from the expressions (30), (31) with the accuracy in terms $\propto$ $\omega_{\chi} \sim

\begin{displaymath}\gamma_{1} \approx - \frac{\omega_{\rm T} Q^{2}+ \omega_{\nu}
\omega_{\rm g}^{2}}{Q^{2} +
\omega_{\rm g}^{2}} ,
\end{displaymath} (32)

\begin{displaymath}\gamma_{2,3} \approx \pm i \sqrt{Q^{2} +
\omega_{\rm g}^{2}} ...
...}+ 2 \omega_{\nu}
(Q^{2} + \omega_{\rm g}^{2})} \cdot
\end{displaymath} (33)

Note that the second term on the right hand side of Eq. (33) is small compared to the first one. Which one of these three modes is unstable depends on the conditions (17) and (18) and the sign of $Q^{2}+\omega_{\rm g}^{2}$. This quantity is negative if

\begin{displaymath}\left( \frac{\partial \Omega}{\partial z} \right)^{2} >
...\partial \Omega}{\partial z} \cos \theta
\sin \theta \right),
\end{displaymath} (34)

and positive in the opposite case. If the condition (18) is fulfilled but condition (17) is not then the third mode is unstable at $Q^{2} + \omega_{\rm g}^{2} < 0$, and the first one at $Q^{2} + \omega_{\rm g}^{2} > 0$. If the condition (17) is satisfied but the condition (18) is not then the first and third modes are unstable at $Q^{2} + \omega_{\rm g}^{2} < 0$, and the second and third modes are unstable at $Q^{2} + \omega_{\rm g}^{2} > 0$. Finaly, if the both conditions (17) and (18) are fulfilled then again the third mode is unstable at $Q^{2} + \omega_{\rm g}^{2} < 0$, and the first one at $Q^{2} + \omega_{\rm g}^{2} > 0$.

Contrary to the classical Goldreich-Schubert instability that arises on the thermal time scale, the rotational instabilities in neutron stars can grow either on the thermal or dynamical timescale. If $Q^{2} + \omega_{\rm g}^{2} > 0$ then the growth rate is

\begin{displaymath}Re \gamma \sim \omega_{\chi} \sim \omega_{\nu}.
\end{displaymath} (35)

If $Q^{2} + \omega_{\rm g}^{2} < 0$ (the inequality (34) is fulfilled) then instability is typically much faster and grows on the dynamical time scale

\begin{displaymath}Re \gamma \sim \sqrt{\vert Q^{2} + \omega_{\rm g}^{2}\vert}.
\end{displaymath} (36)

The necessary condition of such a regime can be obtained by comparing the inequality (34) from one side and the inequalities (17) and (18) from another side. Then, the dynamical regime of instability exists if

\begin{displaymath}\max (\mu^{2}/2, \eta) > 1
\end{displaymath} (37)

that can generally be fulfilled in neutron stars.

5 Conclusion

We considered hydrodynamic stability of differentially rotating neutron stars during the evolutionary stage when they are non-superfluid and transparent to neutrinos. Differentially rotating neutron stars can be formed either in core collapse or in the merger of binary neutron stars. If the Rayleigh stability criterion is fulfilled then the rotational instabilities can arise only if the angular velocity depends on z. However, the instabilities caused by this dependence in neutron stars turns out to be qualitatively different from those in ordinary stars. The classical Goldreich-Schubert instability that likely operates in ordinary stars can arise at any $\partial \Omega / \partial z \neq 0$. In neutron stars, the rotational instabilities can appear only if $\vert \partial \Omega /
\partial z\vert$ exceeds some threshold value. This value depends on the angular velocity of a neutron star and can be small for rapidly rotating stars with $\Omega > \omega_{\rm BV}$ or large for slowly rotating stars with $\Omega < \omega_{\rm BV}$. The difference to ordinary stars is caused by a particular character of kinetic processes in neutron stars where viscosity is large and generally comparable to the thermal diffusivity. Unfortunately, calculations of the Brunt-Vaisala frequency, $\omega_{\rm BV}$, in young neutron stars have some problems because it depends on a poorly known equation of state. However, most likely this frequency is of the order of 102-103 s-1. Therefore, the critical value of the period, discriminating between the rapid unstable and slow stable rotation in neutron stars, is approximately 0.1-0.01 s.

Contrary to the classical Goldreich-Schubert instability, the rotational instabilities in neutron stars can arise either on the thermal or dynamical timescale. If the condition (37) is fulfilled then the growth time of instability is comparable to the rotation period. Otherwise, the growth rate is determined by the thermal effects, and the corresponding growth time is $\sim $ $10^{-3}
\rho_{14}^{-5/4} T_{9}^{2} \lambda^{2}$. For perturbations with $\lambda =1$ km, the growth time is $\sim $100 days, whereas perturbations with $\lambda =1$ m grow on a timescale $\sim $10 s.

In the present paper, we have addressed the behaviour only of axisymmetric perturbations. It is clear, however, that the obtained results can apply to nonaxisymmetric perturbations with the azimuthal wavelength much longer than the vertical or radial ones, $\min (k_{r},
k_{z}) \gg k_{\varphi}$. Note that, generally,non-axisymmetric perturbations can also undergo other instabilities (e.g., the r-mode instability, see Andersson et al. 1999). Interaction of the neutron star crust with turbulent motions caused by rotational instabilities may result in small irregularities of the measured periods of young pulsars.


This work has been supported by the Spanish Ministery of Science and Technology (grant AYA2001-3490-C02-02).


Copyright ESO 2003