EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 496 / No 3 (March IV 2009) A&A, 496 3 (2009) 609-617 Full HTML
Free Access
Issue
A&A
Volume 496, Number 3, March IV 2009
Page(s) 609 - 617
Section Astrophysical processes
DOI https://doi.org/10.1051/0004-6361/200810793
Published online 09 February 2009

Angular momentum transport during X-ray bursts in neutron stars: a numerical general relativistic hydrodynamical study

A. Hujeirat1 - F.-K. Thielemann2


1 - ZAH, Landessternwarte Heidelberg-Königstuhl, Universität Heidelberg, 69120 Heidelberg, Germany
2 - Departement Physik, Universität Basel, Switzerland

Received 13 August 2008 / Accepted 18 January 2009

Abstract
Aims. The distribution of angular momentum of matter during X-ray bursts in neutron stars is studied by means of axi-symmetric general relativistic hydrodynamics.
Methods. The set of fully general relativistic Navier-Stokes equations is solved implicitly using the implicit solver GR-I-RMHD in combination with a third order spatial and second order temporal advection scheme. The viscous operators are formulated using a Kerr-like metric in the fixed background of a slowly rotating neutron star whose radius coincides with the corresponding last stable orbit. The importance of these operators and their possible simplifications are discussed as well. To verify the consistency and accuracy of the solution procedure, the time-dependent evolution of non-rotating heat bubbles during their rise to the surface of a white dwarf are followed and compared with previous calculations.
Results. In the rotating case and depending on the viscosity parameter, $\alpha _{\rm {tur}}$, it is found that the viscously-initiated fronts at the center of bursts propagate at much faster speed than the fluid motion. These fast fronts act to decouple angular momentum from matter: angular momentum is transported outwards while matter sinks inwards into the deep gravitational well of the neutron star, thereby enhancing the compression of matter necessary to initiate ignition, that subsequently spreads over the whole surface of the neutron star on the viscous time scale. Based on the numerical simulations, we find that a viscosity parameter $\alpha_{\rm {tur}} = \mathcal{O}(0.1)$ is most suitable for fitting observations of neutron stars during X-ray bursts. It is argued that the spin up observed in the cooling tails of X-ray bursts is a transient phase, which eventually should be followed by a spin down phase. This delay can be attributed to a significant lengthening of the viscous time scale due to rapid cooling of matter in the outer layers.

Key words: stars: neutron - X-rays: bursts - relativity - hydrodynamics - methods: numerical - stars: rotation

1 Introduction

Some observed neutron stars (NSs) belong to the family of ultra-compact objects, in which general relativistic effects are prominent (Stergioulas 2003; Psaltis 2008; Shapiro & Teukolsky 1983). Depending on the equation of state, NSs may live even inside their last stable orbits, making the conversion efficiency of gravitational energy into radiation even higher than that of accreting Schwarzschild black holes (BHs, Camenzind 2007).

Neutron stars in low-mass X-ray binaries accrete material from their companions via disks. The angular momentum associated with the accreted matter is capable of spinning up NSs, possibly making them millisecond pulsars (Bhattacharya & van den Heuvel 1991). Aside from spinning, the hydrogen and/or helium-rich material accumulated over hours to days is expected to spread over the surface of the NS and form a thin and highly compressed shell. The matter in this shell is thermonuclear unstable and may easily undergo a runaway ignition that results in liberating energy of the order of 1039-1040 erg in less than 100 s (Lewin et al. 1995). Such thermonuclear explosions are believed to be the origin of the observed X-ray bursts on accreting NSs (for reviews, see Cumming 2004, and the references therein).

However, following the dynamical evolution of the accreted matter from the disk, through the boundary layer down to the mixing layer, its spread over the surface and finally into the ignition phase is a complicated task and would require another decade of additional investigations. Of particular importance is understanding how the disk matter reaches and interacts with the very outer envelopes of strongly or weakly magnetized NSs, the rate of spreading of accreted matter over the surface of the NS, the role and efficiency of the magnetic-rotational instability (MRI) in re-distributing angular momentum in the boundary layer (BL) as well as the dissipation rate of rotational kinetic energy. Also, the way the shear in the low-density BL viscous-heats and chemically mixes the matter with the matter in the deeper layers, taking into account magnetic-generated turbulence, heat and radiation transfer is of vital importance for understanding the mechanisms leading to X-ray burst events (Piro & Bildsten 2007).

BLs, on the other hand, are considered to be optimal regions where the bulk of the kinetic energy of the inflowing matter is re-directed into powering outflows and jets. Strong shearing and enhanced dissipation in combination with appropriate topological changes of the magnetic fields in the BL may provide the mechanisms required for formation and acceleration of powerful jets from around compact objects (Chakrabarti 2001). Indeed, numerous highly relativistic jets have been observed to also emanate from around accreting NSs with bulk Lorentz factors that are comparable to those emanating from microquasars (Migliari 2008,2006).

At the simulation site, numerous general relativistic hydrodynamical calculations have been carried out to study different aspects of NSs, such as formation, merger, inner structure, accretion or jets around NSs (Thielemann 1990; Abdikamalov et al. 2008; Liebendörfer et al. 2002; Özel & Psaltis 2003; Marti & Müller 2003; Shibata 2003; Shibata et al. 2006; McKinney 2006), see also the references therein.

Duez et al. (2004) carried out general relativistic calculations to study the formation of hypermassive NSs, taking into account the effect of viscosity. The authors found that viscosity drives the NS's inner core into rigid rotation and simultaneously transports angular momentum outwards into the outer envelope. As a consequence, the core is found to contract in a quasi-stationary manner while the outer layers expand to form a differentially rotating torus.

