Turbulent Mixing in the Surface Layers of Accreting Neutron Stars

TL;DR

This paper investigates how turbulent mixing driven by magnetic dynamo effects in accreting neutron stars influences X-ray burst behavior, revealing new regimes of ignition and stable burning that depend on accretion rate and spin.

Contribution

It introduces models showing how magnetic dynamo-induced turbulence affects mixing and ignition regimes in neutron star surface layers, explaining observed burst variability.

Findings

Mixing can trigger short recurrence time bursts (~5-30 min).

Strong mixing leads to stable helium burning, quenching X-ray bursts.

Different accretion rates and spin frequencies influence mixing and burst regimes.

Abstract

During accretion a neutron star (NS) is spun up as angular momentum is transported through its surface layers. We study the resulting differentially rotating profile, focusing on the impact this has for type I X-ray bursts. The predominant viscosity is likely provided by the Tayler-Spruit dynamo. The radial and azimuthal magnetic field components have strengths of ~10^5 G and ~10^10 G, respectively. This leads to nearly uniform rotation at the depths of interest for X-ray bursts. A remaining small shear transmits the accreted angular momentum inward to the NS interior. Though this shear gives little viscous heating, it can trigger turbulent mixing. Detailed simulations will be required to fully understand the consequences of mixing, but our models illustrate some general features. Mixing has the greatest impact when the buoyancy at the compositional discontinuity between accreted matter and ashes is overcome. This occurs at high accretion rates, at low spin frequencies, or may depend on the ashes from the previous burst. We then find two new regimes of burning. The first is ignition in a layer containing a mixture of heavier elements from the ashes. If ignition occurs at the base of the mixed layer, recurrence times as short as ~5-30 minutes are possible. This may explain the short recurrence time of some bursts, but incomplete burning is still needed to explain these bursts' energetics. When mixing is sufficiently strong, a second regime is found where accreted helium mixes deep enough to burn stably, quenching X-ray bursts. We speculate that the observed change in X-ray burst properties near one-tenth the Eddington accretion rate is from this mechanism. The carbon-rich material produced by stable helium burning would be important for triggering and fueling superbursts. (abridged)

Figures

Equations

$$\frac{d}{dt}(R^2\Omega )=\frac{1}{R^2\rho }\frac{\partial }{\partial z}\left(\rho \nu R^4\frac{\partial \Omega }{\partial z}\right),$$

$$-\frac{\dot{M}}{4\pi R^2}\frac{d}{dz}(R^2\Omega )=\frac{1}{R^2}\frac{d}{dz}\left(\rho \nu R^4\frac{d\Omega }{dz}\right).$$

$$-\dot{M}{R}^2{\Omega }_{\mathrm{K}}+\dot{M}{R}^2\Omega =-4\pi \rho \nu R^{4}d\Omega /dz.$$

$$4\pi R^3\rho \nu q\Omega =\dot{M}{R}^2{\Omega }_{\mathrm{K}},$$

$$Ri\equiv \frac{N^2}{q^2{\Omega }^2},$$

$$N^2=N_T^2+N_{\mu }^2.$$

$$N_T^2=\frac{g}{H}\frac{{\chi }_T}{{\chi }_{\rho }}\left[{\nabla }_{\mathrm{ad}}-{\left(\frac{d\ln T}{d\ln P}\right)}_{\ast }\right],$$

$$H=P/\rho g=33.1\mathrm{cm}{\mu }_{1.33}^{-1}{T}_8,$$

$$N_{\mu }^2=-\frac{g}{H{\chi }_{\rho }}\left[{\chi }_{{\mu }_e}{\left(\frac{d\ln {\mu }_e}{d\ln P}\right)}_{\ast }+{\chi }_{{\mu }_i}{\left(\frac{d\ln {\mu }_i}{d\ln P}\right)}_{\ast }\right],$$

$$K>{\nu }_{k}R{e}_c.$$

$${\nu }_k=1.4\times 10^{-3}\mathrm{c}{\mathrm{m}}^2{\mathrm{s}}^{-1}{\rho }_6^{-1}{T}_8^{5/2},$$

$$K=48.8\mathrm{c}{\mathrm{m}}^2{\mathrm{s}}^{-1}{\mu }_{1.33}{\kappa }_{0.04}^{-1}{\rho }_6^{-2}{T}_8^3,$$

$$R{i}_s\equiv \frac{{\nu }_{k}R{e}_c}{K}\frac{N_T^2}{q^2{\Omega }^2}.$$

$$N_T={\left(\frac{3}{20}\frac{g}{H}\right)}^{1/2}=9.2\times 10^5{\mathrm{s}}^{-1}{\mu }_{1.33}^{1/2}{T}_8^{-1/2}.$$

$$q_{\mathrm{KH}}=223{\kappa }_{0.04}^{1/2}{\rho }_6^{1/2}{T}_8^{-3/4}{\Omega }_{0.1}^{-1},$$

$$R{i}_{\mathrm{BC}}\equiv 4{\left(\frac{R}{H}\right)}^2{\left(\frac{\Omega }{N}\right)}^2=8.5\times 10^3{\mu }_{1.33}{T}_8^{-1}{\Omega }_{0.1}^2.$$