This behavior is similar to accretion of matter via rotating disks. Here the viscosity acts to decouple matter from angular momentum, in that it transports angular momentum outwards, while forcing the matter to sink deeper into the gravitational well of the central object (Pringle 1981). In the absence of viscosity, angular momentum as well as magnetic fields in ideal MHD are frozen-in to the matter. Thus, while strong magnetic fields are essential to enable rigid rotation, viscosity, on the other hand, drives the outer layers into differential rotation.

The type of rotation in the outer layers, differential or rigid rotation, may have profound effects on the conditions leading to the X-ray bursts observed on NSs (Spitkovsky et al. 2002; Bildsten & Strohmayer 1999).

Indeed, recent observations of X-ray bursts revealed so called burst oscillations, in which a spin-up or spin-down of the NSs in their cooling tails have been detected, reaching a plateau on the asymptotic limit (see Strohmayer 2001a; Strohmayer & Markwardt 1999; Strohmayer et al. 1997, for a detailed discussion). It has been argued that the increase/decrease of the spin of NSs during bursts is connected to the redistribution of angular momentum of the thermonuclear shell (Strohmayer 2001b). Accordingly, when a thermonuclear shell starts to expand at the burst onset, the moment of inertia increases while its spin decreases. When the shell starts to contract and subsequently recouples to the NS, the inertia decreases and the spin increases.

Also, several X-ray bursts on NSs display a spin down rather than spin up in their cooling tails (Strohmayer 1999b). In this case, however, it was suggested that the spin down probably begins episode of thermonuclear energy release, most likely in layers underlying those responsible for the initial runaway.

Cumming & Bildsten (2003), however, investigated in detail the hydrostatic expansion during bursts and the expected change of spin due to angular momentum conservation and concluded that a shell expansion/contraction alone cannot explain the mechanisms underlying the observed spin-up/down of the NSs during bursts. The model in which ignition starts at a point and spreads over the whole surface of the NS via burning fronts appears to fit observations, which reveals that the X-ray emitting area increases during the bursts (Strohmayer et al. 1997). However, the role of rotation, the nature of these burning fronts and the manner in which they affect their surrounding are poorly understood.

In this paper we present a first attempt to model the rotational evolution of thermally induced bursts beneath the atmosphere of a rotating NS and to study the role and effects of the viscosity on the redistribution of angular momentum under strong gravitational field conditions. Our investigation relies on employing a general relativistic hydrodynamical solver, in which turbulent-eddies have the effect of friction that gives rise to an enhanced re-distribution of angular momentum.

The paper is organized as follows: In Sect. 2 we describe the additional viscous operators that have been incorporated into the solver to study the viscous-redistribution of angular momentum. The results of several model calculations aimed at studying the distribution of angular momentum during X-ray bursts on NSs are presented and discussed in Sect. 3, while in Sect. 4 the results are summarized.

2 The general relativistic Navier-Stokes equations

The set of the general relativistic hydrodynamical equations and their derivations are well described in Sect. 2 of Hujeirat et al. (2008). In this section we list the viscous operators of the momentum equations, which we have incorporated into the implicit solver.

The stress energy tensor for viscous flows has the following form (Richardson & Chung 2002; Camenzind 2007; Font 2003):

\begin{displaymath}{{\rm T}^{\mu\nu}= \rm {{T_{\rm {PF}}^{\mu\nu}} + \left\{\und... ...ar{\sigma}^{\mu\nu} + \frac{\Theta}{3}h^{\mu\nu}]} \right\}}, \end{displaymath} (1)

where $\mu,\nu$ are indices that correspond to the four coordinates $\{t,r,\theta,\varphi\}$ and $\rm {T_{\rm {PF}}^{\mu\nu}}$, ${\rm T}_{\rm {Vis}}^{\mu \nu}$ denote the stress energy tensor due to perfect and viscous flows, respectively. P, $\eta,~\Theta$, are the pressure, which is calculated from the equation of state corresponding to a polytropic or to an ideal gas, the dynamical viscosity which is assumed to be identical to the shear viscosity, and ${\Theta (\doteq\nabla_{\mu} u ^\mu)},$ which measures the divergence or convergence of the fluid world lines, respectively. ${h^{\mu\nu}= u^\mu u^\nu + g^{\mu\nu}}$ is the spatial projection tensor, whereas $\bar{\sigma}$ corresponds to the symmetric spatial shear tensor: \({ \bar{\sigma}^{\mu\nu} = \nabla_{\varsigma} u^\mu h^{\varsigma\nu} + \nabla_{\varsigma} u^\nu h^{\varsigma\mu}}\).

For the X-ray burst calculations, the general relativistic Navier-Stokes equations are solved using the Boyer-Lindquist coordinates in the background of a slowly rotating NS. We use the Kerr metric to describe the spacetime outside the NS. The spin ``a'' in the Kerr metric is set to describe the slow rotation $\Omega_{\rm NS}$ of the NS, whose radius is set to be located far away outside its corresponding event horizon. Therefore, this validates the use of this metric to describe the dynamics of the matter in the atmosphere which is located outside the NS, especially as the atmosphere has a negligibly small mass compared to the mass of the enclosed degenerate core. We note that the large matter-density in the core forces the geometry to deviate considerably from flat space. Therefore, studying the dynamics inside and outside the core necessitates the construction of a combined metric in a manner similar to the Hartle formalism (Hartle 1967).

The elements of the metric read as follows:

\begin{displaymath}g_{\mu\nu} = \left[ \begin{array}{cccc} g_{\rm tt} & 0 & 0 ... ...rphi {\rm t}}& 0 & 0 & g_{\varphi\varphi} \end{array}\right], \end{displaymath}

  where   

\begin{displaymath}\left\{ \begin{array}{lll} g_{\rm tt} &= & \beta_{\varphi}\b... ...^2\cos{\theta} = \alpha \sqrt{\Upsilon}. \end{array}\right . \end{displaymath} (2)

In this formulation, the parameter `` $\Omega_{\rm NS}$'' denotes the spin of the neutron star, which is taken to be much smaller than the break-up frequency. $\beta^{\varphi}$ is the frame-dragging frequency associated with the rotation of the NS:

\begin{displaymath}\beta^{\varphi} = \Omega_{\rm {FD}} = \left[\frac{2G}{c^2}\ri... ...}\right]\left(\frac{1}{r^3}\right)R^2_{\rm NS}\Omega_{\rm NS}, \end{displaymath} (3)

where $J_* = I_* \Omega_{\rm {NS}}$, and ${I_* = \frac{2}{5} M_{\rm NS} R^2_{\rm NS}}$ is the moment of inertia of the NS. The parameters: ${c,~ M_{\rm {NS}},~ G,~r_{\rm g} (= \frac{GM_{\rm {NS}}}{c^2}),~ \alpha}$ denote the speed of light, mass of the NS, the gravitational constant, the gravitational radius and the lapse function, respectively. In writing these expressions, we made use of the coordinate transformation $\bar{\theta} = \pi/2{-}\theta$, where we use the latitude $\theta$ instead of the polar distance angle $\bar{\theta}$; hence the appearance of ``$\cos$'' instead of ``$\sin$'' in the metric terms.

The set of general relativistic Navier-Stokes equations in axi-symmetry can be written as the residual vector equation:

R = 0. (4)

The components of this vector read as follows:
1.
The continuity equation

\begin{displaymath}R_1 = {\partial_t D} + L1_{r~\theta} D = 0 \end{displaymath} (5)

2.
The radial momentum equation

\begin{displaymath}R_2= {\partial_t M_r} + L1_{r~\theta} M_r - f_r - L2^r_{r~\theta} M_r = 0 \end{displaymath} (6)

3.
The vertical momentum equation

\begin{displaymath}R_3= {\partial_t M_{\theta}} + L1_{r~\theta} M_{\theta} - f_{\theta} - L2^\theta_{r~\theta} M_{\theta} = 0 \end{displaymath} (7)

4.
The angular momentum equation

\begin{displaymath}R_4= {\partial_t M_{\varphi}} + L1_{r~\theta} M_{\varphi} - f_{\varphi} - L2^\varphi_{r~\theta} M_{\varphi} = 0 \end{displaymath} (8)

5.
The internal energy equation

\begin{displaymath}R_5 = {\partial_t {\cal{E}}^{\rm d}} + L1_{r~\theta} {\cal{E}... ...\frac{\partial u^t}{\partial t} + L1_{r\theta} u^t\right] = 0, \end{displaymath} (9)

where L $1_{r~\theta}$ are first order advection operators that have the form:

\begin{displaymath}\begin{array}{lll} {\rm L}1_{r~\theta} ~q~ &=& {(-g)^{-1/2}} ... ...ot q V^r + \bar{\nabla}_{\theta} \cdot q V^\theta. \end{array}\end{displaymath}

$f_{r,~\theta,~\varphi}$ are force terms that include pressure gradients, centrifugal and gravitational forces acting along the radial, polar and azimuthal directions, respectively. ${D} (\doteq \rho u^t)$ is the modified relativistic mass density. ${M}_{\mu}$ are the four-momenta: \((M_t,M_r,M_{\theta},M_{\varphi}) \doteq {{\overline{\rm D}}}(u_t,u_r,u_{\theta},u_{\varphi}),\) where \( {{\overline{\rm D}}} \doteq D h,\) and ut is the time-like velocity, \( V^\mu= u^\mu/u^t\) is the transport velocity and ``h'' denotes the enthalpy. \(L2^{\xi}_{r\theta}\) are the spatial projections of the viscous stress energy tensor $\rm {T^{\mu\nu}_{Vis}}~$ (see Eq. (1)) in the respective direction. These are obtained from the projection of the viscous tensor along the vector normal to the hyperspace, i.e., constant in time:

\begin{displaymath}L2^{\xi}_{r\theta} = \nabla_{\mu} {\rm T}_{\rm {Vis}}^{\mu\xi... ... + \Gamma_{\mu\lambda}^{\xi} {\rm T}_{\rm {Vis}}^{\mu\lambda}, \end{displaymath}

where $\xi=\{r,\theta,\varphi\}$. $\nabla_\mu$ corresponds to the spatial divergence of a tensor taken in the Boyer-Lindquist coordinates and $\Gamma_{\mu\lambda}^{\xi}$ are the Christoffel's symbols of the second kind.

For completeness, we re-write the forms of these second order viscous operators explicitly as follows:

\begin{displaymath}\begin{array}{llll} {L2^{r}_{r\theta}} & = &{\bar{\nabla}_r \... ...ot u^\theta\right)~\left(u_r u^r + 1\right) \Big]}, \end{array}\end{displaymath} (10)


\begin{displaymath}\begin{array}{llll} {L2^{\theta}_{r\theta}} & = &\bar{\nabla... ...heta\right)~(u^\theta)^2 + g^{\theta\theta})\Big]}, \end{array}\end{displaymath} (11)