$$\displaystyle\nu_{\rm BC}=\left\{\begin{array}[]{ccc}\displaystyle\frac{1}{3}\frac{1}{Ri^{1/2}}H^{2}\Omega,&&Ri\leq Ri_{\rm BC},\\ &&\\ \displaystyle\frac{1}{3}\frac{Ri_{\rm BC}}{Ri^{3/2}}H^{2}\Omega,&&Ri>Ri_{\rm BC},\end{array}\right.$$

$${\dot{m}}_{\mathrm{Edd}}$$

$$\displaystyle Ri=\left\{\begin{array}[]{lcc}2.2\times 10^{3}\ \mu_{1.33}^{-3/2}\rho_{6}T_{8}^{3/2}\dot{m}_{0.1}^{-1}\Omega_{0.1},&&Ri\leq Ri_{\rm BC}\\ &&\\ 4.4\times 10^{3}\ \mu_{1.33}^{-1/4}\rho_{6}^{1/2}T_{8}^{1/4}\dot{m}_{0.1}^{-1/2}\Omega_{0.1}^{3/2},&&Ri>Ri_{\rm BC}\end{array}\right.$$

$$\displaystyle q_{\rm BC}=\left\{\begin{array}[]{lcc}14\ \mu_{1.33}^{5/4}\rho_{6}^{-1/2}T_{8}^{-5/4}\dot{m}_{0.1}^{1/2}\Omega_{0.1}^{-3/2},&&Ri\leq Ri_{\rm BC}\\ &&\\ 9.9\ \mu_{1.33}^{5/8}\rho_{6}^{-1/4}T_{8}^{-5/8}\dot{m}_{0.1}^{1/4}\Omega_{0.1}^{-7/4},&&Ri>Ri_{\rm BC}\end{array}\right.$$

$$\frac{B_r{B}_{\varphi }}{4\pi }=\rho {\nu }_{e}q\Omega .$$

$$B_{r,\mathrm{crit}}=(4\pi \rho {)}^{1/2}\frac{H_{\Omega }}{t}=(4\pi \rho {)}^{1/2}\frac{R}{qt},$$

$$t_{\Omega }={\left(\frac{3R^2{\pi }^2}{\eta {\Omega }^2{q}^2}\right)}^{1/3},$$

$$\eta \approx 0.7\mathrm{c}{\mathrm{m}}^2{\mathrm{s}}^{-1}{T}_8^{-3/2}.$$

$$B_{\Omega }$$

$$l_z<R{\omega }_{\mathrm{A}}/N,$$

$$l_z^2>\eta /{\sigma }_B=\eta \Omega /{\omega }_{\mathrm{A}}^2.$$

$${\left(\frac{{\omega }_{\mathrm{A}}}{\Omega }\right)}^4>\frac{\eta }{R^2\Omega }{\left(\frac{N}{\Omega }\right)}^2.$$

$$B_{\varphi }={\sigma }_B^{-1}q\Omega B_r.$$

$$q_{\min }={\left(\frac{N}{\Omega }\right)}^{7/4}{\left(\frac{\eta }{R^{2}N}\right)}^{1/4}$$

Full text

††slugcomment: Accepted for publication in The Astrophysical Journal

1. Introduction

As neutron stars (NSs) in low mass X-ray binaries accrete material from their companions they are expected to be spun up by this addition of angular momentum, possibly becoming millisecond pulsars (Bhattacharya & van den Heuvel, 1991). This suspicion has received support by the discovery of accretion driven millisecond pulsars (Wijnands & van der Klis, 1998), as well as the millisecond oscillations seen during type I X-ray bursts (Chakrabarty et al., 2003), the unstable ignition of the accumulating fuel (for reviews, see Lewin et al., 1995; Strohmayer & Bildsten, 2006; Galloway et al., 2006). The majority of the rotational kinetic energy of the accreted material is dissipated at low densities ($\sim 10^{-4}-10^{-1}\ {\rm g\ cm^{-3}}$) in the boundary layer (for example, Inogamov & Sunyaev, 1999). Nevertheless, angular momentum must be transported into the NS interior if it is to be spun up, even if at times the angular momentum could be radiated as gravitational waves (Bildsten, 1998b). This transport implies a non-zero, albeit small, shear throughout the outer liquid parts of the NS. Such shearing may lead to viscous heating as well as chemical mixing at depths far below the low density boundary layer where the majority of the shearing occurs.

Such a picture was previously investigated by Fujimoto (1988, 1993) for both accreting NSs and white dwarfs. His studies demonstrated that the shearing profile of a NS could be large enough that mixing through baroclinic instabilities may be important at the depths of type I X-ray bursts. Such a result is compelling because it may help to explain some of the remaining discrepancies between the theory and observations of these bursts. Though a simple limit cycle picture has been successful in qualitatively explaining the bursts’ primary characteristics, including their energies ($10^{39}-10^{40}\ {\rm ergs}$), recurrence times (hours to days), and durations ($\sim 10-100\ {\rm s}$) (Fujimoto et al., 1981; Bildsten, 1998a), outstanding problems still remain (Fujimoto et al., 1987; van Paradijs et al., 1988; Bildsten, 2000). Chief among these is the critical accretion rate required for stable accumulation (which we discuss in more detail in §6). Another problem is the occurrence of successive X-ray bursts with recurrence times as short as $\sim 10$ minutes (Galloway et al., 2006, and references therein), for which there has been speculation that this could be due to mixing during accretion (Fujimoto et al., 1987).

In this present work we reassess the importance of angular momentum transport and the resulting mixing. We find that the hydrodynamic instabilities studied by Fujimoto (1993) are insufficient to prevent strong shearing of magnetic fields. This leads to generation of the Tayler-Spruit dynamo (Spruit, 2002), where toroidal field growth is balanced by Tayler instabilities to create a steady-state magnetic field, which provides a torque on the shearing layers that is larger than any purely hydrodynamic mechanism. The result is nearly uniform rotation and little viscous heating, too little to affect either X-ray bursts or superbursts (thermonuclear ignition of ashes from previous X-ray bursts, Cumming & Bildsten, 2001; Strohmayer & Brown, 2002). Turbulent mixing is found to be non-negligible and in some cases may mix fresh material with the ashes of previous bursts. The key is whether the strong buoyancy due to the larger density of the ashes below can be overcome. When mixing occurs, we find two new regimes of burning for an accreting NS. The first is that fuel can be mixed down to depths necessary for premature unstable ignition. The timescale for ignition of such bursts is short enough ($\sim 5-30\ {\rm minutes}$) to explain the short recurrence time bursts (Galloway et al., 2006). The second is that if material is mixed to sufficiently large depths (and therefore temperatures) it can burn stably, ceasing X-ray bursts altogether. Such a mechanism may be responsible for the quenching of X-ray bursts seen at surprisingly low accretion rates for many atoll class NSs (Cornelisse et al., 2003). We derive an analytic formula that estimates what spins and accretion rates necessary to transfer between these two regimes (eq. [87]).

1.1. The Basics of Angular Momentum Transport

Before beginning more detailed analysis, it is useful to present some general equations for angular momentum transport. This demonstrates why we expect the NS interior to be shearing and explicitly relates the torque from the accreting matter to this shear. Figure 1 schematically shows the expected rotation profile from the accretion disk, through the boundary layer, and into the NS interior. Material accreted at a rate $\dot{M}$ reaches the NS surface with a nearly Keplerian spin frequency of $\Omega_{\rm K}=(GM/R^{3})^{1/2}=1.4\times 10^{4}\ {\rm s^{-1}}\ M_{1.4}^{1/2}R_{6}^{-3/2}$, where $M_{1.4}\equiv M/1.4M_{\odot}$ and $R_{6}\equiv R/10^{6}\ {\rm cm}$, which has a kinetic energy per nucleon of $\approx 200\ {\rm MeV\ nuc^{-1}}$. The majority of this energy is dissipated in a boundary layer of thickness $H_{\rm BL}\ll R$ (as studied by Inogamov & Sunyaev, 1999) and never reaches far into the NS surface. Nevertheless, angular momentum is being being added at a rate of $\dot{M}(GMR)^{1/2}$, so that a torque of this magnitude must be communicated into the NS. This implies a non-zero shear rate in the interior liquid layers, down to the solid crust. In the present work we are interested in the shear at the depths where X-ray bursts ignite, near $\rho\approx 10^{6}\ {\rm g\ cm^{-3}}$, which is well below the boundary layer. Though we focus on a magnetic mechanism for angular momentum transport, we expect the boundary layer and densities up to $\approx 5\times 10^{3}\ {\rm g\ cm^{-3}}$ to be dominated by hydrodynamic instabilities (which we discuss in more detail in §3).

The pressure scale height at the depth of X-ray bursts is $H\approx 30\ {\rm cm}\ll R$, which allows us to use a Cartesian coordinate system with $z$ as the radial coordinate. This is far below the depths considered by Inogamov & Sunyaev (1999) when they investigated the uneven covering of the NS surface by accreted fuel. Furthermore, all of the transport mechanisms we consider work most efficiently in directions perpendicular to gravity (because no work is performed), thus it is a good approximation to consider the surface layers as concentric spheres, each with constant $\Omega$. Transfer of angular momentum is reduced to a one-dimensional diffusion equation (Fujimoto, 1993)

[TABLE]

where $\Omega$ is the NS spin frequency and $\nu$ is the viscosity. The total time derivative is given by $d/dt=\partial/\partial t+V_{\rm adv}\partial/\partial z$, where $V_{\rm adv}$ is the advecting velocity of the fluid in an Eulerian frame. In steady-state, we take $\partial/\partial t=0$ and set $V_{\rm adv}=-\dot{M}/(4\pi R^{2}\rho)$ so that

[TABLE]

Integrating from the surface where the spin is $\Omega_{\rm K}$ (due to the disk) down to a depth $z$ where the local spin is $\Omega$, and assuming $d\Omega/dz=0$ at the surface,

[TABLE]

Taking the limit $\Omega\ll\Omega_{\rm K}$ we find,

[TABLE]

where $q\equiv d\log\Omega/d\log z$ is the shear rate. This equation shows that $q>0$ when angular momentum is transported inward. In general, we will find that $q$ is rather small ($\lesssim 1$) at the depths of interest. Nevertheless $q$ is large enough to activate instabilities that help to transport angular momentum, as well as mix material.

1.2. Outline of Paper

We begin by comparing and contrasting some well-known hydrodynamic instabilities in §2 and estimate the resulting shear rates. In §3, we discuss the consequences that this shearing has on a magnetic field, which motivates implementation of the Tayler-Spruit dynamo. In §4, we calculate accumulating NS envelopes without the effects of viscous angular momentum transport. This enables us to judge when such effects must be incorporated. We calculate models including mixing in §5. We conclude in §6 with a summary of our results and a discussion of type I X-ray burst and superburst observations.

2. Hydrodynamic Viscosity Mechanisms

In the following sections we discuss hydrodynamic instabilities. This is not an exhaustive survey (for further details, see Heger et al., 2000), but rather is meant to summarize those instabilities that are most crucial to our problem, so as to set the context for the magnetic transport mechanism we study later.

2.1. Kelvin-Helmholtz Instability

The Kelvin-Helmholtz instability (also referred to as the dynamical shear instability) is governed by the Richardson number

[TABLE]

where $N$ is the Brunt-Väisälä frequency, which is composed of contributions from thermal and compositional buoyancy,

[TABLE]

The thermal buoyancy is given by

[TABLE]

where $g=GM/R^{2}=1.87\times 10^{14}\ {\rm cm\ s^{-2}}M_{1.4}R_{6}^{-2}$ is the surface gravitational acceleration (ignoring the effects of general relativity), $\chi_{Q}=\partial\ln P/\partial\ln Q$, with all other intensive variables set constant, $\nabla_{\rm ad}=(\partial\ln T/\partial\ln P)_{\rm ad}$ is the adiabatic temperature gradient, the asterisk refers to derivatives of the envelope’s profile, and the pressure scale height is,

[TABLE]

where $T_{8}\equiv T/10^{8}\ {\rm K}$. For these analytic estimates we assume an ideal gas equation of state and use a pure helium composition with a mean molecular weight of $\mu=1.33\mu_{1.33}$. We omit the scalings with mass and radius to simplify presentation. The compositional buoyancy is (Bildsten & Cumming, 1998)

[TABLE]

where $\mu_{e}$ ($\mu_{i}$) is the mean molecular weight per electron (ion). Linear analysis shows that Kelvin-Helmholtz instability occurs when $Ri<1/4$, which develops into strong turbulence that readily transports angular momentum. This result assumes that thermal diffusion can be ignored for the unstable fluid perturbations, in other words, that the perturbations are adiabatic.

Fluid perturbations with a characteristic size $L$ and speed $V$ become non-adiabatic when the timescale for thermal diffusion, $L^{2}/K$, where $K$ is the thermal diffusivity, is less than the timescale of the perturbation, $L/V$. The ratio of these two timescales is the Péclet number, $Pe\equiv VL/K$ (Townsend, 1958). The restoring force provided by thermal buoyancy is weakened when $Pe<1$, which requires the substitution of $N_{T}^{2}\rightarrow PeN^{2}_{T}$ and promotes instability for regions where $N_{T}>N_{\mu}$. Thermal diffusion is most efficient at small lengthscales, which motivates setting $LV/\nu_{k}=Re_{c}$ (Zahn, 1992), where $\nu_{k}$ is the kinematic viscosity and $Re_{c}$ is the critical Reynolds number for turbulence, which is of order 1000. This gives the Péclet number approximately related to the Prandtl number, $Pr$, by $Pe\approx Re_{c}Pr$. The turbulent perturbations are thus non-adiabatic when (Zahn, 1992)

[TABLE]

In the non-degenerate surface layers the kinematic viscosity is dominated by ions, and has a value of (Spitzer, 1962)

[TABLE]

where $\rho_{6}\equiv\rho/10^{6}\ {\rm g\ cm^{-3}}$, and we assume a Coulomb logarithm of $\ln\Lambda=20$. Setting $K=16\sigma_{\rm SB}T^{3}/(3c_{p}\kappa\rho^{2})$, where $\sigma_{\rm SB}$ the Stefan-Boltzmann constant, $c_{p}$ the specific heat, and $\kappa$ the opacity, the thermal diffusivity is

[TABLE]

where we approximate $c_{p}=5k_{\rm B}/2\mu m_{p}$ and scale the opacity to $\kappa_{0.04}\equiv\kappa/0.04\ {\rm cm^{2}\ g^{-1}}$ (the opacity is largely given by electron scattering, but is decreased by degeneracy effects, see Paczyński, 1983; Bildsten, 1998a). Substituting equations (11) and (12) into equation (10), we find that the perturbations are non-adiabatic at depths of \rho~{}\lesssim~{}\mbox{3\times 10^{7}}\ {\rm g\ cm^{-3}}\ T_{8}^{1/2}. The new “secular” Richardson number associated with this limit is,

[TABLE]

When $Ri_{s}<1/4$, the so-called “secular shear instability” arises.

The competing effects of accretion increasing $q$ versus turbulence developing when $Ri_{s}<1/4$ (and decreasing $q$) drive the surface layers toward marginally satisfying $Ri_{s}=1/4$ (assuming at this moment that the sole viscous mechanism is the Kelvin-Helmholtz instability). This expectation is borne out in the white dwarf studies of Yoon & Langer (2004). Thus we can trivially estimate the $q$ due to this mechanism. The thermal buoyancy is (Bildsten, 1998a)

[TABLE]

We substitute $Ri_{s}=1/4$ into equation (13), and assuming $Re_{c}=1000$, solve for a shear rate of

[TABLE]

where $\Omega_{0.1}=\Omega/0.1\Omega_{\rm K}$. A shear rate this large would promote prodigious viscous heating, as well as ample mixing, but as we soon show, such a large shear is prohibited by other instabilities.

2.2. Baroclinic Instability

Another important hydrodynamic instability that has been studied extensively for accreting degenerate stars is the baroclinic instability (Fujimoto, 1993), which we quickly summarize here. The interested reader should consult Fujimoto (1987, 1988) for further details (also see the discussion in Cumming & Bildsten, 2000).

The baroclinic instability arises because surfaces of constant pressure and density no longer coincide if hydrostatic balance is to be maintained when differential rotation is present. In such a configuration, fluid perturbations along nearly horizontal directions are unstable, though with a sufficient radial component to allow mixing of angular momentum and material. The instability can roughly be broken into two limits, depending on a critical baroclinic Richardson number (Fujimoto, 1987),

[TABLE]

When $Ri>Ri_{\rm BC}$, Coriolis effects limit the horizontal scale of perturbations. This results in two parametrizations for viscosity estimated from linear theory (Fujimoto, 1993),

[TABLE]

where a factor of order unity is usually included in these prescriptions, called $\alpha_{\rm BC}$, to account for uncertainty in how linear theory relates to the saturated amplitudes of the instability. For simplicity, we set $\alpha_{\rm BC}=1$ in our analysis below.

By substituting $\nu_{\rm BC}$ into the angular momentum equation (eq. [4]), we solve for the shearing profile. Since we are interested in how the shear relates to the bursting properties, it useful to write these results in terms of the local accretion rate $\dot{m}=\dot{M}/(4\pi R^{2})$, which is typically parametrized in terms of the local Eddington rate

[TABLE]

where $X$ is the hydrogen mass fraction and $\sigma_{\rm Th}$ is the Thomson cross-section. The Richardson number in each case is then

[TABLE]

where $\dot{m}_{0.1}=\dot{m}/0.1\dot{m}_{\rm Edd}$ and we have used the hydrogen deficient ($X=0$) value for $\dot{m}_{\rm Edd}$ (eq. [21]). The transition to the case $Ri\gtrsim Ri_{\rm BC}$ occurs roughly at depths $\rho\gtrsim 4\times 10^{6}\ {\rm g\ cm^{-3}}\ T_{8}^{-5/2}$. The shear rate for each case is

[TABLE]

This demonstrates that generally $q_{\rm BC}\ll q_{\rm KH}$, so that the baroclinic instability triggers before the Kelvin-Helmholtz instability. This prevents the shear rate from ever becoming large enough for the Kelvin-Helmholtz instability to operate at depths of $\rho\gtrsim 6\times 10^{4}\ {\rm g\ cm^{-3}}$.

2.3. Other Hydrodynamic Instabilities

In addition to the two instabilities described above, there are a number of other possibilities including, but not limited to, Eddington-Sweet circulation (von Zeipel, 1924a, b; Baker & Kippenhahn, 1959), Solberg-Høiland instability (Wasiutyński, 1946; Tassoul, 1978; Endal & Sofia, 1978), Goldreich-Schubert-Fricke instability (Goldreich & Schubert, 1967; Fricke, 1968) and Ekman pumping (Pedlosky, 1987). At this time we avoid assessing each of these individually. As we show below, magnetic fields are likely to dominate and are interesting since they have received the least attention in past works.

3. The Importance of Magnetic Effects

Given the larger than order unity shear rates derived above, we estimate the consequences this has for a magnetic field. The point we want to emphasize is that even a reasonably small field will be wrapped by the shear flow until it becomes dynamically important.

Assuming shellular rotation, a component of radial field is stretched to have a toroidal component, $B_{\phi}=nB_{r}$, where $n$ is the number of differential turns, given as $n=q\Omega t$, and $t$ is the duration of the shearing. The toroidal field growth is very fast. For the Kelvin-Helmholtz case, $B_{\phi}\sim B_{r}$ in merely $\sim 10^{-6}\ {\rm s}$ ($\sim 10^{-4}\ {\rm s}$ for the baroclinic case). As the toroidal field becomes larger it exerts an azimuthal stress on the shearing layer equal to $B_{r}B_{\phi}/4\pi$, which can be written as an effective viscosity, $\nu_{e}$,

[TABLE]

In a timescale $t\approx H_{\Omega}^{2}/\nu_{e}$ this torquing significantly decreases the spin of the layer, where $H_{\Omega}=dz/d\ln\Omega=R/q$ is the shearing scale height. Setting $B_{\phi}=q\Omega tB_{r}$, we solve for the critical initial radial field needed to affect the shearing,

[TABLE]

so that $t$ is basically the Alfvén travel time through a shearing scale height. Using $\rho\approx 10^{6}\ {\rm g\ cm^{-3}}$ as burning density and $t\approx 1\ {\rm hr}$, a fiducial timescale for accumulation, our estimates for $q_{\rm KH}$ and $q_{\rm BC}$ imply $B_{r,\rm crit}\sim 10^{4}\ {\rm G}$ and $\sim 10^{5}\ {\rm G}$, respectively Any initial field larger than this gets so wound up that Lorentz forces alter the shearing profile.

It is possible that the intrinsic magnetic field of the NS may be large enough that magnetic stresses never allow the shear rates to become so large. Spruit (1999) argued that this depends on the “rotational smoothing time,” the timescale for non-axisymmetric components of the magnetic field to be expelled as differential rotational brings field with opposite polarities together. This timescale is estimated as

[TABLE]

where $\eta$ is the magnetic diffusivity, which in the non-degenerate limit is given by (Spitzer, 1962)

[TABLE]

If $t_{\Omega}$ is larger than the Alfvén travel time through the layer (given by $t$ in eq. [32]), then there is sufficient time for magnetic torques to act and the shearing is merely a perturbation on the magnetic field that is quickly damped away. On the other hand, if $t_{\Omega}$ is sufficiently small then the shear dominates and the magnetic field is made axisymmetric on a timescale $t_{\Omega}$. Substituting $t_{\Omega}$ into equation (32), we estimate the critical field, below which shear dominates,

[TABLE]

The lack of persistent pulsations from accreting NSs implies a dipole field strength $\lesssim 5\times 10^{7}\ {\rm G}$ (Piro & Bildsten, 2005). We therefore consider it plausible that they may have intrinsic magnetic fields $<B_{\Omega}$. In this case, we expect a very axisymmetric magnetic field to be created in a timescale of $t_{\Omega}\sim 200\ {\rm s}$, which is then wrapped until it becomes dynamically important.

3.1. The Tayler-Spruit Dynamo

As the toroidal magnetic field continues to wrap, it becomes increasingly important to the dynamics of the shearing and is also subject to magnetohydrodynamic instabilities. The combination of these effects have been shown to give rise to the Tayler-Spruit dynamo (Spruit, 2002). In this picture, shearing grows the toroidal field, which then initiates Tayler instabilities (non-axisymmetric, pinch-like instabilities including stratification, Tayler, 1973; Spruit, 1999). This turbulently creates poloidal field components that once again shear to be toroidal. This cycle continues, creating a steady-state field.

The minimum shear needed for this process to operate can be argued simply. (See Spruit 2006 for a mathematical derivation that uses the dispersion relation from Acheson 1978.) We note that Denissenkov & Pinsonneault (2007) have recently given an alternate prescription for this same mechanism using solely heuristic arguments. Since their results have not been shown consistent with a more rigorous mathematical analysis we consider Spruit’s conclusions more reliable at this time. A vertical perturbation, $l_{z}$, is limited by buoyancy forces to be (eq. [6] from Spruit 2002)

[TABLE]

where $\omega_{\rm A}=B/[(4\pi\rho)^{1/2}R]$ is the Alfvén frequency. At small lengthscales, magnetic diffusion damps out perturbations. In the limit of $\Omega\gg\omega_{\rm A}$, which we are considering, the Tayler instability growth rate is $\sigma_{B}=\omega_{\rm A}^{2}/\Omega$, so that

[TABLE]

Combining these two relations gives the minimum $\omega_{\rm A}$ needed for the dynamo to act,

[TABLE]

During the timescale for Tayler instability, $\sigma_{B}^{-1}$, $B_{r}$ is stretched into $B_{\phi}$ by an amount

[TABLE]

The largest amplification is achieved for magnetic fields that extend the largest radial lengthscale available, so that assuming equation (36) is marginally satisfied along with the induction equation we find $B_{r}/B_{\phi}=l_{z}/R=\omega_{\rm A}/N$. Combining this with equation (39) we obtain $q=(N/\Omega)(\omega_{\rm A}/\Omega)$. Substituting this into equation (38), we find (Spruit, 2002)

[TABLE]

as the minimum shear needed for the Tayler instability to operate. Though this result is consistent with more rigorous analysis (Spruit, 2006), it should be viewed with some caution as we apply it to accreting NSs. The thin shell geometry we consider is quite different than the spherical geometries typically used when invoking the Tayler-Spruit dynamo. Simulations by Braithwaite (2006) demonstrate that the Tayler instability is strongest along the rotation axis, which is only realized at the poles in the NS case. Since we find the dynamo to be so much stronger than any hydrodynamic transport mechanisms, we consider it to be a reasonable approximation for these magnetic effects, even if its strength is decreased due to geometry.

The value of $N$ used to evaluate equation (40) depends on what is supplying the buoyancy. We follow Spruit (2002) and separately consider cases of $N_{\mu}\gg N_{T}$ and $N_{\mu}\ll N_{T}$, denoted as case 0 and 1, respectively. In case 1 non-adiabatic effects become important when $\eta/K<1$, and we must take $N_{T}^{2}\rightarrow(\eta/K)N_{T}^{2}$ (analogous to the above analysis of the secular shear instability). For this case, which dominates for most of the envelope, we find

[TABLE]

Both $q_{\rm KH}$ and $q_{\rm BC}$ are considerably above this, thus the Tayler-Spruit dynamo activates long before the onset of purely hydrodynamic instabilities.

By assuming that equation (38) is marginally satisfied Spruit (2002) derived the steady-state field strengths. We summarize the relevant prescriptions needed for our study. The steady-state azimuthal and radial field components are

[TABLE]

for case 0. For case 1,

[TABLE]

The effective viscosities, as defined by equation (31), are

[TABLE]

Although these viscosities are appropriate for angular momentum transport, mixing of material is less efficient since it requires expending work to exchange fluid elements (versus just exerting shear stresses). The mixing diffusivities are

[TABLE]

which is just equal to the effective turbulent magnetic diffusivity. In the Appendix we show that these prescriptions are consistent with energy conservation considerations.

3.2. Shearing Profile

The Tayler-Spruit dynamo transports angular momentum, causes viscous heating, and mixes material. We follow the procedure we used for the baroclinic case and assume steady-state angular momentum transport. This is a good approximation when the viscous timescale, $t_{\rm visc}=H^{2}/\nu$, is less than the timescale of accretion. For these initial estimates we focus on the $N_{T}\gg N_{\mu}$ limit (case 1) since this dominates except at compositional boundaries (which we revisit in §4). Substituting $\nu_{1}$ into equation (4), the shear rate for the Tayler-Spruit dynamo is

[TABLE]

Scaling to values appropriate for accreting NSs,

[TABLE]

Note the scalings with $\dot{m}$ and $\Omega$. At high $\dot{m}$ angular momentum is fed into the star faster, creating more shear. At low $\Omega$ the shear is greater because of the larger relative angular speed between the accreted material and the NS. These are generic feature we expect for any viscous mechanism (compare the scaling of $q_{\rm TS}$ with $q_{\rm KH}$ and $q_{\rm BC}$ from eqs. [15] and [30], respectively). The viscosity that gives the smallest shear rate will likely be the most important at a given depth. Using this criterion, we find that the Tayler-Spruit dynamo is dominant for densities $\rho\gtrsim 3\ {\rm g\ cm^{-3}}\ T_{8}^{-5/2}\dot{m}_{0.1}^{2}\Omega_{0.1}^{-3}$, or $\rho\gtrsim 5\times 10^{3}\ {\rm g\ cm^{-3}}$ for $T\approx 5\times 10^{6}\ {\rm K}$. At shallower depths, Kelvin-Helmholtz instabilities damp the shear, which is consistent with the use of hydrodynamic instabilities by Inogamov & Sunyaev (1999) for understanding the initial spreading of accreted material in the boundary layer.

The steady-state magnetic field components are

[TABLE]

Cumming & Bildsten (2000) argue that if $B_{r}\gtrsim 10^{6}\ {\rm G}$ it would become dynamically important to determining the drift of X-ray burst oscillations (Muno et al., 2002). The interaction of such fields with the shearing of an X-ray burst from an accreting millisecond pulsar has been explored by Lovelace et al. (2007). At high accretion rates, $B_{r}$ increases and comes close to this limit, suggesting that these magnetic fields may be important for understanding the dynamics of X-ray burst oscillations, even from NSs that do not show persistent pulsations (which are preferentially seen at high accretion rates, Muno et al., 2004).

3.3. Viscous Heating

Viscous shearing heats the surface layer, which can also be thought of as the rate of magnetic energy destruction as the dynamo builds and destroys magnetic field (Maeder & Meynet, 2004). The heating rate per unit mass is

[TABLE]

which for the Tayler-Spruit dynamo becomes

[TABLE]

To put this number into perspective we express it in terms of the energy released per accreted nucleon. This is found by multiplying the above result by the total mass per unit area down to depth of interest, $y\approx H\rho$ (the column depth), and dividing by $\dot{m}$. We then find

[TABLE]

In comparison, the thermal energy per nucleon at $10^{8}\ {\rm K}$ is $\approx~{}10\ {\rm keV\ nuc^{-1}}$ and burning of helium into carbon releases $\approx~{}0.6\ {\rm MeV\ nuc^{-1}}$. We conclude that viscous processes do not heat the layer sufficiently to alter X-ray bursts.

3.4. Turbulent Mixing

It is not immediately clear how much of an observational impact is provided by the magnetic fields and viscous heating estimated above. In contrast, as we shall show, shear mixing produced by the dynamo could have important ramifications for the structure and composition of the surface layers. For this reason, we devote most of the remainder of our study to investigating the consequences of mixing.

We first consider some estimates that highlight mixing’s importance. For this, we parametrize the mixing diffusivity as $\alpha_{\rm TS}D$, where $D$ is given by equations (48) or (49) and $\alpha_{\rm TS}$ is a factor of order unity that accounts for uncertainties in the Tayler-Spruit prescription. The features we find are general enough that other viscosities can be incorporated by increasing or decreasing $\alpha_{\rm TS}$.

The material becomes fully mixed over a scale height, $H$, in a time $t_{\rm mix}=H^{2}/(\alpha_{\rm TS}D)$, which gives for case 1 ($N_{T}\gg N_{\mu}$),

[TABLE]

Accretion is also advecting material downward, which happens in a timescale $t_{\rm acc}=y/\dot{m}\approx\rho H/\dot{m}$,

[TABLE]

The ratio of these timescales is

[TABLE]

The scaling with $\rho$ shows that mixing is generally more important than advection at shallow depths ($\rho\lesssim 3\times 10^{7}\ {\rm g\ cm^{-3}}$), while the scaling with $\Omega$ shows that mixing is more important for slower spinning NSs (as expected from our discussion of shear rates in §3.2).

4. Accumulating, Non-mixed Models

We now use numerical calculations to consider the angular momentum transport through the surface layers. In this section we calculate accumulating models without directly incorporating the effects from mixing. This verifies many of the analytic estimates derived above and motivates when mixing must be included (which is done for models in §5).

4.1. Shear Profile Calculations

We calculate the envelope profile for helium accumulating on an iron ocean (which represents the iron-peak ashes from previous X-ray bursts). As discussed in §2.2, we approximate the surface as having a constant gravitational acceleration, and plane-parallel geometry. In addition to using $z$ as our radial coordinate, we find it useful to use the column depth, $y$, defined as $dy=-\rho dz$, giving a pressure $P=gy$ from hydrostatic balance. We solve for $\rho$ using the analytic equation of state from Paczyński (1983). For the liquid phase, when $1\leq\Gamma\leq 173$ where $\Gamma=(4\pi n_{i}/3)^{1/3}Z^{2}e^{2}/(k_{\rm B}T)$, $Z$ is the charge per ion, and $n_{i}$ is the ion density, we include the ionic free energy of Chabrier & Potekhin (1998).

During the accumulating phase, a negligible amount of helium burning takes place, so we assume that the flux is constant and set by heating from the crust. Previous studies of accreting NSs have shown that the interior thermal balance is set by electron captures, neutron emissions, and pycnonuclear reaction in the inner crust (Miralda-Escudé et al., 1990; Zdunik et al., 1992; Bildsten & Brown, 1997; Brown & Bildsten, 1998; Brown, 2000, 2004) which release $\approx 1\textrm{ MeV}/m_{p}\approx 10^{18}\ {\rm ergs\ g}^{-1}$ (Haensel & Zdunik, 1990, 2003). Depending on the accretion rate and thermal structure of the crust, this energy is either conducted into the core or released into the ocean such that for an Eddington accretion rate up to $\approx 92\%$ of the energy is lost to the core and exits as neutrinos (Brown, 2000). We therefore set the heating to $150\ {\rm keV\ nuc^{-1}}$, giving a flux $2.2\times 10^{21}\ {\rm ergs\ cm^{-2}\ s^{-1}}\langle\dot{m}\rangle_{0.1}$, where $\langle\dot{m}\rangle_{0.1}$ is the time average accretion rate in units of $0.1\dot{m}_{\rm Edd}$ and the average is over timescales of order the thermal time of the crust (many years). For simplicity we assume $\langle\dot{m}\rangle=\dot{m}$. Recent calculations by Gupta et al. (2006) suggest that heating is stronger than previously thought, but not sufficiently high to qualitatively change our results. We ignore the additional flux from compressional heating because it only contributes $\sim c_{p}T\sim 10\ {\rm keV\ nuc^{-1}}$ (Bildsten, 1998a).

In a one-zone estimate ignition occurs at the base of the helium layer when

[TABLE]

(Fujimoto et al., 1981) where $\epsilon_{3\alpha}$ is the heating rate from triple-$\alpha$ reactions (for which we use the rate from Fushiki & Lamb 1987 for numerical calculations),

[TABLE]

and the derivatives are both taken at constant pressure. We consider models where the column of accreted helium $y_{\rm acc}$ that is just at this ignition point. This allows the maximum amount of time for angular momentum redistribution during the accumulating phase. The He/Fe boundary is assumed sharp, since the timescale for diffusion between the layers is much longer than the timescale of accumulation (Brown et al., 2002).

We solve for the temperature profile by integrating the radiative diffusion equation,

[TABLE]

The opacity is set using electron-scattering (Paczyński, 1983), free-free, and conductive opacities (Schatz et al., 1999). There is a sharp change in opacity at the He/Fe boundary (because the high $Z$ of the iron makes it more opaque), which means that the top of the iron layer is convectively unstable. The characteristic convective velocity estimated from $F\approx\rho V_{\rm conv}^{3}$ is $V_{\rm conv}\sim 10^{5}\ {\rm cm\ s^{-1}}$, which is much less than the local sound speed of $\sim 10^{8}\ {\rm cm\ s^{-1}}$. This means that the convection is very efficient, and we therefore simply set $(d\ln T/d\ln P)_{*}=\nabla_{\rm ad}$ in the convective region. This convective region has not been noted in previous studies, though a super-adiabatic temperature gradient is apparent in the figures shown in both Brown et al. (2002) and Cumming (2003). It is not yet clear how this impacts the bursting properties of NSs, and we delay exploring this in detail for a future study. Nevertheless, we must include convection so as to have accurate differentially rotating profiles.

In Figure 2 we plot the temperature profile, $K$ and $\eta$, and $N$ for three accumulating models. The magnetic diffusivity is set using the conductivity from Schatz et al. (1999). The base of the accumulating helium is set where unstable helium ignition occurs, as designated by the thick long-dashed line in the top panel (eq. [60]). The convective zone begins just below this and is bracketed by solid circles at its top and bottom. The convective zone is also seen in the plot of $N$ since $N_{T}$ is effectively zero in this region. The change in composition at the He/Fe boundary gives a large buoyancy contribution, which we estimate as

[TABLE]

where $\Delta\ln\mu\approx 0.44$ is the logarithmic change in the mean molecular weight at the boundary. We denote this by a solid square for each model in the bottom panel. For the majority of the profile the Tayler-Spruit dynamo is given by case 1 ($N_{T}\gg N_{\mu}$). Since $K\gg\eta$, the perturbations creating the dynamo are non-adiabatic for case 1, as we described in §3.1. The only place case 0 ($N_{\mu}\gg N_{T}$) is important is at the He/Fe boundary.

We next solve for the shear rates by solving equation (4) with either $\nu_{0}$ or $\nu_{1}$, depending on which case is appropriate. We assume angular momentum transport does not affect the thermal or compositional structure, in other words, that the transport is just happening “in the background.” This allows us to assess which effects are crucial for subsequent iterations that include angular momentum transport in the actual structure calculation.

In Figure 3 lines denoted the profiles calculated assuming case 1 of the Tayler-Spruit dynamo. We assume that the convection is effectively instantaneous in transporting material and angular momentum, thus all the quantities become very small in the convective region. This is appropriate since the convective overturn timescale, $H/V_{\rm conv}\sim 10^{-4}\ {\rm s}$, is much less than the timescales for mixing or accretion. The squares indicate the corresponding values due to the compositional discontinuity found by using case 0. Since the viscous timescale, $t_{\rm visc}=H^{2}/\nu$, is much less than the timescale it took to accrete to this base column $t_{\rm acc}=y_{\rm acc}/\dot{m}\sim 10^{3}-10^{5}\ {\rm s}$ (for the range of $\dot{m}$ considered) our steady-state assumption is valid.

These calculations highlight the importance of the compositional jump since $t_{\rm visc}$, $q$, and $dE/d\ln y$ are all amplified here. This is because the large buoyancy reduces angular momentum transport across this boundary. The viscous time is nearly independent of accretion rate in regions where $N_{T}\gg N_{\mu}$ because the viscosity in case 1 is independent of $q$. The energy deposition is always much smaller than the heat coming the crust, and it falls off somewhat faster with depth than expected from the estimate presented in §3.3 due to electron degeneracy effects decreasing $N_{T}$. New heat sources at a depth of $\approx 10^{12}\ {\rm g\ cm^{-2}}$ could ease the difficulty calculations have in recreating the low ignition columns needed to explain superburst recurrence times (Cumming et al., 2006), but this viscous heating is not nearly enough to correct this problem.

In Figure 4 we plot the magnetic fields found within the radiative zones for the models from Figure 3. We assume that within the convective zone the dynamo is not able to operate, and do not calculate a magnetic field here. It is also interesting to compare these fields to those derived by Cumming et al. (2001), who calculated the steady-state fields expected when Ohmic diffusion balances advection through accretion. Their principal result was that the steady-state horizontal field drops by $\approx\dot{m}/0.02\dot{m}_{\rm Edd}$ orders of magnitude from the crust up through the ocean. In contrast, the magnetic fields we find are nearly constant with depth. We therefore do not expect these fields to persist if the accretion ceases for a time and instead to be expelled on an Ohmic diffusion time ($\sim\ {\rm days}$ near the top of the ocean). When the dynamo is active, the steady-state is reached quickly enough that Ohmic diffusion can be ignored.

In Figure 5 we compare what happens as the spin is changed by plotting $\Omega=0.03$, 0.1, and $0.3\Omega_{\rm K}$ ($67$, $220$, and $670\ {\rm Hz}$), all for $\dot{m}=0.1\dot{m}_{\rm Edd}$. Note that the shearing is most dramatic at smaller $\Omega$. In fact, the shearing and heating profiles are very sensitive to the value of $\Omega$, as was demonstrated by the analytic estimates. We do not plot the associated magnetic field for these models since we have already plotted some examples in Figure 4 and the analysis of §3.2 provides adequate estimates.

These plots of $q$ show that very little shearing is present. To emphasize this fact, in Figure 6 we plot the actual spin frequency as a function of depth found by integrating $q$ for a range of accretion rates. The discontinuity at the He/Fe boundary complicates this estimate. Noting that

[TABLE]

we approximate the spin change at this boundary as

[TABLE]

The change of spin is generally $\lesssim 0.5\ {\rm Hz}$ across the accumulating layer, with the majority of the spin change occurring at the compositional boundary. The layer is very nearly in uniform rotation.

4.2. Mixing and the Compositional Barrier

The above calculations confirm that the viscosity is too large for either appreciable shearing or viscous heating. As demonstrated in §3.4, mixing should be important, but Figures 3 and 5 argue that we must take into account the large $N_{\mu}$ at the He/Fe boundary. Assuming that $N_{\mu}$ scales like equation (63), we estimate

[TABLE]

We substitute this into $\nu_{0}$ to solve for $q$ using equation (4), which is used to estimate a mixing timescale at the boundary,

[TABLE]

Using equation (58) we find a ratio of

[TABLE]

Unlike in equation (59), this new ratio depends on $\dot{m}$. This is because the viscosity in case 0 has a different dependence on $q$. Now as $\dot{m}$ is increased $t_{\rm mix}$ decreases faster than $t_{\rm acc}$. Above a critical accretion rate of

[TABLE]

mixing can no longer be ignored. For densities and temperatures expected at the base of the accumulating layer, this critical accretion rate lies in the range of $0.1-1.0\dot{m}_{\rm Edd}$. Note that this depends very strongly on the spin rate, $\dot{m}_{\rm crit,1}\propto\Omega^{3}$, thus we expect the slower spinning NSs to be considerably more affected by mixing.

To test these analytic estimates we compare the mixing and accretion timescales in Figure 7. From the top panel to the bottom panel we increase $\dot{m}$ (fixing $\Omega=0.1\Omega_{\rm K}$). At low $\dot{m}$, the He/Fe boundary at $y_{\rm acc}$ (shown by the filled squares) acts as a barrier to mixing since $t_{\rm mix}>t_{\rm acc}$ at this depth. When this happens it is a good approximation to ignore mixing and assume two separate layers during the accumulation phase. As $\dot{m}$ increases, $t_{\rm mix}$ at $y_{\rm acc}$ becomes less and less until finally $t_{\rm mix}<t_{\rm acc}$, so that material should be mixed past $y_{\rm acc}$. When this occurs, our accumulating model can no longer ignore mixing. The mixing between helium and iron occurs down to a depth where $t_{\rm mix}$ is equal to the length of time accretion has been taking place, $t_{\rm acc}=y_{\rm acc}/\dot{m}$. The key point we want to emphasize is that because of the buoyancy barrier, the effect of mixing turns on abruptly, and when it does, mixing will occur well past $y_{\rm acc}$.

5. The Effects of Turbulent Mixing

Once $t_{\rm mix}<t_{\rm acc}$ at the buoyancy barrier, the compositional profile of the NS is very different, which we now explore. We treat the mixing as complete, which we diagram in Figure 8 and summarize here. Material accretes at a rate $\dot{m}$ for a time $t_{\rm acc}$ with a helium mass fraction $Y_{0}$, supplying a column of material $y_{\rm acc}=\dot{m}t_{\rm acc}$. The total column of helium that accretes during this time is therefore $Y_{0}y_{\rm acc}$. Mixing causes this newly accreted material to mix past $y_{\rm acc}$ to the mixing depth $y_{\rm mix}$ defined as where $t_{\rm mix}=t_{\rm acc}$. The helium is fully mixed down to this depth, resulting in a diluted mass fraction $Y_{\rm mix}=Y_{0}y_{\rm acc}/y_{\rm mix}$ within the mixed layer. In the following sections we consider two scenarios that can result from the mixing: (1) unstable ignition of the mixed fuel, and (2) stable burning when the material mixes to sufficient depths.

5.1. Numerical Calculations of Mixed Ignition

As the column of accreted material grows, it can still reach the correct conditions for unstable ignition, but mixing causes two changes: (1) the ignited layer has a diluted helium fraction, $Y_{\rm mix}<Y_{0}$, and (2) the ignition takes place at the base of the mixed layer, at a depth $y_{\rm mix}>y_{\rm acc}$, resulting in a recurrence time for ignition much less than when mixing is not included.

Both of these effects are easiest to explore using a semi-analytic model. In this section all calculations use $Y_{0}=1$. We consider other values of $Y_{0}$ for our analytic estimates in the next section. We solve for the mixed accumulating structure by first assuming an amount of accretion $y_{\rm acc}$, which for a given $\dot{m}$ implies an accretion time $y_{\rm acc}/\dot{m}$. We integrate the radiative diffusion equation (eq. [62]) down to a depth where $t_{\rm mix}=t_{\rm acc}$ giving $y_{\rm mix}$, where we assume a constant flux profile as was discussed in §4 since little helium burning takes place during accumulation. We then estimate the mixed helium fraction as $Y_{\rm mix}=y_{\rm acc}/y_{\rm mix}$. This estimate is improved by substituting $Y_{\rm mix}$ back into our envelope integration, and iterating until $Y_{\rm mix}$ converges. The non-helium component in the layer is taken to be iron.

In Figure 9 we plot the resulting profiles for a NS accreting at $\dot{m}=0.1\dot{m}_{\rm Edd}$ and $\Omega=0.1\Omega_{\rm K}$. In the four panels we consider values of $y_{\rm acc}$ of $10^{6}$, $3\times 10^{6}$, $10^{7}$, and $4\times 10^{7}\ {\rm g\ cm^{-2}}$ (from left to right and then up to down, denoted by the filled circles), which is meant to mimic the accumulation of fuel on the NS surface. For each integration the envelope profile continues down to a depth $y_{\rm mix}$, which also gives $Y_{\rm mix}$ as displayed in each panel. The last model reaches the conditions necessary for ignition at the base of the mixed layer. The recurrence time for these mixed-ignition models is much less than for those without mixing. For the plotted model, the recurrence time is $t_{\rm rec}=4\times 10^{7}\ {\rm g\ cm^{-2}}/1.5\times 10^{3}\ {\rm g\ cm^{-2}\ s^{-1}}\approx 44\ {\rm minutes}$. In contrast the model shown in Figure 2 with the same accretion rate of $\dot{m}=0.1\dot{m}_{\rm Edd}$ has $t_{\rm rec}\approx 1.5\ {\rm days}$! The shorter recurrence time is not only due to mixing carrying helium down to deeper depth, but also a change in the thermal profile. Since significant iron is mixed up into the accumulating material, the free-free opacity, which scales $\kappa_{\rm ff}\propto Z^{2}/A$, where $Z$ and $A$ are the charge and mass per nucleon, respectively, is now the dominate opacity mechanism. The accumulating layer is more opaque and therefore hotter for a given flux in comparison to the pure-helium models considered before, which contributes to the shallow ignition depths.

We calculate $t_{\rm rec}$ for a grid of models with various $\dot{m}$ and $\Omega$ in Figure 10. The recurrence time is shorter for stronger mixing, which occurs at high $\dot{m}$ or low $\Omega$, and can be as short as $\approx 5-30\ {\rm minutes}$. We warn though that all of these calculations assume that complete mixing can occur, and as we already showed in §4, buoyancy may prevent this (eq. [69]). Nevertheless, it is interesting to calculate the mixed-ignition conditions for a wide range of parameter space because the conditions left from previous bursts may vary, and quantities such as $\Delta\ln\mu$ may be smaller at times if, for example, there is incomplete burning in a previous burst. At sufficiently high $\dot{m}$ or low $\Omega$ the ignition takes place at high enough temperatures that the envelope does not ignite unstably. We consider this case in more depth in the following section.

The short recurrence times that we find are similar to some seen for multiple bursts (Galloway et al., 2006, and references therein). However, our model cannot explain the energetics of these bursts. The energetics are typically quantified in a distance independent measure, the so-called $\alpha$-value, which is the ratio of energy released in the persistent emission between bursts to the energy of the burst itself. For pure helium ignition $\alpha\sim 100$. In contrast, the $\alpha$-values for the short recurrence time bursts are typically $\sim 10$ (Boirin et al., 2007). Since the only helium that will burn in our mixed ignition models is that accreted since the last outburst, we always find $\alpha\sim 100$. Therefore we still require an additional nuclear energy source, such as incomplete burning from the previous X-ray burst, to explain such a low $\alpha$-value. This in fact may not be a problem because incomplete burning may naturally explain the short recurrence times, since this would lead to smaller compositional gradients and therefore stronger mixing.

We plot the $Y_{\rm mix}$ as a function of $\dot{m}$ and $\Omega$ in Figure 11. The crosses correspond to models that ignite unstably. The filled circles are stable accreting models that are discussed on §5.3 and §5.4. The implications of this mixed ignition for the unstable burning during the X-ray burst can be easily tested by more sophisticated numerical simulations by just considering a mixed accumulating column, for example, by artificially setting $Y\approx 0.1-0.6$ as shown in this Figure 11. Since viscous heating is negligible, these initial tests do not need to resolve the shearing profiles.

5.2. Analytic Estimates of Mixed Ignition

We now estimate the properties of the mixed ignition models. These solutions directly show how mixing depends on the properties of the accumulating layer, in particular the prefactor $\alpha_{\rm TS}$ and $Y_{0}$, without having to consider a multitude of models. The effects of a free-free opacity are important in deriving the correct atmospheric conditions. To include this analytically, we use the free-free opacity from Schatz et al. (1999), simplified to a one-component plasma and with the dimensionless Gaunt factor set to unity,

[TABLE]

where we have estimated $\mu_{e}\approx 2$ as is correct within 10% for the any of the elements of interest. Integrating the radiative diffusion equation (eq. [62]), assuming a constant flux of $F=10^{21}\ {\rm ergs\ cm^{-2}\ s^{-1}}F_{21}$ and an ideal gas equation of state, the temperature as a function of column $y$ is

[TABLE]

where $y_{8}\equiv y/10^{8}\ {\rm g\ cm^{-2}}$. Using equation (57), we find the mixing timescale as a function of $y$,

[TABLE]

Note that $t_{\rm mix}\propto y^{1.37}$ is a higher power than $t_{\rm acc}=y/\dot{m}\propto y$. This explicitly shows that mixing dominates at lower $y$, but accretion always wins at some depth. Setting $t_{\rm mix}$ equal to $t_{\rm acc}=y_{\rm acc}/\dot{m}$ gives the depth where mixing can extend to for a given column of accreted material $y_{\rm acc}$,

[TABLE]

where $y_{\rm acc,8}=y_{\rm acc}/10^{8}\ {\rm g\ cm^{-2}}$. The mixed helium fraction down to this depth is

[TABLE]

Since $Y_{\rm mix}\propto y_{\rm acc}^{0.27}$ the strength of mixing decreases ($Y_{\rm mix}$ gets larger) as more material accreted, which was demonstrated by the four panels in Figure 9.

We next estimate what ignition depth is expected for this fully mixed accumulating layer. The energy generation rate for triple-$\alpha$ burning is approximated as

[TABLE]

where $f$ is factor that accounts for screening effects. To make progress analytically we expand the exponential as $\exp\left(-44/T_{8}\right)\approx 7.95\times 10^{-10}(T_{8}/2.1)^{21}$. Using our temperature profile (eq. [71]), the condition that $d\epsilon_{3\alpha}/dT=d\epsilon_{\rm cool}/dT$ (eq. [60]) implies an ignition depth of

[TABLE]

Setting $y_{\rm ign}=y_{\rm mix}$ from equation (73), we solve for $y_{\rm acc}$, the critical column of material that must be accreted to cause ignition. This is then substituted back into equations (73) and (74) to find that ignition occurs at a depth

[TABLE]

with a composition of

[TABLE]

Comparing this with the numerical calculation is easiest if we assume $F$ is set by $\dot{m}$ ($F_{21}=2.2\dot{m}_{0.1}$) as well as scaling $Z^{2}/A\approx 12$ and $\mu\approx 2.1$ as appropriate for the iron-rich composition. This gives

[TABLE]

The recurrence time is $t_{\rm rec}=Y_{\rm mix,ign}y_{\rm mix,ign}/\dot{m}$, resulting in

[TABLE]

Equations (79) and (80) confirm the scalings found for the numerical calculations in Figures 10 and 11. We have plotted the $Y_{\rm mix}\propto\dot{m}^{-0.21}$ scaling in Figure 11 to emphasize this. These analytic results also show how strongly these results depend on the parameter $\alpha_{\rm TS}$, for which $Y_{\rm mix,ign}$ is especially sensitive.

5.3. Steady-State Mixing and Burning

Another possibility is that the helium is mixed and burned by triple-$\alpha$ reactions in steady-state, leading to stable burning. The basic idea is similar to that described above for mixed ignition, except now the depth of the accumulating layer is set by the helium burning timescale, $t_{3\alpha}=E_{3\alpha}/(Y\epsilon_{3\alpha})$, where $E_{3\alpha}=5.84\times 10^{17}\ {\rm ergs\ g^{-1}}$ is the energy per mass released from this burning. As shown in Figure 12, material is mixed to sufficient depths where $t_{3\alpha}$ (which decreases with depth) is equal to $t_{\rm mix}$ (which increases with depth), which defines a mixing (or burning) depth $y_{\rm mix}=y(t_{\rm mix}=t_{3\alpha})$. During a mixing timescale, the amount of material that is accreted is $y_{\rm acc}=\dot{m}t_{\rm mix}(y_{\rm mix})$, so that the total column of helium that has been accreted is $Y_{0}\dot{m}t_{\rm mix}(y_{\rm mix})$. This is diluted over a depth $y_{\rm mix}$, so that the mixed helium fraction is $Y_{\rm mix}=Y_{0}y_{\rm acc}/y_{\rm mix}$. This is all occurring in steady-state, material moves through the mixed layer at a rate $\dot{m}$, but this layer does not move up or down in pressure (column) coordinates, as the burning is stable.

The basic results of steady-state mixing and burning are best shown using a simple numerical model. The equation that describes helium continuity, including depletion by triple-$\alpha$ burning, becomes, in the plane-parallel limit (Fujimoto, 1993),

[TABLE]

We assume changes in $\rho$ and $D$ with depth are small in comparison to changes in $Y$, and, following our derivation of the angular momentum equation in §2.2, we take the steady state limit to derive

[TABLE]

In the limit where $t_{\rm acc}\gg t_{\rm mix}$ that we are interested in, the term on the left-hand side can be dropped. Finally, making the approximation that $\rho Dd/dz~{}\approx~{}y/t_{\rm mix}$, we find

[TABLE]

This equation mimics the properties we expect from mixing. When mixing is strong $t_{\rm mix}/t_{3\alpha}\ll 1$, and $dY/dy\approx 0$, the composition does not change with depth. At the depth where $t_{\rm mix}\approx t_{3\alpha}$, $dY/dy\approx-Y/y$ and the helium is depleted exponentially. All the helium burns into carbon, so that carbon has a mass fraction $X_{12}=1-Y$.

The envelope profiles are found from simultaneously integrating three differential equations: (1) radiative transfer, equation (62), (2) the entropy equation, $dF/dy=-\epsilon_{3\alpha}$, and (3) continuity of helium, equation (83). We integrate using a shooting method, but first we must set three boundary conditions for the flux, temperature, and helium mass fraction. Since all of the accreted helium must burn if the envelope is in steady-state, we set the surface flux to $F=Y_{0}E_{3\alpha}\dot{m}+F_{c}$, where $F_{c}=150\ {\rm keV\ nuc^{-1}}\langle\dot{m}\rangle$ (as discussed in §4.1). The surface temperature is set from the radiative zero solution (Schwarzschild, 1958). The helium abundance in the mixed region, $Y_{\rm mix}$, is an eigenvalue. It is varied until shooting gives the correct base flux of $F_{c}$. This is easily found through iteration since when $Y_{\rm mix}$ is set too large, too much burning occurs, and the base flux is too small (and vice versa for small $Y_{\rm mix}$).

In Figure 13 we plot a steady-state envelope using $Y_{0}=1$, $\dot{m}=0.3\dot{m}_{\rm Edd}$, $\Omega=0.1\Omega_{\rm K}$, and $\alpha_{\rm TS}=1$. The shooting method demonstrates that the initial helium abundance within the mixed layer is $Y_{\rm mix}=0.097$. The top panel shows the temperature profile (thin solid line), and the critical curve for stability where $d\epsilon_{3\alpha}/dT=d\epsilon_{\rm cool}/dT$ (thick dashed line). The bottom panel shows that the majority of triple-$\alpha$ burning occurs at $\approx 3\times 10^{8}\ {\rm g\ cm^{-2}}$ (thick solid line), which is where the helium is depleted. Comparing the middle and bottom panel shows that the majority of the burning takes place near where $t_{\rm mix}=t_{3\alpha}$ (as required by construction), which is deeper than where $t_{\rm acc}=t_{3\alpha}$ (the normal condition for steady burning).

To make sure the steady-state models we find are physically realizable, we must check the thermal stability of the helium burning. In Figure 14 we compare the quantities $d\epsilon_{3\alpha}/d\ln T$ and $d\epsilon_{\rm cool}/d\ln T$ for three different accretion rates. If the cooling derivative is always larger, then the model is thermally stable. We find stable accretion at $\dot{m}$’s considerably less than the stable accretion rate of $\dot{m}\approx 10\dot{m}_{\rm Edd}\approx 2\times 10^{6}\ {\rm g\ cm^{-2}\ s^{-1}}$ expected for pure helium accretion estimated without mixing (Bildsten, 1995, 1998a). Models that are found to be stable in this way are plotted as filled circles in Figure 11 (from §5.1). Only at large $\dot{m}$ and small $\Omega$ are the models found to be stable. If we where to increase $\alpha_{\rm TS}$, a wider range of the models in Figure 11 become stable.

5.4. The Critical $\dot{m}$ for Stability

Since we have found a set of models that can stably accrete, mix, and then burn helium, it is interesting to ask what accretion rates and spins are required for this to occur, and how does it depend on parameters such as $\alpha_{\rm TS}$.

First we must derive the correct condition for stability including the fact that free-free opacity is important in setting the radiative profile. The one-zone condition for stability at the base of the mixed layer is

[TABLE]

where the derivatives are taken at constant pressure. For $\epsilon_{3\alpha}\propto T^{\zeta}\rho^{\chi}$, where $\zeta=44/T_{8}-3$, stability requires (Bildsten, 1998a)

[TABLE]

Substituting the scalings for a free-free opacity and ideal gas equation of state we find $T_{8}>3.26$ is required for stability.

By substituting the analytic form we found for the mixed ignition depth (eq. [77]) into the temperature profile (eq. [71]) we can find the temperature at the base of the mixed layer,

[TABLE]

By simply asking when $T_{\rm ign}>3.26\times 10^{8}\ {\rm K}$, we derive a stabilizing $\dot{m}$ of

[TABLE]

This is in reasonable agreement with Figure 11, which shows stability occurs for $0.9\dot{m}_{\rm Edd}\lesssim\dot{m}\lesssim 1.0\dot{m}_{\rm Edd}$ for $\Omega=0.1\Omega_{\rm K}$. It is interesting that $\dot{m}_{\rm crit,2}$ is near (within an order of magnitude) a value where bursts are observed to change. It is conceivable that this mechanism may act to stabilize X-ray bursts for the hydrogen-rich accreting systems, and given the strong scaling $\dot{m}_{\rm crit,2}$ has with $\alpha_{\rm TS}$ it is possible that more detailed calculations could give results that agree even better with the critical accretion rates that are observed. We discuss this idea in more detail in the following section.

6. Discussion and Conclusion

We have revisited the problem of angular momentum transport in the surface layers of accreting NSs. We found that the hydrodynamic instabilities used by Fujimoto (1993) in a previous study are dwarfed by the magnetic effects of the Tayler-Spruit dynamo. The large viscosity provided by this process results in a very small shear rate and negligible viscous heating. The turbulent mixing is sufficiently large to have important consequences for X-ray bursts. We constructed simple models, both analytic and numerical, to explore mixing for pure helium accretion. From these models we can make a few conclusions that are likely general enough to apply to most viscous mechanisms. As a guide, we show the different burning regimes we find in Figure 15. These can be summarized as follows:

• •

Mixing is strongest at large $\dot{m}$ (when angular momentum is being added at greater rates) and small $\Omega$ (which gives a larger relative angular momentum between the NS and accreted material).

• •

Mixing has trouble overcoming the buoyancy barrier at chemical discontinuities. But once mixing breaks through this, it extends down to a depth where $t_{\rm mix}=t_{\rm acc}$, which is generally much deeper than the accreted column. This means that mixing should turn on abruptly (as a function of either $\dot{m}$ or $\Omega$). It also means that the importance of mixing depends on the particular ashes left over from previous bursts. If, for example, incomplete burning results in small compositional gradients, mixing would be important in subsequent bursts.

• •

Mixing of freshly accreted material with the ashes from previous X-ray bursts can lead to two new effects. First, the layer may ignite, but now in a mixed environment with a short recurrence time of $\sim 5-30\ {\rm minutes}$. Second, if the mixing is strong enough, accreted helium can mix and burn in steady-state, quenching X-ray bursts. Both of these regimes have observed analogs, namely the short recurrence time bursts (for example in, Boirin et al., 2007) and the stabilization of bursting seen at $\approx 0.1\dot{m}_{\rm Edd}$ (Cornelisse et al., 2003).

The mixed ignition case can be studied easily using the current numerical experiments (e.g., Woosley et al., 2004) by just artificially accreting fuel with a mixed composition of $Y_{\rm mix}\approx 0.1-0.6$. These calculations are simplified by our conclusion that shearing and heating can be ignored, at least for initial studies. We next conclude by speculating about some of the other ramifications of turbulent mixing.

6.1. Superbursts

An ongoing mystery in the study of bursting NSs is the recurrence times for superbursts, thermonuclear ignition of carbon in the X-ray burst ashes at columns of $\approx 10^{11}-10^{12}\ {\rm g\ cm^{-2}}$ (Cumming & Bildsten, 2001; Strohmayer & Brown, 2002). This problem could be alleviated by enhanced heating from the core on the order of $1\ {\rm MeV\ nuc^{-1}}$ (Cumming et al., 2006), but this is more than is expected theoretically, even in the newest calculations (Gupta et al., 2006). Shear heating is not large enough to solve this problem, as demonstrated in Figures 3 and 5.

The regime of stable helium burning we have found may, however, create a carbon rich ocean that would assist in the ignition of superbursts. Carbon fractions of greater than 10% are needed to reproduce the lightcurves and recurrence times of superbursts (Keek et al., 2006). Calculations of rp-process burning show that unstable burning cannot give carbon fractions this high (Schatz et al., 2003). Observationally, in’t Zand et al. (2003) showed that the $\alpha$-value (see §5.1) is preferentially large for superbursting systems, indicating that some stable burning is occurring, perhaps due to the turbulent mixing we have studied.

6.2. Hydrogen-rich Accretion

One of the main deficiencies of our calculations are the simplified compositions, since most accreting NSs are expected to be accreting a fuel abundant in hydrogen. If the only crucial burning is triple-$\alpha$, such envelopes can be considered within the framework of our models by using a solar value of $Y_{0}=0.3$. Our models fail, though, when sufficient carbon is produced by triple-$\alpha$ to feedback into hydrogen burning (which burns via the hot-CNO cycle, Hoyle & Fowler, 1965). Such a scenario is interesting thought because it could potentially produce very hydrogen poor bursts.

In the standard theoretical framework for X-ray bursts, flashes should be hydrogen poor at low $\dot{m}$ (when there is sufficient time for the hot-CNO cycle to act), and then become mixed hydrogen-helium fuel at higher $\dot{m}$. Paradoxically, observations show a transition at a seemingly universal luminosity of $\approx 2\times 10^{37}\ {\rm ergs\ s^{-1}}$ (approximately $0.1\dot{m}_{\rm Edd}$), but in the opposite sense (Cornelisse et al., 2003). Cooper & Narayan (2006b) argue that their models in fact show this transition (also see Narayan & Heyl, 2003) because carbon created during helium simmering increases the the hot-CNO rate and the temperature, which decreases the temperature sensitivity of triple-$\alpha$ reactions. To explain the discrepancies between their models and more detailed numerical simulations (Woosley et al., 2004; Heger et al., 2005) requires a decrease in the breakout reactions rate of 15O$(\alpha,\gamma)^{19}$Ne (Cooper & Narayan, 2006a; Fisker et al., 2006). Unfortunately, the most recent experimental results do not support a decreased rate (Tan et al., 2007; Fisker et al., 2007). Furthermore, helium accreting systems show this same transition in bursting properties, which is not explained within their framework.

Another idea that may recreate this trend is that the fraction of the star covered by the fuel increases with the global accretion fast enough that the local accretion rate actually decreases (Bildsten, 2000). This interpretation is supported by the recent work of Heger et al. (2005), who claim that the mHz oscillations observed at around this same universal luminosity of $2\times 10^{37}\ {\rm ergs\ s^{-1}}$ (Revnivtsev et al., 2001) are a probe of the local accretion rate where bursting is occurring. The main problem with this suggestion is that it is difficult to understand how the covering of the surface can remain so anisotropic all the way down to the depths of where ignition occurs. Inogamov & Sunyaev (1999) calculated the spreading of accreted material from the equator, where the viscosity is due to a turbulent boundary layer, and find that spreading occurs orders of magnitude shallower in depth than where ignition takes place.

Since mixing gets stronger with $\dot{m}$, it may be that the observed transition is in fact hydrogen being turbulently mixed and burned analogous to what we have found for pure-helium accretion. A possible implication of such an explanation is that slowly spinning NSs are more likely to have their hydrogen depleted resulting in helium-rich bursts. Assuming that the burst oscillation frequencies are indicative of the NS spin frequency (which is close to true or exactly true for all current explanations Heyl, 2004; Lee, 2004; Piro & Bildsten, 2005; Payne & Melatos, 2006), the systems can be broken into slow spinning ($\sim 300\ {\rm Hz}$) and fast spinning ($\sim 600\ {\rm Hz}$) classes. This slowly spinning class includes 4U $1916-053$ ($270\ {\rm Hz}$, Galloway et al., 2001), 4U $1702-429$ ($330\ {\rm Hz}$, Markwardt et al., 1999) and 4U $1728-34$ ($363\ {\rm Hz}$, Strohmayer & Markwardt, 1999). Do these systems show helium-rich looking bursts, and if so, is this due to turbulent mixing destroying the hydrogen they are accreting? This can be answered by looking at the recent summary of RXTE burst observations by Galloway et al. (2006). 4U $1916-053$ has bursts consistent with helium-rich fuel, but this system also has an orbital period of $\approx 50\ {\rm min}$ (Grindlay et al., 1988). This is an “ultracompact” system in which the donor is too small to be a H-rich star, so the accretion is probably helium-rich to begin with. The other two NSs both have bursts with decay times and $\alpha$-values that suggest helium-rich fuel, which is also supported by model fits to radius expansion bursts from 4U $1728-34$ (Galloway et al., 2006). Furthermore, the bursts of 4U $1728-34$ have look very similar to those of 4U $1820-30$ (Cumming, 2003), a known ultracompact (see discussion in Podsiadlowski et al., 2002, and references therein). While it is possible that 4U $1702-429$ and 4U $1728-34$ have helium-rich donors (and are thus ultracompacts), it may also be that they are accreting hydrogen-rich fuel and show helium-rich bursts because their low spins lead to mixing. As binary parameters of these systems become better known, it will be more clear whether turbulent mixing is indeed needed to explain their burst properties.

We thank Henk Spruit for helpful discussions. This work was supported by the National Science Foundation under grants PHY 99-07949 and AST 02-05956.

Appendix A Mixing and Energy Conservation for the Tayler-Spruit Dynamo

In §3.1 we summarized the prescriptions that Spruit (2002) provides for the Tayler-Spruit dynamo. The mixing diffusivity is assumed to be equal to the turbulent magnetic diffusivity. Although this is plausible, it is not shown rigorously. Below we argue that such a scaling is consistent with energy conservation.

Consider a layer differentially rotating with a speed $\Delta V$. The equation for the energy per unit mass expresses the fact that energy from differential rotation can go into either viscous shearing or vertical mixing

[TABLE]

where $D=H^{2}/t_{\rm mix}$ is the mixing diffusivity. The ratio of the two right hand terms is called the flux Richardson number (Fujimoto 1988)

[TABLE]

where $Ri\equiv N^{2}/(q^{2}\Omega^{2})$ is the Richardson number. The flux Richardson number is interpreted as the ratio of energy that goes into mixing versus heating (and usually taken to be $Rf\sim 0.1-1$). From this result we can solve for the mixing diffusivity

[TABLE]

This matches (up to a factor of order unity, $Rf/4$) what Spruit (2002) gives for the mixing diffusivity. The factor of $Ri$ in the denominator takes into account the difficulty in overcoming buoyancy to do mixing. The larger the buoyancy, the larger $Ri$ is and the smaller the mixing diffusivity.