\begin{displaymath}\begin{array}{llll} {L2^{\varphi}_{r\theta}} & = {\bar{\nabla... ...\varphi u^\varphi +1) \Big]},&\!\!\!\!&\!\!\!\! \\ \end{array}\end{displaymath} (12)

where $g_{\mu\nu}$ is the covariant form of the metric tensor  $g^{\mu\nu}$.

2.1 Simplifying considerations

Most of the above-described collection of viscous terms contains highly non-linear, first and second order operators, some of which are difficult to handle numerically and in most cases enlarge the band-width of the coefficient matrix, while others may decelerate, rather than accelerate, the convergence of the numerical solution procedure. In axi-symmetry, few of these terms can be simplified or even neglected without violating the physical consistency of the numerical scheme.

To outline our simplification strategy, we first mention the following relevant issues:

1.
in most astrophysical fluid flows the molecular viscosity is too small to be relevant on observationally reasonable time scales. This implies that the corresponding Reynolds number, which expresses the ratio of inertial to viscous forces, is too large and cannot be captured by solving the hydrodynamical equations numerically. Thus, in the absence of other forms of viscosity, e.g., turbulent viscosity, the above operators have negligible effects on the dynamical evolution of the flow;
moreover, these operators must vanish asymptotically whenever the fluid velocity approaches the speed of light;
2.
turbulent viscosity is more common in modeling astrophysical fluid dynamics. In rotating astrophysical flows, turbulent viscosity is a fundamentally important mechanism for angular momentum transport. Therefore, the viscous operators of the angular momentum equation are important and should converge to the usual Newtonian form whenever the velocity becomes sub-relativistic;

3.
the viscous operators appearing in the radial and vertical momentum equations act, in general, to diffuse and smooth strong velocity gradients. However, the mixed derivatives[*] appearing in these two equations are numerically difficult to handel and their inclusion in the implicit solution procedure may significantly lower the efficiency of the solver. In particular, they enlarge the band width of the corresponding Jacobian, make it difficult to find an appropriate and easy-to-invert pre-conditioners and subsequently increase the computational costs. On the other hand, the physical effect of mixed derivatives in fluid dynamics is generally small or even negligible, depending on the properties of the flow. For example, mixed derivatives have vanishingly small effects in advection-dominated plasmas due to the low viscous interaction. In turbulent dissipative flows however, they act to mainly enhance the viscous coupling between the velocity components. However, as the turbulence in most gravitationally bound astrophysical flows is generally subsonic, the coupling between the momentum equations is predominated by the pressure gradients. This implies that the mixed derivatives in the hydro-equations describing the dynamics of subsonic turbulent flow, can, to first order approximations, be neglected. From the numerical point of view, neglecting such operators is necessary to enhance the sparseness of the Jacobian and to enable the use of pre-conditioners that are based on the directional splitting strategy.
Therefore, our simplification strategy relies on considering only second order, mixed-free and Laplace-like operators. In numerical matrix algebra, such operators generally enhance the diagonal dominance of the coefficient matrix and stabilize its inversion procedure.

Specifically, using our time-implicit formulation, the following viscous operators are used in the numerical solution procedure:

\begin{displaymath}\begin{array}{llll} {\tilde{L2}^{r}_{r\theta}} & = &{\bar{\n... ...ta\right)^2 + g^{\theta\theta}\right)\right\}} \\ \end{array}\end{displaymath} (13)


\begin{displaymath}\begin{array}{llll} {\tilde{L2}^{\theta}_{r\theta}} & = & {\... ...(u^\theta)^2 + g^{\theta\theta}\right]\right\}} \\ \end{array}\end{displaymath} (14)


\begin{displaymath}\begin{array}{llll} {\tilde{L2}^{\varphi}_{r\theta}} & = &{\... ...u^\varphi (u_\varphi u^\varphi +1) ] }. \!\!\!\!\\ \end{array}\end{displaymath} (15)

We note that in carrying out these simplifications, care should be taken to still recover the classical non-relativistic form of the Navier-Stokes equations.

Indeed, it can be easily verified that in the non-relativistic regime the radial component of the diffusion operator ${\tilde{L2}^{\varphi}_{\rm {r\varphi}}}$ reduces to the classical Newtonian form:

\begin{displaymath}{\tilde{L2}^{\varphi}_{r} = \frac{1}{r^2}\frac{\partial}{\partial r} r^4 \eta \frac{\partial\Omega}{\partial r}}, \end{displaymath} (16)

where $\eta = \rho \nu$ and $\nu$ denotes the kinematic viscosity.

2.2 Viscosity prescription

Similar to classical accretion disks, we assume that molecular viscosity is too small to have a significant effect on the angular momentum distribution on a short time scale such as the thermonuclear one. Therefore, we adopt the turbulent viscosity prescription:

\begin{displaymath}\nu_{\rm {tur}} \doteq \langle V_{\rm {tur}}\rangle \langle \... ...{\rm {tur}} V_{\rm {S}} \times \alpha_{\rm {2}} R_{\rm {NS}}, \end{displaymath} (17)

where $\langle V_{\rm {tur}}\rangle,~~ \langle \ell_{\rm {tur}}\rangle$ correspond to mean values of velocity and length scale of eddies in a turbulent medium, respectively. These are set to be respectively smaller than the sound speed $V_{\rm {S}}$ and smaller than the radius of the NS. Thus, $\alpha_{\rm {tur}},~\alpha_{\rm {2}}$ are constants that are set to be smaller than unity. In the present paper, all model calculations assume $\langle \ell_{\rm {tur}}\rangle = 0.1 R_{\rm {NS}}$, i.e., $\alpha_{\rm {2}}= 0.1$, whereas the parameter $\alpha _{\rm {tur}}$ may differ from one model calculation to another.

3 Heat bubble calculations

3.1 Numerical solution method

The set of hydrodynamical equations are solved using a pre-conditioned defect-correction iteration procedure. The matrix equation to be solved in each iteration is:

\begin{displaymath}{\tilde{A}}\mu = d, \end{displaymath} (18)

where $\tilde{A}$ is a preconditioner, i.e., a coefficient matrix that is similar to the Jacobian J, but which is much easier to invert.

J is obtained by calculating the entries resulting from $\partial R/\partial q,$ where R denotes the vector of equations (Eq. (4)) and q the vector of variables.

In this formulation $\mu = q^{i+1}-q^{i}$ corresponds to the correction between two successive iterations and d is the defect (for a detailed description of the method see Hujeirat et al. 2008; Hujeirat 2005).

We mention that the solver employed here relies on the conservative formulation of the hydrodynamical equations, using the finite volume formulation. For strongly time-dependent simulations, an advection scheme of third order spatial and second order temporal accuracies is used (Hujeirat 2005). As a pre-conditioner we use the approximate factorization method (AFM), which is proven to be most appropriate for modeling low and high Mach number flows (Hujeirat et al. 2007a). Low Mach number flows in astrophysics are generally encountered in the interior of stars. As these flows are gravitationally bound and pressure supported, their motions are extremely subsonic, hence their corresponding Mach number is rather low. On the other hand, the velocities in high Mach number flows are supersonic and they may easily turn into shock-dominated flows.

3.2 Heat bubble propagation in the atmosphere of a non-rotating white dwarf

Rising bubbles in stellar environments of white dwarfs have been extensively studied by Almgren et al. (2006, see also the references therein). Although white dwarfs are degenerate stars, their typical radii are approximately three orders of magnitude larger than their corresponding horizons, so that the dynamics of their atmosphere can be safely treated within the Newtonian regime.

To test the capability of our general relativistic solver at capturing the propagation of strongly time-dependent heat bubbles in strong gravitational fields, such as X-ray bursts on neutron stars, we adopt the same setup of the heat bubble problem as described in Sect. 4.3 of Almgren et al. (2006). Although we use spherical geometry and solve the equations using a general relativistic formulation, we expect our solver to provide solutions that agree well with the results of Almgren et al. (2006), that have been obtained using Cartesian coordinates, uniform grid distribution and a Newtonian solver.

The domain of calculation is restricted to the first quadrant:

\begin{displaymath}\begin{array}{lll} \mathcal{D} &=& [R_{\rm {in}}\leq r \leq R... ...leq r \leq 1.35 ]\times[0 \leq \theta \leq \pi/2], \end{array}\end{displaymath}

where length scales are measured in units of the radius of the star's core.

The domain $\mathcal{D}$ is divided into non-uniformly distributed finite volume cells: 100 in the radial and 170 in the polar direction, where the minimum grid spacings is set to coincide with the center of the initial heat bubble.

As initial conditions, a constant density and temperature are assumed. As such a configuration is dynamically unstable, the flow is set to evolve hydrodynamically until a hydrostatic equilibrium is recovered. Such a strategy for constructing initial conditions is recommended when rotation, heat diffusion or a complicated equation of state are used. In the present paper however, we neglect chemical composition and use the equation of state which corresponds to an ideal gas.

A small region is then chosen (Fig. 1), where the matter is replaced by a thinner but much hotter plasma, while keeping it in pressure equilibrium with the surrounding media.

 \begin{figure} \par\includegraphics[width=5.5cm,clip]{0793Fig1.ps} \par\end{figure} Figure 1:

The starting distribution of the temperature in the domain of calculation. Low (high) values correspond to blue (yellow) color. The red color corresponds to the hot bubble. The radius is in $R_{\rm NS}$ units (see Eq. (19)).
Open with DEXTER

In Fig. 2 we show several snapshots of the rising bubble. The location in time, internal structures and surface morphology of the bubble during its rise to the surface of the white dwarfs are quite similar and agree well with those obtained using the unsplit and low Mach number schemes reported by Almgren et al. (2006).

 \begin{figure} \par\includegraphics*[width=9cm,clip]{0793Fig2.ps} \end{figure} Figure 2:

A rising bubble in the white dwarf atmosphere. The domain of calculation is covered by $100 \times 170$ non-uniformly distributed finite volume cells ( top left panel). The colored images are snapshots of the temperature of the rising bubble after 0.25, 0.5. 0.75 and 1.0 s after the flash. Low-to-high values of temperature correspond to blue-to-red. On the left panel, snapshot of the density distribution at 0.5, 0.75 and 1.0 s in black and white are shown. Radii are given in $5\times 10^8$ cm units.
Open with DEXTER

3.3 Rotating heat bubbles in deep gravitational fields of neutron stars

Similar to Sect. 3.2, we apply our general relativistic solver to model the rise of a rotating bubble starting from below the atmosphere of a rotating neutron star. The radius of the neutron star's core $R_{\rm NS}$ is set to be equal to the following radius,

\begin{displaymath}R_{\rm {NS}} = 3~ r_{\rm {g}}~\left(1 + \sqrt{1-\Omega^2_{\rm {NS}}}\right), \end{displaymath} (19)

which is smaller than the classical last stable orbit of a Schwarzschild black hole. The inner boundary of the domain of calculation is taken to be the radius of the core, whereas the outer boundary is located at $R_{\rm {out}} = 1.35~R_{\rm {NS}}$. The core is set to rotate at 40% the break-up velocity, whereas $\Omega $ at the outer radius is 10 times smaller. The matter in the domain is set to adjust its rotation to the boundaries through viscous interaction.

The initial distribution of the variables is constructed by solving the general relativistic Navier-Stokes equations to obtain stationary and differentially rotating flow configurations. The procedure runs as follows: for a given distribution of the angular velocity $\Omega $, quasi-stationary solutions for the density, temperature and poloidal components of the momentum equations are sought. Using these distributions as initial conditions, we then re-run the calculations to seek quasi-stationary solutions for the equations describing the time-evolution of the density, temperature, poloidal component and also for the toroidal component of the momentum equations.

As in the previous section, the heat bubble is injected into the domain by replacing the medium at a certain location by a hot and tenuous plasma, while keeping it in pressure and rotational equilibrium with the surrounding media.

To assure consistency of the numerical scheme with the physical problem and minimize possible numerical errors, we use an advection scheme of third order accuracy in space and second order in time. The grid spacing, $({\rm d}x_j,~r {\rm d}\theta_k),$ at the center of the burst is taken to be 0.001. This means that numerical errors resulting from the spatial discretization of the transport operators are of the order of 10-9. Errors may result also from the time-discretization of the equations. However, as the problem is strongly time-dependent, the time step size has been taken to be equivalent to CFL = 0.5. This implies a time step size of the order $\delta t \sim 2 \times 10^{-4}$. In this case, the maximum possible temporal error is of order 10-7. Consequently, the combined maximum numerical error resulting from the advection scheme would be of order 10-7. An additional source of numerical errors is the discretization of second order operators, i.e., viscous operators in the Navier-Stokes equations. Within the context of finite volume formulation, the maximum possible numerical errors here cannot be larger than dxj2, which is of the order of 10-6. This is three orders of magnitude smaller than the smallest turbulent viscosity coefficient used in the present investigations.

Unlike the calculations in the previous section, the purpose of the present calculations is to unveil the response of the surrounding region to violent events associated with a dramatic change in the distribution of the angular momentum.

 \begin{figure} \par\includegraphics[width=9cm,clip]{0793Fig3.ps} \par\end{figure} Figure 3:

The distribution of the angular velocity in the bursting region for four different viscosity parameters $\alpha _{\rm {tur}}$ shortly after burst events $(\approx $0.1 ms). In the left panel 30 uniformly distributed isolines of the angular frequency $\Omega $ are shown. In the right panel we display the radial profiles of $\Omega $ across the bursts. The solid line in the second-right plot corresponds to a relaxed $\Omega $-profile after 10 ms. Radii are given in $R_{\rm NS}$ units (see Eq. (19)).
Open with DEXTER

Therefore we run several numerical calculations with $\alpha_{\rm {tur}} = 1.0$, 0.1, 0.01, 0.001. The initial stationary configurations have been obtained using the corresponding value of $\alpha _{\rm {tur}}$. The results are shown in Fig. 3 and can be summarized as follows:

1.
all model calculations show a pronounced deficiency of angular momentum at the central part of the burst, accompanied by a significant increase at the boundary of the bubble. Thus, the burst leads to the formation of a dynamically unstable flow-configuration: a shell of slow rotating matter is bounded both from below and from above by relatively fast rotating matter;
2.
the viscous-induced fronts of angular momentum are found to propagate outwards at much faster speed than inwards. This implies that the matter in the deeper layers adapts its conditions to the inner boundary much faster than at the outer boundary. This is a consequence of the adopted prescription of the viscosity $\eta$, which is more effective in hotter and denser regions of the plasma. Also, the outwards-oriented velocity is obviously higher than the inward one (Fig. 4);
3.
the difference between the rotational velocity at the center of the bubble and that at its boundaries becomes more significant, the smaller the $\alpha _{\rm {tur}}$-parameter chosen.
In the $\alpha_{\rm {tur}} = 1.0$ case (see Fig.  3a1), the rise time of the bubble is found to be extremely long compared to those obtained with smaller $\alpha_{\rm {tur}}.$

 \begin{figure} \par\includegraphics[width=8.5cm,clip]{0793Fig4.ps} \par\end{figure} Figure 4:

The profile of the radial velocity across the heat bubble shortly after the burst (0.3 ms), using $\alpha _{\rm {tur}}=0.1$. Obviously, the outward-oriented fluid motion is much stronger than the inward one. The radius is in $R_{\rm NS}$ units (see Eq. (19)).
Open with DEXTER

Here the viscous pressure in the radial direction $P^r_{\rm {visc}}\sim\alpha_{\rm {tur}}\frac{\partial u^r }{\partial r}$ has an opposite sign and acts to reduce the effects of the thermal pressure  $P_{\rm {th}}$. In the extreme case, when $P_{\rm {th}} + P_{\rm {visc}} \rightarrow 0$, the effective sound speed, $V^{\rm eff}_{\rm {S}} \sim \delta(P_{\rm {the}} + P_{\rm {visc}})/\delta \rho \rightarrow 0$, hence the propagation time $\tau_{\rm {pro}}\sim R/V^{\rm eff}_{\rm {S}} \rightarrow \infty$. The same effect would have occurred if weak magnetic fields were included. Since the media is magnetic-rotation unstable, the generation and dissipation of turbulence would yield an effective viscosity coefficient comparable to $\alpha _{\rm {tur}}$. However, this effect becomes irrelevant when magnetic fields are beyond equipartition.

 \begin{figure} \par\includegraphics[width=12cm,clip]{0793Fig5.ps} \par\end{figure} Figure 5:

25 isolines of the angular velocity ( left, $\Omega_{\rm{min}}=0.03, \Omega_{\rm{\max}}=0.45$) and of the Lorentz factor ( right panel, $\Gamma_{\rm{min}}=1.1, \Gamma_{\rm{\max}}=1.43$) overlayed on their corresponding colored images after 0.5 ms are shown, using $\alpha _{\rm {tur}}=0.01$. Blue-to-red color corresponds to low-to-high values of $\Omega $ an $\Gamma $. Radii are given in $R_{\rm NS}$ units (see Eq. (19)).
Open with DEXTER

In this case turbulence will be diminished and the strong magnetic flux tubes will quickly float up to the surface, probably prior to burst events.

On the other hand, $\tau_{\rm {pro}}$ becomes of the order of one second, when using $\alpha_{\rm {tur}} \sim \mathcal{O}(10^{-1})$, which fits well into the observed duration of burst events on NSs.

In the low $\alpha _{\rm {tur}}$ cases, (see Figs. 3b1,c1,d1 and Fig. 5), the significant increase of the rotational velocity at the outer $\Omega $-fronts is obvious, but unreasonably large. In the Fig. 3d1, the matter at the outer front is found to be gravitationally unbound to the central NS, giving rise to strong outflows.

However, as outflows during X-ray bursts can be excluded on observational grounds, we conclude that $\alpha _{\rm {tur}}$ must acquire much larger values, and that specifically $\alpha_{\rm {tur}} \sim \mathcal{O}(10^{-1})$.

To study the viscous-reaction of the matter in the adjusting layers to the sudden increase of rotation induced by a burst, we have run separate calculations in the following manner. A solution for the hydrodynamical equations including rotation has been hydrodynamically calculated. As a second step, we modified the $\Omega-$profile by including a Gaussian perturbation of the form $\tau_{\rm {0}}$ depicted in Fig. 6. We then followed the time evolution of this profile on the viscous time scale. The profiles $\tau_{\rm {1}}, \tau_{\rm {2}} , \tau_{\rm {10}}$ correspond to $\tau_{\rm {visc}}/10,\tau_{\rm {visc}}/5,\tau_{\rm {visc}}$.

These calculations show that the effect of turbulent viscosity is to mainly transport angular momentum outwards. As a consequence, the deficiency in the rotational support in the deep layers enhances the compression of matter and give rise to an additional burst in the neighboring shell. This chain of reactions may run away to spread over the whole surface of the NS on the viscous time scale, which is of the order of one second, assuming $\alpha_{\rm {tur}} \approx 0.1.$

When the outer layers cool, the turbulent viscosity decreases and the corresponding viscous time scale increases as $\tau_{\rm {visc}} \sim 1/\surd{T}$. This implies that after the burst, the time scale needed for the outer layers to adjust their rotation to the bulk of the star might lengthen by an additional order of magnitude. As a consequence, the observed spin up of NSs in their cooling tail is a manifestation of the increased rotational velocity of the outer layers caused by the burst events, but which, eventually, should decrease at later times.

We note that, if magnetic fields were included, the dynamical evolution of the rotating bubbles would proceed in a similar manner as presented here, provided they are in sub-equipartition with the thermal energy. The effective time scale of X-ray burst on NSs is much shorter than the dynamical time scale and therefore is too fast for the magnetic-rotational instability (MRI; Hawley et al. 1996) to fully develop during the burst. On the other hand, a highly turbulent pre-burst media with a fully developed MRI on the verge of entering an enhanced-dissipation phase cannot be excluded as a possible mechanism that could lead to runawy ignition and subsequently to X-ray bursts.

4 Summary and conclusions

In this paper we have presented the set of general relativistic Navier-Stokes equations in axi-symmetry, using the Boyer-Lindquist coordinates in the background of a slowly rotating neutron star.

To make the numerical solution method viable, the collection of viscous operators has been reduced considerably by considering the dominant second order operators only. Such operators generally stabilize the matrix inversion procedure and enhance the convergence rate of the numerical method.

The set of equations are solved numerically using an implicit solution procedure which is based on a pre-conditioned defect-correction iterative method. Similar to Taylor-flows between concentric spheres, we use the ``Approximate Factorization Method'' as a pre-conditioner, which has superior converging properties over other non-symmetric methods, such as the black-white line Gauss-Seidel method.

In the non-rotating case, we have shown that the solver is capable of accurately reproducing the time-evolution of heat bubbles during their rise to the surface of the white dwarfs, as reported by Almgren et al. (2006) and Nonaka (2008).

In the rotating case, it has been shown that viscous-generated fronts inside heat bubbles propagate into the surrounding quite rapidly. The effect of these fronts is mainly to transport angular momentum to the outer layers, leaving the matter in the deeper layer with less rotational support, hence more compressed.

 \begin{figure} \par\includegraphics[width=6.5cm, bb=0 350 474 657,clip]{0793Fig6... ...udegraphics[width=6.5cm, bb=277 512 522 737,clip]{0793Fig7.ps} \par\end{figure} Figure 6:

Successive snapshots of the omega-profiles in the outer layers of a neutron star. Shell I shows the burst-induced propagation of omega-fronts into the neighboring shells. In shell II, the viscous-induced reaction to the omega-fronts is shown. The schematic picture above shows that viscosity transports the excess of angular momentum preferably in the outward direction. The lower figure, which is obtained using hydrodynamical calculations, clearly confirms this behavior. The radius is in $R_{\rm NS}$ units (see Eq. (19)).
Open with DEXTER

Our numerical results show that a viscosity parameter of the order of $\alpha_{\rm {tur}}\sim 0.1$ is the most reliable value for fitting observations of NSs during X-ray bursts. A much larger value yields a propagation time that is much longer than one second, whereas smaller values yield unstable shell configurations and give rise to gravitationally unbound outflows.

In addition, a possible mechanism that may underly the rapid spread of burning ignition fronts has been presented. Accordingly, the viscous-generated fronts inside a heat bubble may transport angular momentum into the adjusting layers in the polar direction. The viscosity then acts to decouple matter from angular momentum, subsequently enhancing compression in the deeper layers and giving rise to thermonuclear runaway. This chain of reactions may run away to spread over the whole surface of the NSs on the viscous time scale $\tau_{\rm {vis}}\sim 1~ \rm {s}$ .

We note that for the process of decoupling of angular momentum from matter to operate efficiently, the matter must be viscous and highly stratified. This is contrary to the case when the flow is ideal and the angular momentum is frozen-in to the matter. Here an excess of angular momentum would, of course, shed the matter away. In the non-ideal case, viscous torques acts to spread angular momentum in the opposite direction to gravity, which in turn attracts the matter and tries to sink it deeper into the gravitational well of the central object. As gravity in gravitationally bound flows is the dominant force, the viscosity coefficient $\alpha _{\rm {tur}}$ must be sufficiently large for the decoupling process to operate efficiently.

Apparently, as our calculations show, $\alpha_{\rm {tur}}\sim 0.1$ is not only the preferable value in accretion disks, but in the outer pre-burst turbulent layers of accreting NSs also.

Finally, the increased rotational velocity of the outer layers may be connected to the observed spin up of neutron stars during X-ray bursts. However, this increase will eventually be followed by a spin down in later times when the outer layers cool down to their pre-burst thermal state.

Concerning frame dragging, while included in the present simulations, its effects are still to be quantified. Also, the effects of magnetic fields and thermal diffusion is the subject of an ongoing work.

Acknowledgements
A.H. thanks Bernhard Keil and Paul P. Hilscher for carefully reading the manuscript and for the valuable discussions. This work is supported by the Klaus-Tschira Stiftung under the project number 00.099.2006.

References

Footnotes

... derivatives[*]
i.e., $\partial_r(f\partial_\theta q)$ and $\partial_\theta(f\partial_r q)$, where f denotes an arbitrary non-linear viscous function.

All Figures

  \begin{figure} \par\includegraphics[width=5.5cm,clip]{0793Fig1.ps} \par\end{figure} Figure 1:

The starting distribution of the temperature in the domain of calculation. Low (high) values correspond to blue (yellow) color. The red color corresponds to the hot bubble. The radius is in $R_{\rm NS}$ units (see Eq. (19)).
Open with DEXTER
In the text

  \begin{figure} \par\includegraphics*[width=9cm,clip]{0793Fig2.ps} \end{figure} Figure 2:

A rising bubble in the white dwarf atmosphere. The domain of calculation is covered by $100 \times 170$ non-uniformly distributed finite volume cells ( top left panel). The colored images are snapshots of the temperature of the rising bubble after 0.25, 0.5. 0.75 and 1.0 s after the flash. Low-to-high values of temperature correspond to blue-to-red. On the left panel, snapshot of the density distribution at 0.5, 0.75 and 1.0 s in black and white are shown. Radii are given in $5\times 10^8$ cm units.
Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=9cm,clip]{0793Fig3.ps} \par\end{figure} Figure 3:

The distribution of the angular velocity in the bursting region for four different viscosity parameters $\alpha _{\rm {tur}}$ shortly after burst events $(\approx $0.1 ms). In the left panel 30 uniformly distributed isolines of the angular frequency $\Omega $ are shown. In the right panel we display the radial profiles of $\Omega $ across the bursts. The solid line in the second-right plot corresponds to a relaxed $\Omega $-profile after 10 ms. Radii are given in $R_{\rm NS}$ units (see Eq. (19)).
Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=8.5cm,clip]{0793Fig4.ps} \par\end{figure} Figure 4:

The profile of the radial velocity across the heat bubble shortly after the burst (0.3 ms), using $\alpha _{\rm {tur}}=0.1$. Obviously, the outward-oriented fluid motion is much stronger than the inward one. The radius is in $R_{\rm NS}$ units (see Eq. (19)).
Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=12cm,clip]{0793Fig5.ps} \par\end{figure} Figure 5:

25 isolines of the angular velocity ( left, $\Omega_{\rm{min}}=0.03, \Omega_{\rm{\max}}=0.45$) and of the Lorentz factor ( right panel, $\Gamma_{\rm{min}}=1.1, \Gamma_{\rm{\max}}=1.43$) overlayed on their corresponding colored images after 0.5 ms are shown, using $\alpha _{\rm {tur}}=0.01$. Blue-to-red color corresponds to low-to-high values of $\Omega $ an $\Gamma $. Radii are given in $R_{\rm NS}$ units (see Eq. (19)).
Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=6.5cm, bb=0 350 474 657,clip]{0793Fig6... ...udegraphics[width=6.5cm, bb=277 512 522 737,clip]{0793Fig7.ps} \par\end{figure} Figure 6:

Successive snapshots of the omega-profiles in the outer layers of a neutron star. Shell I shows the burst-induced propagation of omega-fronts into the neighboring shells. In shell II, the viscous-induced reaction to the omega-fronts is shown. The schematic picture above shows that viscosity transports the excess of angular momentum preferably in the outward direction. The lower figure, which is obtained using hydrodynamical calculations, clearly confirms this behavior. The radius is in $R_{\rm NS}$ units (see Eq. (19)).
Open with DEXTER
In the text

Copyright ESO 2009

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.