ACCEPTED FOR PUBLICATION IN THE ASTROPHYSICAL JOURNAL
Preprint typeset using LATEX style emulateapj v. 08/22/09
TURBULENT MIXING IN THE SURFACE LAYERS OF ACCRETING NEUTRON STARS
ANTHONY L. PIRO1 AND LARS BILDSTEN
Kavli Institute for Theoretical Physics, Kohn Hall, University of California, Santa Barbara, CA 93106;
piro@kitp.ucsb.edu, bildsten@kitp.ucsb.edu
Accepted for publication in The Astrophysical Journal
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, where toroidal magnetic
field growth and Tayler instabilities balance to support a steady-state magnetic field. The radial and azimuthal
components have strengths of ∼ 105 G and ∼ 1010 G, respectively. This field provides a Maxwell stress on the
shearing surface layers, which leads to nearly uniform rotation at the depths of interest for X-ray bursts (near
densities of ≈ 106 g cm−3). 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 (when the spin is small, the
relative speed of the accreted material is larger), 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.
Subject headings: accretion, accretion disks — stars: magnetic fields — stars: neutron — X-rays: bursts —
X-rays: stars
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 ac-
cretion 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 ma-
jority of the rotational kinetic energy of the accreted material
is dissipated at low densities (∼ 10−4 − 10−1 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
1 Current address: Astronomy Department and Theoretical Astrophysics
Center, 601 Campbell Hall, University of California, Berkeley, CA 94720;
tpiro@astro.berkeley.edu
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 observa-
tions of these bursts. Though a simple limit cycle picture has
been successful in qualitatively explaining the bursts’ primary
characteristics, including their energies (1039 − 1040 ergs), re-
currence times (hours to days), and durations (∼ 10 − 100 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 suc-
cessive X-ray bursts with recurrence times as short as ∼ 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 in-
stabilities to create a steady-state magnetic field, which pro-
vides a torque on the shearing layers that is larger than any
purely hydrodynamic mechanism. The result is nearly uni-
form 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
http://arxiv.org/abs/0704.1278v1
2 PIRO & BILDSTEN
be non-negligible and in some cases may mix fresh mate-
rial with the ashes of previous bursts. The key is whether
the strong buoyancy due to the larger density of the ashes be-
low 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 unsta-
ble ignition. The timescale for ignition of such bursts is short
enough (∼ 5−30 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 tempera-
tures) 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 for-
mula that estimates what spins and accretion rates necessary
to transfer between these two regimes (eq. [77]).
1.1. The Basics of Angular Momentum Transport
Before beginning more detailed analysis, it is useful to
present some general equations for angular momentum trans-
port. 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 ex-
pected rotation profile from the accretion disk, through the
boundary layer, and into the NS interior. Material accreted
at a rate Ṁ reaches the NS surface with a nearly Keplerian
spin frequency ofΩK = (GM/R
3)1/2 = 1.4×104 s−1 M
1.4 R
where M1.4 ≡ M/1.4M⊙ and R6 ≡ R/10
6 cm, which has a ki-
netic energy per nucleon of ≈ 200 MeV nuc−1. The major-
ity of this energy is dissipated in a boundary layer of thick-
ness HBL ≪ 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 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 liq-
uid 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 ρ≈ 106 g cm−3, which is well below the boundary layer.
Though we focus on a magnetic mechanism for angular mo-
mentum transport, we expect the boundary layer and densities
up to ≈ 5×103 g cm−3 to be dominated by hydrodynamic in-
stabilities (which we discuss in more detail in §3).
The pressure scale height at the depth of X-ray bursts is
H ≈ 30 cm ≪ R, which allows us to use a Cartesian coordi-
nate 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 perpendic-
ular to gravity (because no work is performed), thus it is a
good approximation to consider the surface layers as concen-
tric spheres, each with constant Ω. Transfer of angular mo-
mentum is reduced to a one-dimensional diffusion equation
(Fujimoto 1993)
(R2Ω) =
, (1)
where Ω is the NS spin frequency and ν is the viscosity. The
total time derivative is given by d/dt = ∂/∂t +Vadv∂/∂z, where
Vadv is the advecting velocity of the fluid in an Eulerian frame.
In steady-state, we take ∂/∂t = 0 and set Vadv = −Ṁ/(4πR2ρ)
FIG. 1.— A diagram demonstrating the rotation profile of material from
the accretion disk, through the boundary layer, and into the NS. Material
reaches the NS surface with a Keplerian spin frequency, ΩK, the majority of
which is dissipated within the boundary layer (gray region) with thickness
HBL ≪ R. Nevertheless, angular momentum is still being added to the NS,
and a torque must be communicated through the NS. This implies a small,
non-zero, amount of shear throughout the NS interior.
so that
(R2Ω) =
. (2)
Integrating from the surface where the spin is ΩK (due to the
disk) down to a depth z where the local spin is Ω, and assum-
ing dΩ/dz = 0 at the surface,
− ṀR2ΩK + ṀR2Ω = −4πρνR4dΩ/dz. (3)
Taking the limit Ω≪ ΩK we find,
4πR3ρνqΩ = ṀR2ΩK, (4)
where q ≡ d logΩ/d logz is the shear rate. This equation
shows that q > 0 when angular momentum is transported in-
ward. In general, we will find that q is rather small (. 1) at
the depths of interest. Nevertheless q is large enough to acti-
vate 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 shear-
ing has on a magnetic field, which motivates implementation
of the Tayler-Spruit dynamo. In §4, we calculate accumulat-
ing NS envelopes without the effects of viscous angular mo-
mentum 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 insta-
bilities. 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.
TURBULENT MIXING IN ACCRETING NEUTRON STARS 3
2.1. Kelvin-Helmholtz Instability
The Kelvin-Helmholtz instability (also referred to as the dy-
namical shear instability) is governed by the Richardson num-
, (5)
where N is the Brunt-Väisälä frequency, which is composed
of contributions from thermal and compositional buoyancy,
N2 = N2T + N
µ. (6)
The thermal buoyancy is given by
N2T =
∇ad −
d lnT
d lnP
, (7)
where g = GM/R2 = 1.87× 1014 cm s−2M1.4R
6 is the surface
gravitational acceleration (ignoring the effects of general rela-
tivity), χQ = ∂ lnP/∂ lnQ, with all other intensive variables set
constant, ∇ad = (∂ lnT/∂ lnP)ad is the adiabatic temperature
gradient, the asterisk refers to derivatives of the envelope’s
profile, and the pressure scale height is,
H = P/ρg = 33.1 cm µ−11.33T8, (8)
where T8 ≡ T/10
8 K. For these analytic estimates we assume
an ideal gas equation of state and use a pure helium composi-
tion with a mean molecular weight of µ = 1.33µ1.33. We omit
the scalings with mass and radius to simplify presentation.
The compositional buoyancy is (Bildsten & Cumming 1998)
N2µ = −
d lnµe
d lnP
d lnµi
d lnP
, (9)
where µe (µi) is the mean molecular weight per electron (ion).
Linear analysis shows that Kelvin-Helmholtz instability oc-
curs when Ri < 1/4, which develops into strong turbulence
that readily transports angular momentum. This result as-
sumes 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 diffu-
sion, L2/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 ≡ VL/K (Townsend
1958). The restoring force provided by thermal buoyancy
is weakened when Pe < 1, which requires the substitution
of N2T → PeN
T and promotes instability for regions where
NT > Nµ. Thermal diffusion is most efficient at small length-
scales, which motivates setting LV/νk = Rec (Zahn 1992),
where νk is the kinematic viscosity and Rec 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 ≈ RecPr. The turbulent perturbations are
thus non-adiabatic when (Zahn 1992)
K > νkRec. (10)
In the non-degenerate surface layers the kinematic viscosity
is dominated by ions, and has a value of (Spitzer 1962)
νk = 1.4× 10
−3 cm2 s−1ρ−16 T
8 , (11)
where ρ6 ≡ ρ/10
6 g cm−3, and we assume a Coulomb loga-
rithm of lnΛ = 20. Setting K = 16σSBT
3/(3cpκρ
2), where σSB
the Stefan-Boltzmann constant, cp the specific heat, and κ the
opacity, the thermal diffusivity is
K = 48.8 cm2 s−1µ1.33κ
0.04ρ
8 , (12)
where we approximate cp = 5kB/2µmp and scale the opac-
ity to κ0.04 ≡ κ/0.04 cm
2 g−1 (the opacity is largely given by
electron scattering, but is decreased by degeneracy effects, see
Paczyński 1983; Bildsten 1998a). Substituting equations (11)
and (12) into equation (10), we find that the perturbations are
non-adiabatic at depths of ρ . 3× 107 g cm−3 T
8 . The new
“secular” Richardson number associated with this limit is,
Ris ≡
νkRec
. (13)
When Ris < 1/4, the so-called “secular shear instability”
arises.
The competing effects of accretion increasing q versus tur-
bulence developing when Ris < 1/4 (and decreasing q) drive
the surface layers toward marginally satisfying Ris = 1/4 (as-
suming 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)
= 9.2× 105 s−1µ
1.33T
8 . (14)
We substitute Ris = 1/4 into equation (13), and assuming
Rec = 1000, solve for a shear rate of
qKH = 223 κ
0.04ρ
0.1, (15)
where Ω0.1 = Ω/0.1Ω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 insta-
bilities.
2.2. Baroclinic Instability
Another important hydrodynamic instability that has been
studied extensively for accreting degenerate stars is the baro-
clinic instability (Fujimoto 1993), which we quickly sum-
marize 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 con-
stant 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 suffi-
cient 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),
RiBC ≡ 4
= 8.5× 103 µ1.33T
0.1. (16)
When Ri > RiBC, Coriolis effects limit the horizontal scale of
perturbations. This results in two parametrizations for viscos-
4 PIRO & BILDSTEN
ity estimated from linear theory (Fujimoto 1993),
νBC =
Ri1/2
H2Ω, Ri ≤ RiBC,
Ri3/2
H2Ω, Ri > RiBC,
where a factor of order unity is usually included in these pre-
scriptions, called αBC, to account for uncertainty in how lin-
ear theory relates to the saturated amplitudes of the instability.
For simplicity, we set αBC = 1 in our analysis below.
By substituting ν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 ṁ = Ṁ/(4πR2), which is typically parametrized in terms
of the local Eddington rate
ṁEdd =
(1 + X)RσTh
1.5× 105 g cm−2 s−1
(1 + X)R6
, (18)
where X is the hydrogen mass fraction and σTh is the Thomson
cross-section. The Richardson number in each case is then
2.2× 103 µ
1.33 ρ6T
0.1Ω0.1, Ri ≤ RiBC
4.4× 103 µ
1.33 ρ
0.1 Ω
0.1 , Ri > RiBC
where ṁ0.1 = ṁ/0.1ṁEdd and we have used the hydrogen
deficient (X = 0) value for ṁEdd (eq. [18]). The transi-
tion to the case Ri & RiBC occurs roughly at depths ρ &
4× 106 g cm−3 T
8 . The shear rate for each case is
qBC =
1.33ρ
0.1 Ω
0.1 , Ri ≤ RiBC
9.9 µ
1.33ρ
0.1 Ω
0.1 , Ri > RiBC
This demonstrates that generally qBC ≪ qKH, so that the baro-
clinic instability triggers before the Kelvin-Helmholtz insta-
bility. This prevents the shear rate from ever becoming large
enough for the Kelvin-Helmholtz instability to operate at
depths of ρ& 6× 104 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 insta-
bility (Wasiutyń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 in-
teresting 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φ = nBr, where n is
the number of differential turns, given as n = qΩt, and t is
the duration of the shearing. The toroidal field growth is very
fast. For the Kelvin-Helmholtz case, Bφ ∼ Br in merely ∼
10−6 s (∼ 10−4 s for the baroclinic case). As the toroidal field
becomes larger it exerts an azimuthal stress on the shearing
layer equal to BrBφ/4π, which can be written as an effective
viscosity, νe,
= ρνeqΩ. (21)
In a timescale t ≈ H2
/νe this torquing significantly decreases
the spin of the layer, where HΩ = dz/d lnΩ = R/q is the shear-
ing scale height. Setting Bφ = qΩtBr, we solve for the critical
initial radial field needed to affect the shearing,
Br,crit = (4πρ)
1/2 HΩ
= (4πρ)1/2
, (22)
so that t is basically the Alfvén travel time through a shearing
scale height. Using ρ≈ 106 g cm−3 as burning density and t ≈
1 hr, a fiducial timescale for accumulation, our estimates for
qKH and qBC imply Br,crit ∼ 10
4 G and ∼ 105 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 de-
pends on the “rotational smoothing time,” the timescale for
non-axisymmetric components of the magnetic field to be ex-
pelled as differential rotational brings field with opposite po-
larities together. This timescale is estimated as
3R2π2
ηΩ2q2
, (23)
where η is the magnetic diffusivity, which in the non-
degenerate limit is given by (Spitzer 1962)
η ≈ 0.7 cm2 s−1T
8 . (24)
If tΩ is larger than the Alfvén travel time through the layer
(given by t in eq. [22]), then there is sufficient time for mag-
netic 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Ω is sufficiently small then the shear dominates
and the magnetic field is made axisymmetric on a timescale
tΩ. Substituting tΩ into equation (22), we estimate the critical
field, below which shear dominates,
BΩ = (4πρ)
q−1/3
= 1.3× 107 G ρ
0.1 q
−1/3. (25)
The lack of persistent pulsations from accreting NSs implies
a dipole field strength . 5× 107 G (Piro & Bildsten 2005).
We therefore consider it plausible that they may have intrin-
sic magnetic fields < BΩ. In this case, we expect a very
axisymmetric magnetic field to be created in a timescale of
tΩ ∼ 200 s, which is then wrapped until it becomes dynami-
cally important.
TURBULENT MIXING IN ACCRETING NEUTRON STARS 5
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 com-
bination of these effects have been shown to give rise to the
Tayler-Spruit dynamo (Spruit 2002). In this picture, shear-
ing grows the toroidal field, which then initiates Tayler in-
stabilities (non-axisymmetric, pinch-like instabilities includ-
ing 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 deriva-
tion 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, lz, is limited by buoyancy forces to be
(eq. [6] from Spruit 2002)
lz < RωA/N, (26)
where ωA = B/[(4πρ)
1/2R] is the Alfvén frequency. At small
lengthscales, magnetic diffusion damps out perturbations. In
the limit of Ω ≫ ωA, which we are considering, the Tayler
instability growth rate is σB = ω
A/Ω, so that
l2z > η/σB = ηΩ/ω
A. (27)
Combining these two relations gives the minimum ωA needed
for the dynamo to act,
. (28)
During the timescale for Tayler instability, σ−1B , Br is stretched
into Bφ by an amount
Bφ = σ
B qΩBr. (29)
The largest amplification is achieved for magnetic fields that
extend the largest radial lengthscale available, so that assum-
ing equation (26) is marginally satisfied along with the induc-
tion equation we find Br/Bφ = lz/R = ωA/N. Combining this
with equation (29) we obtain q = (N/Ω)(ωA/Ω). Substituting
this into equation (28), we find (Spruit 2002)
qmin =
as the minimum shear needed for the Tayler instability to op-
erate. 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 typi-
cally used when invoking the Tayler-Spruit dynamo. Simula-
tions by Braithwaite (2006) demonstrate that the Tayler insta-
bility is strongest along the rotation axis, which is only real-
ized at the poles in the NS case. Since we find the dynamo to
be so much stronger than any hydrodynamic transport mech-
anisms, 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 (30) depends on
what is supplying the buoyancy. We follow Spruit (2002) and
separately consider cases of Nµ ≫ NT and Nµ ≪ NT , denoted
as case 0 and 1, respectively. In case 1 non-adiabatic effects
become important when η/K < 1, and we must take N2T →
(η/K)N2T (analogous to the above analysis of the secular shear
instability). For this case, which dominates for most of the
envelope, we find
qmin = 0.10 κ
0.04ρ
0.1 . (31)
Both qKH and qBC are considerably above this, thus the Tayler-
Spruit dynamo activates long before the onset of purely hy-
drodynamic instabilities.
By assuming that equation (28) 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
Bφ0 = (4πρ)
1/2RqΩ2/Nµ, (32)
Br0 = q
Bφ0, (33)
for case 0. For case 1,
Bφ1 = (4πρ)
1/2RΩq1/2
)1/8(
, (34)
Br1 =
)1/4(
Bφ1. (35)
The effective viscosities, as defined by equation (21), are
ν0 = R
, (36)
ν1 = R
)1/2(
. (37)
Although these viscosities are appropriate for angular mo-
mentum 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
D0 = R
, (38)
D1 = R
)3/4(
. (39)
which is just equal to the effective turbulent magnetic diffu-
sivity. 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 approxima-
tion when the viscous timescale, tvisc = H
2/ν, is less than the
timescale of accretion. For these initial estimates we focus
on the NT ≫ Nµ limit (case 1) since this dominates except at
compositional boundaries (which we revisit in §4). Substitut-
ing ν1 into equation (4), the shear rate for the Tayler-Spruit
6 PIRO & BILDSTEN
dynamo is
qTS =
)1/2(
. (40)
Scaling to values appropriate for accreting NSs,
qTS = 0.38 κ
0.04T
8 ṁ0.1Ω
0.1 . (41)
Note the scalings with ṁ and Ω. At high ṁ angular momen-
tum is fed into the star faster, creating more shear. At low
Ω 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 (com-
pare the scaling of qTS with qKH and qBC from eqs. [15]
and [20], respectively). The viscosity that gives the small-
est shear rate will likely be the most important at a given
depth. Using this criterion, we find that the Tayler-Spruit dy-
namo is dominant for densities ρ & 3 g cm−3 T
or ρ & 5 × 103 g cm−3 for T ≈ 5 × 106 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
Bφ = 1.3× 10
10 G κ
0.04 ρ
0.1 Ω
0.1 , (42)
Br = 2.1× 10
5 G κ
0.04 ρ
0.1 Ω
0.1 . (43)
Cumming & Bildsten (2000) argue that if Br & 10
6 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 accret-
ing millisecond pulsar has been explored by Lovelace et al.
(2007). At high accretion rates, Br increases and comes close
to this limit, suggesting that these magnetic fields may be
important for understanding the dynamics of X-ray burst os-
cillations, even from NSs that do not show persistent pulsa-
tions (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 destruc-
tion as the dynamo builds and destroys magnetic field
(Maeder & Meynet 2004). The heating rate per unit mass is
ν (qΩ)2 , (44)
which for the Tayler-Spruit dynamo becomes
ǫTS = 5.6× 10
10 ergs g−1 s−1κ
0.04ρ
0.1 .
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 ≈ Hρ (the column depth), and
dividing by ṁ. We then find
d lny
= 0.13 keV nuc−1µ−11.33κ
0.04T
8 ṁ0.1Ω
0.1 . (46)
In comparison, the thermal energy per nucleon at 108 K is
≈ 10 keV nuc−1 and burning of helium into carbon releases
≈ 0.6 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 im-
portance. For this, we parametrize the mixing diffusivity as
αTSD, where D is given by equations (38) or (39) and α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 increas-
ing or decreasing αTS.
The material becomes fully mixed over a scale height, H, in
a time tmix = H
2/(αTSD), which gives for case 1 (NT ≫ Nµ),
tmix = 4.3× 10
2 s α−1TSµ
1.33κ
0.04ρ
6 T8ṁ
0.1 , (47)
Accretion is also advecting material downward, which hap-
pens in a timescale tacc = y/ṁ ≈ ρH/ṁ,
tacc = 2.2× 10
3 s µ−11.33ρ6T8ṁ
0.1. (48)
The ratio of these timescales is
= 0.19 α−1TSµ
1.33κ
0.04ρ
0.1 . (49)
The scaling with ρ shows that mixing is generally more im-
portant than advection at shallow depths (ρ. 3×107 g cm−3),
while the scaling with Ω 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 angu-
lar momentum transport through the surface layers. In this
section we calculate accumulating models without directly in-
corporating 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 approxi-
mate the surface as having a constant gravitational accelera-
tion, 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 = −ρdz, giving a pressure P = gy from hy-
drostatic balance. We solve for ρ using the analytic equa-
tion of state from Paczyński (1983). For the liquid phase,
when 1 ≤ Γ≤ 173 where Γ = (4πni/3)
1/3Z2e2/(kBT ), Z is the
charge per ion, and ni 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. Previ-
ous studies of accreting NSs have shown that the in-
terior thermal balance is set by electron captures, neu-
TURBULENT MIXING IN ACCRETING NEUTRON STARS 7
tron emissions, and pycnonuclear reaction in the in-
ner crust (Miralda-Escudé et al. 1990; Zdunik et al. 1992;
Bildsten & Brown 1997; Brown & Bildsten 1998; Brown
2000, 2004) which release ≈ 1 MeV/mp ≈ 10
18 ergs g−1
(Haensel & Zdunik 1990, 2003). Depending on the accre-
tion rate and thermal structure of the crust, this energy is ei-
ther conducted into the core or released into the ocean such
that for an Eddington accretion rate up to ≈ 92% of the en-
ergy is lost to the core and exits as neutrinos (Brown 2000).
We therefore set the heating to 150 keV nuc−1, giving a flux
2.2× 1021 ergs cm−2 s−1〈ṁ〉0.1, where 〈ṁ〉0.1 is the time aver-
age accretion rate in units of 0.1ṁEdd and the average is over
timescales of order the thermal time of the crust (many years).
For simplicity we assume 〈ṁ〉 = ṁ. Recent calculations by
Gupta et al. (2006) suggest that heating is stronger than previ-
ously thought, but not sufficiently high to qualitatively change
our results. We ignore the additional flux from compressional
heating because it only contributes ∼ cpT ∼ 10 keV nuc
(Bildsten 1998a).
In a one-zone estimate ignition occurs at the base of the
helium layer when
dǫcool
, (50)
(Fujimoto et al. 1981) where ǫ3α is the heating rate from
triple-α reactions (for which we use the rate from Fushiki &
Lamb 1987 for numerical calculations),
ǫcool =
4σSBT
, (51)
and the derivatives are both taken at constant pressure. We
consider models where the column of accreted helium yacc that
is just at this ignition point. This allows the maximum amount
of time for angular momentum redistribution during the accu-
mulating 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 ra-
diative diffusion equation,
16σSBT
. (52)
The opacity is set using electron-scattering (Paczyń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 ≈ ρV 3conv
is Vconv ∼ 10
5 cm s−1, which is much less than the local sound
speed of ∼ 108 cm s−1. This means that the convection is very
efficient, and we therefore simply set (d lnT/d lnP)∗ = ∇ad
in the convective region. This convective region has not
been noted in previous studies, though a super-adiabatic tem-
perature 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 η, and N
for three accumulating models. The magnetic diffusivity is set
using the conductivity from Schatz et al. (1999). The base of
FIG. 2.— Accumulating models for accretion rates of 0.1ṁEdd (dotted line),
0.3ṁEdd (short-dashed line), and 1.0ṁEdd (solid line). For each model the
base of the helium layer is taken to be at the unstable triple-α ignition depth
where dǫ3α/dT = dǫcool/dT (thick long-dashed line in top panel). The top
panel plots the temperature, and the solid circles mark the top and bottom
of the convective zone. The middle panel plots the thermal diffusivity, K,
and magnetic diffusivity, η. The bottom panel shows the Brunt-Väisälä fre-
quency N for both thermal (NT , lines) and compositional (Nµ , solid squares)
contributions.
the accumulating helium is set where unstable helium ignition
occurs, as designated by the thick long-dashed line in the top
panel (eq. [50]). 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 NT is effec-
tively zero in this region. The change in composition at the
He/Fe boundary gives a large buoyancy contribution, which
we estimate as
∆ lnµ
, (53)
where ∆ lnµ ≈ 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 (NT ≫
Nµ). Since K ≫ η, the perturbations creating the dynamo are
non-adiabatic for case 1, as we described in §3.1. The only
place case 0 (Nµ ≫ NT ) is important is at the He/Fe boundary.
We next solve for the shear rates by solving equation (4)
with either ν0 or ν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 struc-
ture calculation.
8 PIRO & BILDSTEN
FIG. 3.— Angular momentum transport in the NS surface layers, for a spin
of Ω = 0.1ΩK , and accretion rates of 0.1ṁEdd (dotted line), 0.3ṁEdd (dashed
line), and 1.0ṁEdd (solid line). The panels display (from top to bottom) the
viscous timescale for angular momentum transport across a scale height, tvisc,
the shear rate, q, and the viscous energy deposition per logarithm column,
dE/d lny. The squares indicate the corresponding values due to the compo-
sitional discontinuity at the base of the accumulating layer.
In Figure 3 lines denoted the profiles calculated assuming
case 1 of the Tayler-Spruit dynamo. We assume that the con-
vection 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/Vconv ∼ 10
−4 s, is much less
than the timescales for mixing or accretion. The squares in-
dicate the corresponding values due to the compositional dis-
continuity found by using case 0. Since the viscous timescale,
tvisc = H
2/ν, is much less than the timescale it took to accrete
to this base column tacc = yacc/ṁ ∼ 10
3 − 105 s (for the range
of ṁ considered) our steady-state assumption is valid.
These calculations highlight the importance of the composi-
tional jump since tvisc, q, and dE/d lny are all amplified here.
This is because the large buoyancy reduces angular momen-
tum transport across this boundary. The viscous time is nearly
independent of accretion rate in regions where NT ≫ Nµ be-
cause 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 degener-
acy effects decreasing NT . New heat sources at a depth of
≈ 1012 g cm−2 could ease the difficulty calculations have in
recreating the low ignition columns needed to explain super-
burst 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
FIG. 4.— Azimuthal and radial field components for a spin of Ω = 0.1ΩK .
Lines have the same meaning as in Fig. 3.
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 hori-
zontal field drops by ≈ ṁ/0.02ṁ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
(∼ days near the top of the ocean). When the dynamo is ac-
tive, 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 Ω = 0.03, 0.1, and 0.3ΩK (67, 220, and 670 Hz),
all for ṁ = 0.1ṁEdd. Note that the shearing is most dramatic at
smaller Ω. In fact, the shearing and heating profiles are very
sensitive to the value of Ω, as was demonstrated by the ana-
lytic estimates. We do not plot the associated magnetic field
for these models since we have already plotted some exam-
ples 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
qΩd lnr = −
d lny, (54)
we approximate the spin change at this boundary as
∆Ω≈ qΩH/R. (55)
The change of spin is generally . 0.5 Hz across the accumu-
lating layer, with the majority of the spin change occurring
at the compositional boundary. The layer is very nearly in
uniform rotation.
TURBULENT MIXING IN ACCRETING NEUTRON STARS 9
FIG. 5.— Same as Fig. 3, but for spins of 0.03ΩK (solid line), 0.1ΩK
(dotted line), and 0.3ΩK (dashed line), all for ṁ = 0.1ṁEdd .
FIG. 6.— The spin frequency as a function of column depth for Ω = 0.1ΩK ,
with ṁ = 0.1 (dashed line), 0.3 (dotted line), and 1.0ṁEdd (solid line). The
dot-dashed line shows a constant spin frequency for comparison.
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 Fig-
ures 3 and 5 argue that we must take into account the large Nµ
at the He/Fe boundary. Assuming that Nµ scales like equation
(53), we estimate
Nµ = 1.6× 10
6 s−1 µ
1.33T
8 (∆ lnµ/0.44)
1/2. (56)
We substitute this into ν0 to solve for q using equation (4),
which is used to estimate a mixing timescale at the boundary,
tmix = H
2/(αTSD0)
= 1.7× 103s α−1TSµ
1.33 ρ
0.1 Ω0.1
∆ lnµ
Using equation (48) we find a ratio of
≈ 0.77 α−1TSµ
1.33 ρ
0.1 Ω0.1
∆ lnµ
Unlike in equation (49), this new ratio depends on ṁ. This is
because the viscosity in case 0 has a different dependence on
q. Now as ṁ is increased tmix decreases faster than tacc. Above
a critical accretion rate of
ṁcrit,1 = 4.6× 10
−2 ṁEdd α
1.33ρ6T
∆ lnµ
mixing can no longer be ignored. For densities and tempera-
tures expected at the base of the accumulating layer, this crit-
ical accretion rate lies in the range of 0.1 − 1.0ṁEdd. Note that
this depends very strongly on the spin rate, ṁcrit,1 ∝ Ω
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 ṁ (fixing Ω = 0.1ΩK). At low ṁ,
the He/Fe boundary at yacc (shown by the filled squares) acts
as a barrier to mixing since tmix > tacc at this depth. When
this happens it is a good approximation to ignore mixing and
assume two separate layers during the accumulation phase.
As ṁ increases, tmix at yacc becomes less and less until fi-
nally tmix < tacc, so that material should be mixed past yacc.
When this occurs, our accumulating model can no longer ig-
nore mixing. The mixing between helium and iron occurs
down to a depth where tmix is equal to the length of time ac-
cretion has been taking place, tacc = yacc/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 yacc.
5. THE EFFECTS OF TURBULENT MIXING
Once tmix < tacc 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 Fig-
ure 8 and summarize here. Material accretes at a rate ṁ for
a time tacc with a helium mass fraction Y0, supplying a col-
umn of material yacc = ṁtacc. The total column of helium that
accretes during this time is therefore Y0yacc. Mixing causes
this newly accreted material to mix past yacc to the mixing
depth ymix defined as where tmix = tacc. The helium is fully
mixed down to this depth, resulting in a diluted mass fraction
Ymix = Y0yacc/ymix 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
10 PIRO & BILDSTEN
FIG. 7.— Comparison of the mixing timescale, tmix (solid line) versus the
accretion timescale tacc (dashed line) for ṁ = 0.1,0.3, and 1.0ṁEdd (top to
bottom panel; all with Ω = 0.1ΩK and αTS = 1). The solid square denotes
tmix due to the compositional discontinuity. At low ṁ, the squares are above
the dashed line, demonstrating that tmix > tacc at the base of the accumulating
layer, which prevents mixing to larger depths. At sufficiently high ṁ, tmix <
tacc at depths below the accumulation depth, so that mixing between helium
and iron can occur.
FIG. 8.— Diagram demonstrating the main features of turbulent mixing.
Material mixes down to where tmix = tacc, which defines the depth ymix. The
total amount of material that accretes is yacc = ṁtacc , giving a total accreted
column of helium of Y0yacc. Within the mixed layer the helium mass fraction
is diluted to new mass fraction of Ymix = Y0yacc/ymix .
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 frac-
tion, Ymix < Y0, and (2) the ignition takes place at the base
of the mixed layer, at a depth ymix > yacc, resulting in a re-
currence 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 Y0 = 1. We
consider other values of Y0 for our analytic estimates in the
next section. We solve for the mixed accumulating structure
by first assuming an amount of accretion yacc, which for a
given ṁ implies an accretion time yacc/ṁ. We integrate the
radiative diffusion equation (eq. [52]) down to a depth where
tmix = tacc giving ymix, where we assume a constant flux pro-
file as was discussed in §4 since little helium burning takes
place during accumulation. We then estimate the mixed he-
lium fraction as Ymix = yacc/ymix. This estimate is improved by
substituting Ymix back into our envelope integration, and iter-
ating until Ymix 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 ṁ = 0.1ṁEdd and Ω = 0.1ΩK. In the four panels we consider
values of yacc of 10
6, 3× 106, 107, and 4× 107 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 ymix, which also gives Ymix 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
trec = 4× 10
7 g cm−2/1.5× 103 g cm−2 s−1 ≈ 44 minutes. In
contrast the model shown in Figure 2 with the same accretion
rate of ṁ = 0.1ṁEdd has trec ≈ 1.5 days! The shorter recur-
rence 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 κff ∝ 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 compari-
son to the pure-helium models considered before, which con-
tributes to the shallow ignition depths.
We calculate trec for a grid of models with various ṁ and
Ω in Figure 10. The recurrence time is shorter for stronger
mixing, which occurs at high ṁ or low Ω, and can be as short
as ≈ 5 − 30 minutes. We warn though that all of these cal-
culations assume that complete mixing can occur, and as we
already showed in §4, buoyancy may prevent this (eq. [59]).
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 ∆ lnµ may be smaller at times if, for example, there is
incomplete burning in a previous burst. At sufficiently high ṁ
or low Ω 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 α-value, which
TURBULENT MIXING IN ACCRETING NEUTRON STARS 11
FIG. 9.— The four panels show how the fully mixed accumulating layer
evolves in time until it reaches conditions necessary for unstable ignition.
The parameters of the NS are ṁ = 0.1ṁEdd and Ω = 0.1ΩK . In each panel,
the column of helium that has been accreted is denoted by a filled circle,
which is from left to right, and up to down yacc = 106 , 3 × 106, 107, and
4× 107 g cm−2. Mixing takes place down to the column reached by the thin
solid line. The mixed helium fraction, Ymix is displayed in the upper left-hand
corner of each panel. The ignition curve associated with each Ymix is shown
as a thick dashed line.
FIG. 10.— The recurrence time for mixed-ignition models as a function of
ṁ. The symbols denoted different spins, as shown in the key. Models that
are at sufficiently high ṁ or low Ω do not ignite unstably, and thus are not
plotted.
is the ratio of energy released in the persistent emission be-
tween bursts to the energy of the burst itself. For pure helium
ignition α∼ 100. In contrast, the α-values for the short recur-
rence time bursts are typically ∼ 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 α ∼ 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 α-value. This in fact may not be a problem
because incomplete burning may naturally explain the short
FIG. 11.— The mixed helium mass fraction, Ymix as a function of ṁ and Ω.
The crosses are mixed ignition models. The filled circles are for stably burn-
ing models and give the amount of helium present in the steady-state mixing
and burning layer. Steady-state burning therefore require either high ṁ or
low Ω. The scaling of Ymix ∝ ṁ
−0.21 is derived in equation (69) and is consis-
tent with these numerical results. Without mixing all of the these considered
models ignite unstably since stable accretion requires ṁ & 10ṁEdd .
recurrence times, since this would lead to smaller composi-
tional gradients and therefore stronger mixing.
We plot the Ymix as a function of ṁ and Ω 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 un-
stable 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 set-
ting Y ≈ 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 mod-
els. These solutions directly show how mixing depends on the
properties of the accumulating layer, in particular the pref-
actor αTS and Y0, 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,
κff ≈ 3.77 cm
2 g−1
, (60)
where we have estimated µe ≈ 2 as is correct within 10% for
the any of the elements of interest. Integrating the radiative
diffusion equation (eq. [52]), assuming a constant flux of F =
1021 ergs cm−2 s−1F21 and an ideal gas equation of state, the
temperature as a function of column y is
T (y) = 1.8× 108 K
µ1.33F21Z
)2/17
8 , (61)
12 PIRO & BILDSTEN
where y8 ≡ y/10
8 g cm−2. Using equation (47), we find the
mixing timescale as a function of y,
tmix = 3.5× 10
3 s α−1TSµ
−0.44
)−0.19
ṁ−10.1Ω
0.1 y
Note that tmix ∝ y
1.37 is a higher power than tacc = y/ṁ ∝ y.
This explicitly shows that mixing dominates at lower y, but
accretion always wins at some depth. Setting tmix equal to
tacc = yacc/ṁ gives the depth where mixing can extend to for a
given column of accreted material yacc,
ymix(yacc) = 1.6× 10
8 g cm−2α0.73TS µ
)0.14
−0.55
0.1 y
acc,8, (63)
where yacc,8 = yacc/10
8 g cm−2. The mixed helium fraction
down to this depth is
Ymix(yacc) =
Y0yacc
ymix(yacc)
= 0.63 α−0.73TS Y0µ
−0.32
1.33 (F21Z
2/A)−0.14Ω0.550.1 y
acc,8.
Since Ymix ∝ y
acc the strength of mixing decreases (Ymix 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-α burning is approximated as
ǫ3α = 5.3× 10
23 ergs g−1 s−1 f
, (65)
where f is factor that accounts for screening effects. To
make progress analytically we expand the exponential as
−44/T8
≈ 7.95× 10−10(T8/2.1)
21. Using our tempera-
ture profile (eq. [61]), the condition that dǫ3α/dT = dǫcool/dT
(eq. [50]) implies an ignition depth of
yign = 9.4× 10
7 g cm−2 f −0.15µ−0.571.33 F
−0.13
×(Z2/A)−0.28Y −0.44mix . (66)
Setting yign = ymix from equation (63), we solve for yacc, the
critical column of material that must be accreted to cause ig-
nition. This is then substituted back into equations (63) and
(64) to find that ignition occurs at a depth
ymix,ign = 1.2× 10
8 g cm−2 α0.38TS f
−0.13µ−0.331.33 Y
−0.38
−0.039
×(Z2/A)−0.17Ω−0.290.1 . (67)
with a composition of
Ymix,ign = 0.57 α
−0.86
−0.047µ−0.561.33 Y
−0.21
×(Z2/A)−0.25Ω0.650.1 . (68)
Comparing this with the numerical calculation is easiest if we
assume F is set by ṁ (F21 = 2.2ṁ0.1) as well as scaling Z
12 and µ ≈ 2.1 as appropriate for the iron-rich composition.
This gives
Ymix,ign ≈ 0.20 α
−0.86
TS ṁ
−0.21
0.1 Ω
0.1 . (69)
The recurrence time is trec = Ymix,ignymix,ign/ṁ, resulting in
trec ≈ 950 s α
−0.48ṁ−1.250.1 Ω
0.1 . (70)
FIG. 12.— Diagram demonstrating the main features of steady-state mix-
ing and burning. Material can mix further down to where tmix = t3α, which
defines the depth ymix. During a timescale tmix(ymix) the amount of material
which has accreted is yacc = ṁtmix(ymix). Within the mixed layer the helium
mass fraction is diluted to new mass fraction of Ymix = Y0yacc/ymix.
Equations (69) and (70) confirm the scalings found for the
numerical calculations in Figures 10 and 11. We have plotted
the Ymix ∝ ṁ
−0.21 scaling in Figure 11 to emphasize this. These
analytic results also show how strongly these results depend
on the parameterαTS, for which Ymix,ign is especially sensitive.
5.3. Steady-State Mixing and Burning
Another possibility is that the helium is mixed and burned
by triple-α 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, t3α = E3α/(Yǫ3α), where
E3α = 5.84× 10
17 ergs g−1 is the energy per mass released
from this burning. As shown in Figure 12, material is mixed
to sufficient depths where t3α (which decreases with depth)
is equal to tmix (which increases with depth), which defines
a mixing (or burning) depth ymix = y(tmix = t3α). During a
mixing timescale, the amount of material that is accreted is
yacc = ṁtmix(ymix), so that the total column of helium that has
been accreted is Y0ṁtmix(ymix). This is diluted over a depth
ymix, so that the mixed helium fraction is Ymix = Y0yacc/ymix.
This is all occurring in steady-state, material moves through
the mixed layer at a rate 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-
α burning, becomes, in the plane-parallel limit (Fujimoto
1993),
. (71)
We assume changes in ρ and D with depth are small in com-
parison to changes in Y , and, following our derivation of the
angular momentum equation in §2.2, we take the steady state
TURBULENT MIXING IN ACCRETING NEUTRON STARS 13
limit to derive
. (72)
In the limit where tacc ≫ tmix that we are interested in, the
term on the left-hand side can be dropped. Finally, making
the approximation that ρDd/dz ≈ y/tmix, we find
. (73)
This equation mimics the properties we expect from mixing.
When mixing is strong tmix/t3α ≪ 1, and dY/dy≈ 0, the com-
position does not change with depth. At the depth where
tmix ≈ t3α, dY/dy ≈ −Y/y and the helium is depleted expo-
nentially. All the helium burns into carbon, so that carbon has
a mass fraction X12 = 1 −Y .
The envelope profiles are found from simultaneously in-
tegrating three differential equations: (1) radiative transfer,
equation (52), (2) the entropy equation, dF/dy = −ǫ3α, and
(3) continuity of helium, equation (73). We integrate using
a shooting method, but first we must set three boundary con-
ditions 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 = Y0E3αṁ + Fc,
where Fc = 150 keV nuc
−1〈ṁ〉 (as discussed in §4.1). The
surface temperature is set from the radiative zero solution
(Schwarzschild 1958). The helium abundance in the mixed
region, Ymix, is an eigenvalue. It is varied until shooting gives
the correct base flux of Fc. This is easily found through itera-
tion since when Ymix is set too large, too much burning occurs,
and the base flux is too small (and vice versa for small Ymix).
In Figure 13 we plot a steady-state envelope using Y0 = 1,
ṁ = 0.3ṁEdd, Ω = 0.1ΩK, and αTS = 1. The shooting method
demonstrates that the initial helium abundance within the
mixed layer is Ymix = 0.097. The top panel shows the tem-
perature profile (thin solid line), and the critical curve for sta-
bility where dǫ3α/dT = dǫcool/dT (thick dashed line). The
bottom panel shows that the majority of triple-α burning oc-
curs at ≈ 3×108 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
tmix = t3α (as required by construction), which is deeper than
where tacc = t3α (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ǫ3α/d lnT
and dǫcool/d lnT for three different accretion rates. If the
cooling derivative is always larger, then the model is thermally
stable. We find stable accretion at ṁ’s considerably less than
the stable accretion rate of ṁ ≈ 10ṁEdd ≈ 2× 10
6 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 ṁ and small Ω are the models found to be stable.
If we where to increase αTS, a wider range of the models in
Figure 11 become stable.
5.4. The Critical 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 accre-
tion rates and spins are required for this to occur, and how
does it depend on parameters such as αTS.
FIG. 13.— An example steady-state mixing and burning envelope model
using Y0 = 1, ṁ = 1.0ṁEdd , Ω = 0.03ΩK , and αTS = 1. The material within the
mixed layer has Ymix = 0.097. The top panel shows the temperature profile
(solid line), as well as the critical curve ignition using a helium mass fraction
of 0.052 (thick dashed line, which is the helium mass fraction at the burning
depth). The middle panel compares the key timescales. The bottom panel
shows the helium (dashed line) and carbon abundances (dotted line), as well
as the energy generation rate for helium burning (thick solid line).
First we must derive the correct condition for stability in-
cluding 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
d lnT
dǫcool
d lnT
, (74)
where the derivatives are taken at constant pressure. For ǫ3α ∝
T ζρχ, where ζ = 44/T8 −3, stability requires (Bildsten 1998a)
ζ − 4 +
∂ lnκ
∂ lnT
∂ lnρ
∂ lnT
∂ lnκ
∂ lnρ
< 0. (75)
Substituting the scalings for a free-free opacity and ideal gas
equation of state we find T8 > 3.26 is required for stability.
By substituting the analytic form we found for the mixed
ignition depth (eq. [67]) into the temperature profile (eq. [61])
we can find the temperature at the base of the mixed layer,
Tign = 2.54× 10
8 K α0.09TS ṁ
0.1 Ω
−0.067
0.1 . (76)
By simply asking when Tign > 3.26× 10
8 K, we derive a sta-
bilizing ṁ of
ṁcrit,2 = 1.0ṁEddα
−0.83
0.1 . (77)
This is in reasonable agreement with Figure 11, which shows
stability occurs for 0.9ṁEdd . ṁ . 1.0ṁEdd for Ω = 0.1ΩK. It
is interesting that ṁcrit,2 is near (within an order of magnitude)
a value where bursts are observed to change. It is conceivable
14 PIRO & BILDSTEN
FIG. 14.— A comparison of dǫ3α/d lnT (solid lines) and dǫcool/d ln T
(dashed lines) as a function of depth for ṁ = 0.1, 0.3, and 1.0ṁEdd (top to
bottom panel). All the models use Ω = 0.03ΩK , αTS = 1, and Y0 = 1. This
demonstrates that only the ṁ = 1.0ṁEdd model is stable out of these three.
that this mechanism may act to stabilize X-ray bursts for the
hydrogen-rich accreting systems, and given the strong scaling
ṁcrit,2 has with αTS it is possible that more detailed calcula-
tions 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 trans-
port in the surface layers of accreting NSs. We found that the
hydrodynamic instabilities used by Fujimoto (1993) in a pre-
vious 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 heat-
ing. The turbulent mixing is sufficiently large to have impor-
tant 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 burn-
ing regimes we find in Figure 15. These can be summarized
as follows:
• Mixing is strongest at large ṁ (when angular momen-
tum is being added at greater rates) and small Ω (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
FIG. 15.— Summary of the three regimes of burning found for models in-
cluding mixing. The boundaries between each regime are shown as shaded
regions to emphasize possible uncertainty in the strength of mixing (we con-
sider αTS = 0.7 − 1.5). The heavy shaded region divides where the buoyancy
barrier is overcome (ṁcrit,1, eq. [59]), and the light shaded region divides
between stable and unstable mixed burning (ṁcrit,2, eq. [77]). An additional
uncertainty in ṁcrit,1 is the value of ∆ lnµ, which could vary depending on
the results of previous bursts.
through this, it extends down to a depth where tmix = tacc,
which is generally much deeper than the accreted col-
umn. This means that mixing should turn on abruptly
(as a function of either ṁ or Ω). It also means that the
importance of mixing depends on the particular ashes
left over from previous bursts. If, for example, incom-
plete 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 ∼ 5 − 30 minutes. Sec-
ond, 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 ≈ 0.1ṁEdd (Cornelisse et al. 2003).
The mixed ignition case can be studied easily using the cur-
rent numerical experiments (e.g., Woosley et al. 2004) by just
artificially accreting fuel with a mixed composition of Ymix ≈
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 ig-
nition of carbon in the X-ray burst ashes at columns
of ≈ 1011 − 1012 g cm−2 (Cumming & Bildsten 2001;
Strohmayer & Brown 2002). This problem could be alle-
viated by enhanced heating from the core on the order of
TURBULENT MIXING IN ACCRETING NEUTRON STARS 15
1 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 car-
bon fractions this high (Schatz et al. 2003). Observationally,
in’t Zand et al. (2003) showed that the α-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 sim-
plified compositions, since most accreting NSs are expected
to be accreting a fuel abundant in hydrogen. If the only cru-
cial burning is triple-α, such envelopes can be considered
within the framework of our models by using a solar value
of Y0 = 0.3. Our models fail, though, when sufficient carbon
is produced by triple-α to feedback into hydrogen burning
(which burns via the hot-CNO cycle, Hoyle & Fowler 1965).
Such a scenario is interesting thought because it could poten-
tially produce very hydrogen poor bursts.
In the standard theoretical framework for X-ray bursts,
flashes should be hydrogen poor at low ṁ (when there is suf-
ficient time for the hot-CNO cycle to act), and then become
mixed hydrogen-helium fuel at higher ṁ. Paradoxically, ob-
servations show a transition at a seemingly universal luminos-
ity of ≈ 2× 1037 ergs s−1 (approximately 0.1ṁEdd), but in the
opposite sense (Cornelisse et al. 2003). Cooper & Narayan
(2006b) argue that their models in fact show this transi-
tion (also see Narayan & Heyl 2003) because carbon cre-
ated during helium simmering increases the the hot-CNO
rate and the temperature, which decreases the temperature
sensitivity of triple-α reactions. To explain the discrep-
ancies 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(α,γ)19Ne
(Cooper & Narayan 2006a; Fisker et al. 2006). Unfortu-
nately, the most recent experimental results do not support
a decreased rate (Tan et al. 2007; Fisker et al. 2007). Further-
more, helium accreting systems show this same transition in
bursting properties, which is not explained within their frame-
work.
Another idea that may recreate this trend is that the frac-
tion of the star covered by the fuel increases with the global
accretion fast enough that the local accretion rate actually de-
creases (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 lu-
minosity of 2× 1037 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 spread-
ing 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 ig-
nition takes place.
Since mixing gets stronger with ṁ, it may be that the ob-
served 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 expla-
nation is that slowly spinning NSs are more likely to have
their hydrogen depleted resulting in helium-rich bursts. As-
suming that the burst oscillation frequencies are indicative
of the NS spin frequency (which is close to true or ex-
actly true for all current explanations Heyl 2004; Lee 2004;
Piro & Bildsten 2005; Payne & Melatos 2006), the systems
can be broken into slow spinning (∼ 300 Hz) and fast spin-
ning (∼ 600 Hz) classes. This slowly spinning class in-
cludes 4U 1916 − 053 (270 Hz, Galloway et al. 2001), 4U
1702 − 429 (330 Hz, Markwardt et al. 1999) and 4U 1728 − 34
(363 Hz, Strohmayer & Markwardt 1999). Do these systems
show helium-rich looking bursts, and if so, is this due to turbu-
lent 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 ≈ 50 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 α-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 dis-
cussion 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
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 ∆V . The equation for the energy per unit mass expresses the fact that
16 PIRO & BILDSTEN
energy from differential rotation can go into either viscous shearing or vertical mixing
(∆V )2
ν(qΩ)2 −
2∆ρgH
ρtmix
ν(qΩ)2 − 2N2D, (A1)
where D = H2/tmix is the mixing diffusivity. The ratio of the two right hand terms is called the flux Richardson number (Fujimoto
1988)
R f =
ν(qΩ)2
Ri, (A2)
where Ri ≡ N2/(q2Ω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 R f ∼ 0.1 − 1). From this result we can solve for the mixing diffusivity
ν. (A3)
This matches (up to a factor of order unity, R f/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.
REFERENCES
Acheson, D. J. 1978, Philos. Trans. R. Soc. Lond. A, 289, 459
Baker, N., & Kippenhahn, R. 1959, Z. Astrophys., 48, 140
Bhattacharya, D., & van den Heuvel, E.P.J., 1991, Physics Reports, 203, 1
Bildsten, L. 1995, ApJ, 438, 852
Bildsten, L. 1998a, in The Many Faces of Neutron Stars, ed. R. Buccheri, J.
van Paradijs & A. Alpar (Dordrecht: Kluwer), 419
Bildsten, L. 1998b, ApJ, 501, L89
Bildsten, L. 2000, in AIP Conf. Proc. 522, Cosmic Explosions, ed. S. S. Holt
& W. W. Zhang (New York: AIP), 359
Bildsten, L. & Brown, E. F. 1997, ApJ, 477, 897
Bildsten, L., & Cumming, A. 1998, ApJ, 506, 842
Boirin, L., Keek, L., Mendez, M., Cumming, A., in ’t Zand, J. J. M., Cottam,
J., Paerels, F., & Lewin, W. H. G. 2007, accepted for publication in
A&A(astro-ph/0701575)
Braithwaite, J. 2006, A&A, 449, 451
Brown, E. F. 2000, ApJ, 536, 915
Brown, E. F. 2004, ApJ, 614, L57
Brown, E. F. & Bildsten, L. 1998, ApJ, 496, 915
Brown, E. F., Bildsten, L., & Chang, P. 2002, ApJ, 574, 920
Chabrier, G., & Potekhin, A. Y. 1998, Phys. Rev. E, 58, 4941
Chakrabarty, D., Morgan, E. H., Muno, M. P., Galloway, D. K., Wijnands,
R., van der Klis, M., & Markwardt, C. B. 2003, Nature, 424, 42
Cooper, R. L. & Narayan, R. 2006a, ApJ, 648, L123
Cooper, R. L. & Narayan, R. 2006b, ApJ, 652, 584
Cornelisse, R., in’t Zand, J. J. M., Verbunt, F., Kuulkers, E., Heise, J., den
Hartog, P. R., Cocchi, M., Natalucci, L., Bazzano, A., & Ubertini, P. 2003,
A&A, 405, 1033
Cumming, A. 2003, ApJ, 595, 1077
Cumming, A. & Bildsten, L. 2000, ApJ, 544, 453
Cumming, A., & Bildsten, L. 2001, ApJ, 559, L127
Cumming, A., Zweibel, E., & Bildsten, L. 2001, ApJ, 557, 958
Cumming, A., Macbeth, J., in’t Zand, J. M. M., & Page, D. 2006, ApJ, 646,
Denissenkov, P. A. & Pinsonneault, M. 2007, ApJ, 655, 1157
Endal, A. S. & Sofia, S. 1978, ApJ, 220, 279
Fisker, J. L., Görres, J., Wiescher, M., & Davids, B. 2006, ApJ, 650, 332
Fisker, J. L., Tan, W., Görres, J., Wiescher, M., & Cooper, R. L. 2007,
submitted to ApJ(astro-ph/0702412)
Fricke, K. 1968, Z. Astrophys., 68, 317
Fujimoto, M. Y. 1987, A&A, 176, 53
Fujimoto, M. Y. 1988, A&A, 198, 163
Fujimoto, M. Y. 1993, ApJ, 419, 768
Fujimoto, M. Y., Hanawa, T., & Miyaji, S. 1981, ApJ, 247, 267
Fujimoto, M. Y., Sztajno, M., Lewin, W. G. H., & van Paradijs, J. 1987, ApJ,
319, 902
Fushiki, I., & Lamb, D. Q. 1987, ApJ, 317, 368
Galloway D., Chakrabarty D., Muno M., & Savov, P. 2001, ApJ, 549, L85
Galloway, D. K., Muno, M. P., Hartman, J. M., Savov, P., Psaltis, D., &
Chakrabarty, D. 2006, submitted to ApJS(astro-ph/0608259)
Goldreich, P. & Schubert, G. 1967, ApJ, 150, 571
Grindlay, J. E., Bailyn, C. D., Cohn, H., Lugger, P. M., Thorstensen, J. R., &
Wegner, G. 1988, ApJ, 334, L25
Gupta, S. Brown, E. F., Schatz, H., Möller, P, & Kratz, K.-L. 2006,
submitted to ApJ(astro-ph/0609826)
Haensel, P. & Zdunik, J. L. 1990, A&A, 227, 431
Haensel, P. & Zdunik, J. L. 2003, A&A, 404, L33
Heger, A., Cumming, A., & Woosley, S. E. 2005, submitted to
ApJ(astro-ph/0511292)
Heger, A., Langer, N., & Woosley, S. E. 2000, ApJ, 528, 368
Heyl, J. S. 2004, ApJ, 600, 939
Hoyle, R. & Fowler, W. A. 1965 in Quasi-Stellar Sources and Gravitational
Collapse, ed. I. Robinson, A. Schild, & E. L. Schucking (Chicago:
Univesity of Chicago Press)
in’t Zand, J. J. M., Kuulkers, E., Verbunt, F., Heise, J., & Cornelisse, R.
2003, A&A, 411, L487
Inogamov, N. A. & Sunyaev, R. A. 1999, Astron. Lett., 25, 269
Keek, L., in’t Zand, J. J. M., & Cumming, A. 2006, A&A, 455, 1031
Lee, U. 2004, ApJ, 600, 914
Lewin, W. H. G., van Paradijs, J., & Taam, R. E. 1995, in X-ray Binaries, ed.
W. H. G. Lewin, J. van Paradijs, E. P. J. van den Heuvel (Cambridge:
Cambridge Univ. Press), 175
Lovelace, R. V. E., Kulkarni, A. K., & Romanova, M. M. 2007, ApJ, 656,
Maeder, A. & Meynet, G. 2004, A&A, 422, 225
Markwardt, C. B., Strohmayer, T., & Swank, J. 1999, ApJ, 512, L125
Miller, M. C. 2000, ApJ, 537, 342
Miralda-Escudé, J., Paczyński, B., & Haensel, P. 1990, ApJ, 362, 572
Muno, M. P., Chakrabarty, D., Galloway, D. K., & Psaltis, D. 2002, ApJ,
580, 1048
Muno, M. P., Galloway, D. K., & Chakrabarty, D. 2004, ApJ, 608, 930
Narayan, R. & Heyl, J. S. 2003, ApJ, 599, 419
Paczyński, B. 1983, ApJ, 267, 315
Payne, D. J. B. & Melatos, A. 2006, 652, 597
Pedlosky, J. 1987, Geophysical Fluid Dynamics (New York: Springer)
Piro, A. L. & Bildsten, L. 2005, ApJ, 629, 438
Podsiadlowski, P., Rappaport, S., & Pfahl, E. D. 2002, ApJ, 565, 1107
Revnivtsev, M., Churazov, E., Gilfanov, M., & Sunyaev, R. 2001, A&A,
372, 138
Schatz, H., Bildsten, L., Cumming, A. & Ouelette, M. 2003, Nucl. Phys. A,
718, 247
Schatz, H., Bildsten, L., Cumming, A. & Wiescher, M. 1999, ApJ, 524, 1014
Schwarzschild, M. 1958, Structure and Evolution of the Stars (New York:
Dover Publications)
Spitzer, L. J. 1962, Physics of fully Ionized Gases (New York: Interscience
Publ.)
Spruit, H. C. 1999, A&A, 349, 189
Spruit, H. C. 2002, A&A, 381, 923
Spruit, H. C. 2006 (astro-ph/0607164)
Strohmayer, T. E. & Bildsten, L. 2006, in Compact Stellar X-Ray Sources,
eds. W. H. G. Lewin and M. van der Klis (Cambridge University Press)
(astro-ph/0301544)
Strohmayer, T. E., & Brown, E. F. 2002, ApJ, 566, 1045
Strohmayer, T. E., & Markwardt C. B. 1999, ApJ, 516, L81
Tan, W. P., Fisker, J. L., Görres, J., Couder, M., & Wiescher, M. 2007,
submitted to Phys. Rev. Lett.
Tassoul, J.-L. 1978, Theory of Rotating Stars (Princeton : Princeton Univ.
Press)
Tayler R.J., 1973, MNRAS, 161, 365
Townsend, A. A. 1958, J. Fluid Mech., 4, 361
van Paradijs, J., Penninx, W., & Lewin, W. H. G. 1988, MNRAS, 233, 437
http://arxiv.org/abs/astro-ph/0701575
http://arxiv.org/abs/astro-ph/0702412
http://arxiv.org/abs/astro-ph/0608259
http://arxiv.org/abs/astro-ph/0609826
http://arxiv.org/abs/astro-ph/0511292
http://arxiv.org/abs/astro-ph/0607164
http://arxiv.org/abs/astro-ph/0301544
TURBULENT MIXING IN ACCRETING NEUTRON STARS 17
von Zeipel, H. 1924a, MNRAS, 84, 684
von Zeipel, H.. 1924b, MNRAS, 84, 665
Wasiutyński, J. 1946, Astrophys. Norvegica, 4, 1
Wijnands, R. & van der Klis, M., 1998, Nature, 394, 344
Woosley, S. E., et al. 2004, ApJS, 151, 75
Yoon, S.-C. & Langer, N. 2004, A&A, 419, 623
Zahn, J.-P. 1992, A&A, 265, 115
Zdunik, J. L., Haensel, P., Paczyński, B., & Miralda-Escude, J. 1992, ApJ,
384, 129

  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)

ACCEPTED FOR PUBLICATION IN THE ASTROPHYSICAL JOURNAL
Preprint typeset using LATEX style emulateapj v. 08/22/09
TURBULENT MIXING IN THE SURFACE LAYERS OF ACCRETING NEUTRON STARS
ANTHONY L. PIRO1 AND LARS BILDSTEN
Kavli Institute for Theoretical Physics, Kohn Hall, University of California, Santa Barbara, CA 93106;
piro@kitp.ucsb.edu, bildsten@kitp.ucsb.edu
Accepted for publication in The Astrophysical Journal
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, where toroidal magnetic
field growth and Tayler instabilities balance to support a steady-state magnetic field. The radial and azimuthal
components have strengths of ∼ 105 G and ∼ 1010 G, respectively. This field provides a Maxwell stress on the
shearing surface layers, which leads to nearly uniform rotation at the depths of interest for X-ray bursts (near
densities of ≈ 106 g cm−3). 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 (when the spin is small, the
relative speed of the accreted material is larger), 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.
Subject headings: accretion, accretion disks — stars: magnetic fields — stars: neutron — X-rays: bursts —
X-rays: stars
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 ac-
cretion 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 ma-
jority of the rotational kinetic energy of the accreted material
is dissipated at low densities (∼ 10−4 − 10−1 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
1 Current address: Astronomy Department and Theoretical Astrophysics
Center, 601 Campbell Hall, University of California, Berkeley, CA 94720;
tpiro@astro.berkeley.edu
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 observa-
tions of these bursts. Though a simple limit cycle picture has
been successful in qualitatively explaining the bursts’ primary
characteristics, including their energies (1039 − 1040 ergs), re-
currence times (hours to days), and durations (∼ 10 − 100 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 suc-
cessive X-ray bursts with recurrence times as short as ∼ 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 in-
stabilities to create a steady-state magnetic field, which pro-
vides a torque on the shearing layers that is larger than any
purely hydrodynamic mechanism. The result is nearly uni-
form 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
http://arxiv.org/abs/0704.1278v1
2 PIRO & BILDSTEN
be non-negligible and in some cases may mix fresh mate-
rial with the ashes of previous bursts. The key is whether
the strong buoyancy due to the larger density of the ashes be-
low 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 unsta-
ble ignition. The timescale for ignition of such bursts is short
enough (∼ 5−30 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 tempera-
tures) 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 for-
mula that estimates what spins and accretion rates necessary
to transfer between these two regimes (eq. [77]).
1.1. The Basics of Angular Momentum Transport
Before beginning more detailed analysis, it is useful to
present some general equations for angular momentum trans-
port. 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 ex-
pected rotation profile from the accretion disk, through the
boundary layer, and into the NS interior. Material accreted
at a rate Ṁ reaches the NS surface with a nearly Keplerian
spin frequency ofΩK = (GM/R
3)1/2 = 1.4×104 s−1 M
1.4 R
where M1.4 ≡ M/1.4M⊙ and R6 ≡ R/10
6 cm, which has a ki-
netic energy per nucleon of ≈ 200 MeV nuc−1. The major-
ity of this energy is dissipated in a boundary layer of thick-
ness HBL ≪ 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 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 liq-
uid 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 ρ≈ 106 g cm−3, which is well below the boundary layer.
Though we focus on a magnetic mechanism for angular mo-
mentum transport, we expect the boundary layer and densities
up to ≈ 5×103 g cm−3 to be dominated by hydrodynamic in-
stabilities (which we discuss in more detail in §3).
The pressure scale height at the depth of X-ray bursts is
H ≈ 30 cm ≪ R, which allows us to use a Cartesian coordi-
nate 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 perpendic-
ular to gravity (because no work is performed), thus it is a
good approximation to consider the surface layers as concen-
tric spheres, each with constant Ω. Transfer of angular mo-
mentum is reduced to a one-dimensional diffusion equation
(Fujimoto 1993)
(R2Ω) =
, (1)
where Ω is the NS spin frequency and ν is the viscosity. The
total time derivative is given by d/dt = ∂/∂t +Vadv∂/∂z, where
Vadv is the advecting velocity of the fluid in an Eulerian frame.
In steady-state, we take ∂/∂t = 0 and set Vadv = −Ṁ/(4πR2ρ)
FIG. 1.— A diagram demonstrating the rotation profile of material from
the accretion disk, through the boundary layer, and into the NS. Material
reaches the NS surface with a Keplerian spin frequency, ΩK, the majority of
which is dissipated within the boundary layer (gray region) with thickness
HBL ≪ R. Nevertheless, angular momentum is still being added to the NS,
and a torque must be communicated through the NS. This implies a small,
non-zero, amount of shear throughout the NS interior.
so that
(R2Ω) =
. (2)
Integrating from the surface where the spin is ΩK (due to the
disk) down to a depth z where the local spin is Ω, and assum-
ing dΩ/dz = 0 at the surface,
− ṀR2ΩK + ṀR2Ω = −4πρνR4dΩ/dz. (3)
Taking the limit Ω≪ ΩK we find,
4πR3ρνqΩ = ṀR2ΩK, (4)
where q ≡ d logΩ/d logz is the shear rate. This equation
shows that q > 0 when angular momentum is transported in-
ward. In general, we will find that q is rather small (. 1) at
the depths of interest. Nevertheless q is large enough to acti-
vate 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 shear-
ing has on a magnetic field, which motivates implementation
of the Tayler-Spruit dynamo. In §4, we calculate accumulat-
ing NS envelopes without the effects of viscous angular mo-
mentum 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 insta-
bilities. 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.
TURBULENT MIXING IN ACCRETING NEUTRON STARS 3
2.1. Kelvin-Helmholtz Instability
The Kelvin-Helmholtz instability (also referred to as the dy-
namical shear instability) is governed by the Richardson num-
, (5)
where N is the Brunt-Väisälä frequency, which is composed
of contributions from thermal and compositional buoyancy,
N2 = N2T + N
µ. (6)
The thermal buoyancy is given by
N2T =
∇ad −
d lnT
d lnP
, (7)
where g = GM/R2 = 1.87× 1014 cm s−2M1.4R
6 is the surface
gravitational acceleration (ignoring the effects of general rela-
tivity), χQ = ∂ lnP/∂ lnQ, with all other intensive variables set
constant, ∇ad = (∂ lnT/∂ lnP)ad is the adiabatic temperature
gradient, the asterisk refers to derivatives of the envelope’s
profile, and the pressure scale height is,
H = P/ρg = 33.1 cm µ−11.33T8, (8)
where T8 ≡ T/10
8 K. For these analytic estimates we assume
an ideal gas equation of state and use a pure helium composi-
tion with a mean molecular weight of µ = 1.33µ1.33. We omit
the scalings with mass and radius to simplify presentation.
The compositional buoyancy is (Bildsten & Cumming 1998)
N2µ = −
d lnµe
d lnP
d lnµi
d lnP
, (9)
where µe (µi) is the mean molecular weight per electron (ion).
Linear analysis shows that Kelvin-Helmholtz instability oc-
curs when Ri < 1/4, which develops into strong turbulence
that readily transports angular momentum. This result as-
sumes 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 diffu-
sion, L2/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 ≡ VL/K (Townsend
1958). The restoring force provided by thermal buoyancy
is weakened when Pe < 1, which requires the substitution
of N2T → PeN
T and promotes instability for regions where
NT > Nµ. Thermal diffusion is most efficient at small length-
scales, which motivates setting LV/νk = Rec (Zahn 1992),
where νk is the kinematic viscosity and Rec 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 ≈ RecPr. The turbulent perturbations are
thus non-adiabatic when (Zahn 1992)
K > νkRec. (10)
In the non-degenerate surface layers the kinematic viscosity
is dominated by ions, and has a value of (Spitzer 1962)
νk = 1.4× 10
−3 cm2 s−1ρ−16 T
8 , (11)
where ρ6 ≡ ρ/10
6 g cm−3, and we assume a Coulomb loga-
rithm of lnΛ = 20. Setting K = 16σSBT
3/(3cpκρ
2), where σSB
the Stefan-Boltzmann constant, cp the specific heat, and κ the
opacity, the thermal diffusivity is
K = 48.8 cm2 s−1µ1.33κ
0.04ρ
8 , (12)
where we approximate cp = 5kB/2µmp and scale the opac-
ity to κ0.04 ≡ κ/0.04 cm
2 g−1 (the opacity is largely given by
electron scattering, but is decreased by degeneracy effects, see
Paczyński 1983; Bildsten 1998a). Substituting equations (11)
and (12) into equation (10), we find that the perturbations are
non-adiabatic at depths of ρ . 3× 107 g cm−3 T
8 . The new
“secular” Richardson number associated with this limit is,
Ris ≡
νkRec
. (13)
When Ris < 1/4, the so-called “secular shear instability”
arises.
The competing effects of accretion increasing q versus tur-
bulence developing when Ris < 1/4 (and decreasing q) drive
the surface layers toward marginally satisfying Ris = 1/4 (as-
suming 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)
= 9.2× 105 s−1µ
1.33T
8 . (14)
We substitute Ris = 1/4 into equation (13), and assuming
Rec = 1000, solve for a shear rate of
qKH = 223 κ
0.04ρ
0.1, (15)
where Ω0.1 = Ω/0.1Ω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 insta-
bilities.
2.2. Baroclinic Instability
Another important hydrodynamic instability that has been
studied extensively for accreting degenerate stars is the baro-
clinic instability (Fujimoto 1993), which we quickly sum-
marize 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 con-
stant 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 suffi-
cient 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),
RiBC ≡ 4
= 8.5× 103 µ1.33T
0.1. (16)
When Ri > RiBC, Coriolis effects limit the horizontal scale of
perturbations. This results in two parametrizations for viscos-
4 PIRO & BILDSTEN
ity estimated from linear theory (Fujimoto 1993),
νBC =
Ri1/2
H2Ω, Ri ≤ RiBC,
Ri3/2
H2Ω, Ri > RiBC,
where a factor of order unity is usually included in these pre-
scriptions, called αBC, to account for uncertainty in how lin-
ear theory relates to the saturated amplitudes of the instability.
For simplicity, we set αBC = 1 in our analysis below.
By substituting ν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 ṁ = Ṁ/(4πR2), which is typically parametrized in terms
of the local Eddington rate
ṁEdd =
(1 + X)RσTh
1.5× 105 g cm−2 s−1
(1 + X)R6
, (18)
where X is the hydrogen mass fraction and σTh is the Thomson
cross-section. The Richardson number in each case is then
2.2× 103 µ
1.33 ρ6T
0.1Ω0.1, Ri ≤ RiBC
4.4× 103 µ
1.33 ρ
0.1 Ω
0.1 , Ri > RiBC
where ṁ0.1 = ṁ/0.1ṁEdd and we have used the hydrogen
deficient (X = 0) value for ṁEdd (eq. [18]). The transi-
tion to the case Ri & RiBC occurs roughly at depths ρ &
4× 106 g cm−3 T
8 . The shear rate for each case is
qBC =
1.33ρ
0.1 Ω
0.1 , Ri ≤ RiBC
9.9 µ
1.33ρ
0.1 Ω
0.1 , Ri > RiBC
This demonstrates that generally qBC ≪ qKH, so that the baro-
clinic instability triggers before the Kelvin-Helmholtz insta-
bility. This prevents the shear rate from ever becoming large
enough for the Kelvin-Helmholtz instability to operate at
depths of ρ& 6× 104 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 insta-
bility (Wasiutyń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 in-
teresting 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φ = nBr, where n is
the number of differential turns, given as n = qΩt, and t is
the duration of the shearing. The toroidal field growth is very
fast. For the Kelvin-Helmholtz case, Bφ ∼ Br in merely ∼
10−6 s (∼ 10−4 s for the baroclinic case). As the toroidal field
becomes larger it exerts an azimuthal stress on the shearing
layer equal to BrBφ/4π, which can be written as an effective
viscosity, νe,
= ρνeqΩ. (21)
In a timescale t ≈ H2
/νe this torquing significantly decreases
the spin of the layer, where HΩ = dz/d lnΩ = R/q is the shear-
ing scale height. Setting Bφ = qΩtBr, we solve for the critical
initial radial field needed to affect the shearing,
Br,crit = (4πρ)
1/2 HΩ
= (4πρ)1/2
, (22)
so that t is basically the Alfvén travel time through a shearing
scale height. Using ρ≈ 106 g cm−3 as burning density and t ≈
1 hr, a fiducial timescale for accumulation, our estimates for
qKH and qBC imply Br,crit ∼ 10
4 G and ∼ 105 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 de-
pends on the “rotational smoothing time,” the timescale for
non-axisymmetric components of the magnetic field to be ex-
pelled as differential rotational brings field with opposite po-
larities together. This timescale is estimated as
3R2π2
ηΩ2q2
, (23)
where η is the magnetic diffusivity, which in the non-
degenerate limit is given by (Spitzer 1962)
η ≈ 0.7 cm2 s−1T
8 . (24)
If tΩ is larger than the Alfvén travel time through the layer
(given by t in eq. [22]), then there is sufficient time for mag-
netic 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Ω is sufficiently small then the shear dominates
and the magnetic field is made axisymmetric on a timescale
tΩ. Substituting tΩ into equation (22), we estimate the critical
field, below which shear dominates,
BΩ = (4πρ)
q−1/3
= 1.3× 107 G ρ
0.1 q
−1/3. (25)
The lack of persistent pulsations from accreting NSs implies
a dipole field strength . 5× 107 G (Piro & Bildsten 2005).
We therefore consider it plausible that they may have intrin-
sic magnetic fields < BΩ. In this case, we expect a very
axisymmetric magnetic field to be created in a timescale of
tΩ ∼ 200 s, which is then wrapped until it becomes dynami-
cally important.
TURBULENT MIXING IN ACCRETING NEUTRON STARS 5
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 com-
bination of these effects have been shown to give rise to the
Tayler-Spruit dynamo (Spruit 2002). In this picture, shear-
ing grows the toroidal field, which then initiates Tayler in-
stabilities (non-axisymmetric, pinch-like instabilities includ-
ing 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 deriva-
tion 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, lz, is limited by buoyancy forces to be
(eq. [6] from Spruit 2002)
lz < RωA/N, (26)
where ωA = B/[(4πρ)
1/2R] is the Alfvén frequency. At small
lengthscales, magnetic diffusion damps out perturbations. In
the limit of Ω ≫ ωA, which we are considering, the Tayler
instability growth rate is σB = ω
A/Ω, so that
l2z > η/σB = ηΩ/ω
A. (27)
Combining these two relations gives the minimum ωA needed
for the dynamo to act,
. (28)
During the timescale for Tayler instability, σ−1B , Br is stretched
into Bφ by an amount
Bφ = σ
B qΩBr. (29)
The largest amplification is achieved for magnetic fields that
extend the largest radial lengthscale available, so that assum-
ing equation (26) is marginally satisfied along with the induc-
tion equation we find Br/Bφ = lz/R = ωA/N. Combining this
with equation (29) we obtain q = (N/Ω)(ωA/Ω). Substituting
this into equation (28), we find (Spruit 2002)
qmin =
as the minimum shear needed for the Tayler instability to op-
erate. 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 typi-
cally used when invoking the Tayler-Spruit dynamo. Simula-
tions by Braithwaite (2006) demonstrate that the Tayler insta-
bility is strongest along the rotation axis, which is only real-
ized at the poles in the NS case. Since we find the dynamo to
be so much stronger than any hydrodynamic transport mech-
anisms, 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 (30) depends on
what is supplying the buoyancy. We follow Spruit (2002) and
separately consider cases of Nµ ≫ NT and Nµ ≪ NT , denoted
as case 0 and 1, respectively. In case 1 non-adiabatic effects
become important when η/K < 1, and we must take N2T →
(η/K)N2T (analogous to the above analysis of the secular shear
instability). For this case, which dominates for most of the
envelope, we find
qmin = 0.10 κ
0.04ρ
0.1 . (31)
Both qKH and qBC are considerably above this, thus the Tayler-
Spruit dynamo activates long before the onset of purely hy-
drodynamic instabilities.
By assuming that equation (28) 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
Bφ0 = (4πρ)
1/2RqΩ2/Nµ, (32)
Br0 = q
Bφ0, (33)
for case 0. For case 1,
Bφ1 = (4πρ)
1/2RΩq1/2
)1/8(
, (34)
Br1 =
)1/4(
Bφ1. (35)
The effective viscosities, as defined by equation (21), are
ν0 = R
, (36)
ν1 = R
)1/2(
. (37)
Although these viscosities are appropriate for angular mo-
mentum 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
D0 = R
, (38)
D1 = R
)3/4(
. (39)
which is just equal to the effective turbulent magnetic diffu-
sivity. 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 approxima-
tion when the viscous timescale, tvisc = H
2/ν, is less than the
timescale of accretion. For these initial estimates we focus
on the NT ≫ Nµ limit (case 1) since this dominates except at
compositional boundaries (which we revisit in §4). Substitut-
ing ν1 into equation (4), the shear rate for the Tayler-Spruit
6 PIRO & BILDSTEN
dynamo is
qTS =
)1/2(
. (40)
Scaling to values appropriate for accreting NSs,
qTS = 0.38 κ
0.04T
8 ṁ0.1Ω
0.1 . (41)
Note the scalings with ṁ and Ω. At high ṁ angular momen-
tum is fed into the star faster, creating more shear. At low
Ω 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 (com-
pare the scaling of qTS with qKH and qBC from eqs. [15]
and [20], respectively). The viscosity that gives the small-
est shear rate will likely be the most important at a given
depth. Using this criterion, we find that the Tayler-Spruit dy-
namo is dominant for densities ρ & 3 g cm−3 T
or ρ & 5 × 103 g cm−3 for T ≈ 5 × 106 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
Bφ = 1.3× 10
10 G κ
0.04 ρ
0.1 Ω
0.1 , (42)
Br = 2.1× 10
5 G κ
0.04 ρ
0.1 Ω
0.1 . (43)
Cumming & Bildsten (2000) argue that if Br & 10
6 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 accret-
ing millisecond pulsar has been explored by Lovelace et al.
(2007). At high accretion rates, Br increases and comes close
to this limit, suggesting that these magnetic fields may be
important for understanding the dynamics of X-ray burst os-
cillations, even from NSs that do not show persistent pulsa-
tions (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 destruc-
tion as the dynamo builds and destroys magnetic field
(Maeder & Meynet 2004). The heating rate per unit mass is
ν (qΩ)2 , (44)
which for the Tayler-Spruit dynamo becomes
ǫTS = 5.6× 10
10 ergs g−1 s−1κ
0.04ρ
0.1 .
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 ≈ Hρ (the column depth), and
dividing by ṁ. We then find
d lny
= 0.13 keV nuc−1µ−11.33κ
0.04T
8 ṁ0.1Ω
0.1 . (46)
In comparison, the thermal energy per nucleon at 108 K is
≈ 10 keV nuc−1 and burning of helium into carbon releases
≈ 0.6 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 im-
portance. For this, we parametrize the mixing diffusivity as
αTSD, where D is given by equations (38) or (39) and α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 increas-
ing or decreasing αTS.
The material becomes fully mixed over a scale height, H, in
a time tmix = H
2/(αTSD), which gives for case 1 (NT ≫ Nµ),
tmix = 4.3× 10
2 s α−1TSµ
1.33κ
0.04ρ
6 T8ṁ
0.1 , (47)
Accretion is also advecting material downward, which hap-
pens in a timescale tacc = y/ṁ ≈ ρH/ṁ,
tacc = 2.2× 10
3 s µ−11.33ρ6T8ṁ
0.1. (48)
The ratio of these timescales is
= 0.19 α−1TSµ
1.33κ
0.04ρ
0.1 . (49)
The scaling with ρ shows that mixing is generally more im-
portant than advection at shallow depths (ρ. 3×107 g cm−3),
while the scaling with Ω 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 angu-
lar momentum transport through the surface layers. In this
section we calculate accumulating models without directly in-
corporating 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 approxi-
mate the surface as having a constant gravitational accelera-
tion, 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 = −ρdz, giving a pressure P = gy from hy-
drostatic balance. We solve for ρ using the analytic equa-
tion of state from Paczyński (1983). For the liquid phase,
when 1 ≤ Γ≤ 173 where Γ = (4πni/3)
1/3Z2e2/(kBT ), Z is the
charge per ion, and ni 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. Previ-
ous studies of accreting NSs have shown that the in-
terior thermal balance is set by electron captures, neu-
TURBULENT MIXING IN ACCRETING NEUTRON STARS 7
tron emissions, and pycnonuclear reaction in the in-
ner crust (Miralda-Escudé et al. 1990; Zdunik et al. 1992;
Bildsten & Brown 1997; Brown & Bildsten 1998; Brown
2000, 2004) which release ≈ 1 MeV/mp ≈ 10
18 ergs g−1
(Haensel & Zdunik 1990, 2003). Depending on the accre-
tion rate and thermal structure of the crust, this energy is ei-
ther conducted into the core or released into the ocean such
that for an Eddington accretion rate up to ≈ 92% of the en-
ergy is lost to the core and exits as neutrinos (Brown 2000).
We therefore set the heating to 150 keV nuc−1, giving a flux
2.2× 1021 ergs cm−2 s−1〈ṁ〉0.1, where 〈ṁ〉0.1 is the time aver-
age accretion rate in units of 0.1ṁEdd and the average is over
timescales of order the thermal time of the crust (many years).
For simplicity we assume 〈ṁ〉 = ṁ. Recent calculations by
Gupta et al. (2006) suggest that heating is stronger than previ-
ously thought, but not sufficiently high to qualitatively change
our results. We ignore the additional flux from compressional
heating because it only contributes ∼ cpT ∼ 10 keV nuc
(Bildsten 1998a).
In a one-zone estimate ignition occurs at the base of the
helium layer when
dǫcool
, (50)
(Fujimoto et al. 1981) where ǫ3α is the heating rate from
triple-α reactions (for which we use the rate from Fushiki &
Lamb 1987 for numerical calculations),
ǫcool =
4σSBT
, (51)
and the derivatives are both taken at constant pressure. We
consider models where the column of accreted helium yacc that
is just at this ignition point. This allows the maximum amount
of time for angular momentum redistribution during the accu-
mulating 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 ra-
diative diffusion equation,
16σSBT
. (52)
The opacity is set using electron-scattering (Paczyń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 ≈ ρV 3conv
is Vconv ∼ 10
5 cm s−1, which is much less than the local sound
speed of ∼ 108 cm s−1. This means that the convection is very
efficient, and we therefore simply set (d lnT/d lnP)∗ = ∇ad
in the convective region. This convective region has not
been noted in previous studies, though a super-adiabatic tem-
perature 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 η, and N
for three accumulating models. The magnetic diffusivity is set
using the conductivity from Schatz et al. (1999). The base of
FIG. 2.— Accumulating models for accretion rates of 0.1ṁEdd (dotted line),
0.3ṁEdd (short-dashed line), and 1.0ṁEdd (solid line). For each model the
base of the helium layer is taken to be at the unstable triple-α ignition depth
where dǫ3α/dT = dǫcool/dT (thick long-dashed line in top panel). The top
panel plots the temperature, and the solid circles mark the top and bottom
of the convective zone. The middle panel plots the thermal diffusivity, K,
and magnetic diffusivity, η. The bottom panel shows the Brunt-Väisälä fre-
quency N for both thermal (NT , lines) and compositional (Nµ , solid squares)
contributions.
the accumulating helium is set where unstable helium ignition
occurs, as designated by the thick long-dashed line in the top
panel (eq. [50]). 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 NT is effec-
tively zero in this region. The change in composition at the
He/Fe boundary gives a large buoyancy contribution, which
we estimate as
∆ lnµ
, (53)
where ∆ lnµ ≈ 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 (NT ≫
Nµ). Since K ≫ η, the perturbations creating the dynamo are
non-adiabatic for case 1, as we described in §3.1. The only
place case 0 (Nµ ≫ NT ) is important is at the He/Fe boundary.
We next solve for the shear rates by solving equation (4)
with either ν0 or ν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 struc-
ture calculation.
8 PIRO & BILDSTEN
FIG. 3.— Angular momentum transport in the NS surface layers, for a spin
of Ω = 0.1ΩK , and accretion rates of 0.1ṁEdd (dotted line), 0.3ṁEdd (dashed
line), and 1.0ṁEdd (solid line). The panels display (from top to bottom) the
viscous timescale for angular momentum transport across a scale height, tvisc,
the shear rate, q, and the viscous energy deposition per logarithm column,
dE/d lny. The squares indicate the corresponding values due to the compo-
sitional discontinuity at the base of the accumulating layer.
In Figure 3 lines denoted the profiles calculated assuming
case 1 of the Tayler-Spruit dynamo. We assume that the con-
vection 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/Vconv ∼ 10
−4 s, is much less
than the timescales for mixing or accretion. The squares in-
dicate the corresponding values due to the compositional dis-
continuity found by using case 0. Since the viscous timescale,
tvisc = H
2/ν, is much less than the timescale it took to accrete
to this base column tacc = yacc/ṁ ∼ 10
3 − 105 s (for the range
of ṁ considered) our steady-state assumption is valid.
These calculations highlight the importance of the composi-
tional jump since tvisc, q, and dE/d lny are all amplified here.
This is because the large buoyancy reduces angular momen-
tum transport across this boundary. The viscous time is nearly
independent of accretion rate in regions where NT ≫ Nµ be-
cause 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 degener-
acy effects decreasing NT . New heat sources at a depth of
≈ 1012 g cm−2 could ease the difficulty calculations have in
recreating the low ignition columns needed to explain super-
burst 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
FIG. 4.— Azimuthal and radial field components for a spin of Ω = 0.1ΩK .
Lines have the same meaning as in Fig. 3.
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 hori-
zontal field drops by ≈ ṁ/0.02ṁ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
(∼ days near the top of the ocean). When the dynamo is ac-
tive, 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 Ω = 0.03, 0.1, and 0.3ΩK (67, 220, and 670 Hz),
all for ṁ = 0.1ṁEdd. Note that the shearing is most dramatic at
smaller Ω. In fact, the shearing and heating profiles are very
sensitive to the value of Ω, as was demonstrated by the ana-
lytic estimates. We do not plot the associated magnetic field
for these models since we have already plotted some exam-
ples 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
qΩd lnr = −
d lny, (54)
we approximate the spin change at this boundary as
∆Ω≈ qΩH/R. (55)
The change of spin is generally . 0.5 Hz across the accumu-
lating layer, with the majority of the spin change occurring
at the compositional boundary. The layer is very nearly in
uniform rotation.
TURBULENT MIXING IN ACCRETING NEUTRON STARS 9
FIG. 5.— Same as Fig. 3, but for spins of 0.03ΩK (solid line), 0.1ΩK
(dotted line), and 0.3ΩK (dashed line), all for ṁ = 0.1ṁEdd .
FIG. 6.— The spin frequency as a function of column depth for Ω = 0.1ΩK ,
with ṁ = 0.1 (dashed line), 0.3 (dotted line), and 1.0ṁEdd (solid line). The
dot-dashed line shows a constant spin frequency for comparison.
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 Fig-
ures 3 and 5 argue that we must take into account the large Nµ
at the He/Fe boundary. Assuming that Nµ scales like equation
(53), we estimate
Nµ = 1.6× 10
6 s−1 µ
1.33T
8 (∆ lnµ/0.44)
1/2. (56)
We substitute this into ν0 to solve for q using equation (4),
which is used to estimate a mixing timescale at the boundary,
tmix = H
2/(αTSD0)
= 1.7× 103s α−1TSµ
1.33 ρ
0.1 Ω0.1
∆ lnµ
Using equation (48) we find a ratio of
≈ 0.77 α−1TSµ
1.33 ρ
0.1 Ω0.1
∆ lnµ
Unlike in equation (49), this new ratio depends on ṁ. This is
because the viscosity in case 0 has a different dependence on
q. Now as ṁ is increased tmix decreases faster than tacc. Above
a critical accretion rate of
ṁcrit,1 = 4.6× 10
−2 ṁEdd α
1.33ρ6T
∆ lnµ
mixing can no longer be ignored. For densities and tempera-
tures expected at the base of the accumulating layer, this crit-
ical accretion rate lies in the range of 0.1 − 1.0ṁEdd. Note that
this depends very strongly on the spin rate, ṁcrit,1 ∝ Ω
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 ṁ (fixing Ω = 0.1ΩK). At low ṁ,
the He/Fe boundary at yacc (shown by the filled squares) acts
as a barrier to mixing since tmix > tacc at this depth. When
this happens it is a good approximation to ignore mixing and
assume two separate layers during the accumulation phase.
As ṁ increases, tmix at yacc becomes less and less until fi-
nally tmix < tacc, so that material should be mixed past yacc.
When this occurs, our accumulating model can no longer ig-
nore mixing. The mixing between helium and iron occurs
down to a depth where tmix is equal to the length of time ac-
cretion has been taking place, tacc = yacc/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 yacc.
5. THE EFFECTS OF TURBULENT MIXING
Once tmix < tacc 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 Fig-
ure 8 and summarize here. Material accretes at a rate ṁ for
a time tacc with a helium mass fraction Y0, supplying a col-
umn of material yacc = ṁtacc. The total column of helium that
accretes during this time is therefore Y0yacc. Mixing causes
this newly accreted material to mix past yacc to the mixing
depth ymix defined as where tmix = tacc. The helium is fully
mixed down to this depth, resulting in a diluted mass fraction
Ymix = Y0yacc/ymix 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
10 PIRO & BILDSTEN
FIG. 7.— Comparison of the mixing timescale, tmix (solid line) versus the
accretion timescale tacc (dashed line) for ṁ = 0.1,0.3, and 1.0ṁEdd (top to
bottom panel; all with Ω = 0.1ΩK and αTS = 1). The solid square denotes
tmix due to the compositional discontinuity. At low ṁ, the squares are above
the dashed line, demonstrating that tmix > tacc at the base of the accumulating
layer, which prevents mixing to larger depths. At sufficiently high ṁ, tmix <
tacc at depths below the accumulation depth, so that mixing between helium
and iron can occur.
FIG. 8.— Diagram demonstrating the main features of turbulent mixing.
Material mixes down to where tmix = tacc, which defines the depth ymix. The
total amount of material that accretes is yacc = ṁtacc , giving a total accreted
column of helium of Y0yacc. Within the mixed layer the helium mass fraction
is diluted to new mass fraction of Ymix = Y0yacc/ymix .
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 frac-
tion, Ymix < Y0, and (2) the ignition takes place at the base
of the mixed layer, at a depth ymix > yacc, resulting in a re-
currence 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 Y0 = 1. We
consider other values of Y0 for our analytic estimates in the
next section. We solve for the mixed accumulating structure
by first assuming an amount of accretion yacc, which for a
given ṁ implies an accretion time yacc/ṁ. We integrate the
radiative diffusion equation (eq. [52]) down to a depth where
tmix = tacc giving ymix, where we assume a constant flux pro-
file as was discussed in §4 since little helium burning takes
place during accumulation. We then estimate the mixed he-
lium fraction as Ymix = yacc/ymix. This estimate is improved by
substituting Ymix back into our envelope integration, and iter-
ating until Ymix 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 ṁ = 0.1ṁEdd and Ω = 0.1ΩK. In the four panels we consider
values of yacc of 10
6, 3× 106, 107, and 4× 107 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 ymix, which also gives Ymix 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
trec = 4× 10
7 g cm−2/1.5× 103 g cm−2 s−1 ≈ 44 minutes. In
contrast the model shown in Figure 2 with the same accretion
rate of ṁ = 0.1ṁEdd has trec ≈ 1.5 days! The shorter recur-
rence 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 κff ∝ 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 compari-
son to the pure-helium models considered before, which con-
tributes to the shallow ignition depths.
We calculate trec for a grid of models with various ṁ and
Ω in Figure 10. The recurrence time is shorter for stronger
mixing, which occurs at high ṁ or low Ω, and can be as short
as ≈ 5 − 30 minutes. We warn though that all of these cal-
culations assume that complete mixing can occur, and as we
already showed in §4, buoyancy may prevent this (eq. [59]).
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 ∆ lnµ may be smaller at times if, for example, there is
incomplete burning in a previous burst. At sufficiently high ṁ
or low Ω 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 α-value, which
TURBULENT MIXING IN ACCRETING NEUTRON STARS 11
FIG. 9.— The four panels show how the fully mixed accumulating layer
evolves in time until it reaches conditions necessary for unstable ignition.
The parameters of the NS are ṁ = 0.1ṁEdd and Ω = 0.1ΩK . In each panel,
the column of helium that has been accreted is denoted by a filled circle,
which is from left to right, and up to down yacc = 106 , 3 × 106, 107, and
4× 107 g cm−2. Mixing takes place down to the column reached by the thin
solid line. The mixed helium fraction, Ymix is displayed in the upper left-hand
corner of each panel. The ignition curve associated with each Ymix is shown
as a thick dashed line.
FIG. 10.— The recurrence time for mixed-ignition models as a function of
ṁ. The symbols denoted different spins, as shown in the key. Models that
are at sufficiently high ṁ or low Ω do not ignite unstably, and thus are not
plotted.
is the ratio of energy released in the persistent emission be-
tween bursts to the energy of the burst itself. For pure helium
ignition α∼ 100. In contrast, the α-values for the short recur-
rence time bursts are typically ∼ 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 α ∼ 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 α-value. This in fact may not be a problem
because incomplete burning may naturally explain the short
FIG. 11.— The mixed helium mass fraction, Ymix as a function of ṁ and Ω.
The crosses are mixed ignition models. The filled circles are for stably burn-
ing models and give the amount of helium present in the steady-state mixing
and burning layer. Steady-state burning therefore require either high ṁ or
low Ω. The scaling of Ymix ∝ ṁ
−0.21 is derived in equation (69) and is consis-
tent with these numerical results. Without mixing all of the these considered
models ignite unstably since stable accretion requires ṁ & 10ṁEdd .
recurrence times, since this would lead to smaller composi-
tional gradients and therefore stronger mixing.
We plot the Ymix as a function of ṁ and Ω 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 un-
stable 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 set-
ting Y ≈ 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 mod-
els. These solutions directly show how mixing depends on the
properties of the accumulating layer, in particular the pref-
actor αTS and Y0, 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,
κff ≈ 3.77 cm
2 g−1
, (60)
where we have estimated µe ≈ 2 as is correct within 10% for
the any of the elements of interest. Integrating the radiative
diffusion equation (eq. [52]), assuming a constant flux of F =
1021 ergs cm−2 s−1F21 and an ideal gas equation of state, the
temperature as a function of column y is
T (y) = 1.8× 108 K
µ1.33F21Z
)2/17
8 , (61)
12 PIRO & BILDSTEN
where y8 ≡ y/10
8 g cm−2. Using equation (47), we find the
mixing timescale as a function of y,
tmix = 3.5× 10
3 s α−1TSµ
−0.44
)−0.19
ṁ−10.1Ω
0.1 y
Note that tmix ∝ y
1.37 is a higher power than tacc = y/ṁ ∝ y.
This explicitly shows that mixing dominates at lower y, but
accretion always wins at some depth. Setting tmix equal to
tacc = yacc/ṁ gives the depth where mixing can extend to for a
given column of accreted material yacc,
ymix(yacc) = 1.6× 10
8 g cm−2α0.73TS µ
)0.14
−0.55
0.1 y
acc,8, (63)
where yacc,8 = yacc/10
8 g cm−2. The mixed helium fraction
down to this depth is
Ymix(yacc) =
Y0yacc
ymix(yacc)
= 0.63 α−0.73TS Y0µ
−0.32
1.33 (F21Z
2/A)−0.14Ω0.550.1 y
acc,8.
Since Ymix ∝ y
acc the strength of mixing decreases (Ymix 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-α burning is approximated as
ǫ3α = 5.3× 10
23 ergs g−1 s−1 f
, (65)
where f is factor that accounts for screening effects. To
make progress analytically we expand the exponential as
−44/T8
≈ 7.95× 10−10(T8/2.1)
21. Using our tempera-
ture profile (eq. [61]), the condition that dǫ3α/dT = dǫcool/dT
(eq. [50]) implies an ignition depth of
yign = 9.4× 10
7 g cm−2 f −0.15µ−0.571.33 F
−0.13
×(Z2/A)−0.28Y −0.44mix . (66)
Setting yign = ymix from equation (63), we solve for yacc, the
critical column of material that must be accreted to cause ig-
nition. This is then substituted back into equations (63) and
(64) to find that ignition occurs at a depth
ymix,ign = 1.2× 10
8 g cm−2 α0.38TS f
−0.13µ−0.331.33 Y
−0.38
−0.039
×(Z2/A)−0.17Ω−0.290.1 . (67)
with a composition of
Ymix,ign = 0.57 α
−0.86
−0.047µ−0.561.33 Y
−0.21
×(Z2/A)−0.25Ω0.650.1 . (68)
Comparing this with the numerical calculation is easiest if we
assume F is set by ṁ (F21 = 2.2ṁ0.1) as well as scaling Z
12 and µ ≈ 2.1 as appropriate for the iron-rich composition.
This gives
Ymix,ign ≈ 0.20 α
−0.86
TS ṁ
−0.21
0.1 Ω
0.1 . (69)
The recurrence time is trec = Ymix,ignymix,ign/ṁ, resulting in
trec ≈ 950 s α
−0.48ṁ−1.250.1 Ω
0.1 . (70)
FIG. 12.— Diagram demonstrating the main features of steady-state mix-
ing and burning. Material can mix further down to where tmix = t3α, which
defines the depth ymix. During a timescale tmix(ymix) the amount of material
which has accreted is yacc = ṁtmix(ymix). Within the mixed layer the helium
mass fraction is diluted to new mass fraction of Ymix = Y0yacc/ymix.
Equations (69) and (70) confirm the scalings found for the
numerical calculations in Figures 10 and 11. We have plotted
the Ymix ∝ ṁ
−0.21 scaling in Figure 11 to emphasize this. These
analytic results also show how strongly these results depend
on the parameterαTS, for which Ymix,ign is especially sensitive.
5.3. Steady-State Mixing and Burning
Another possibility is that the helium is mixed and burned
by triple-α 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, t3α = E3α/(Yǫ3α), where
E3α = 5.84× 10
17 ergs g−1 is the energy per mass released
from this burning. As shown in Figure 12, material is mixed
to sufficient depths where t3α (which decreases with depth)
is equal to tmix (which increases with depth), which defines
a mixing (or burning) depth ymix = y(tmix = t3α). During a
mixing timescale, the amount of material that is accreted is
yacc = ṁtmix(ymix), so that the total column of helium that has
been accreted is Y0ṁtmix(ymix). This is diluted over a depth
ymix, so that the mixed helium fraction is Ymix = Y0yacc/ymix.
This is all occurring in steady-state, material moves through
the mixed layer at a rate 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-
α burning, becomes, in the plane-parallel limit (Fujimoto
1993),
. (71)
We assume changes in ρ and D with depth are small in com-
parison to changes in Y , and, following our derivation of the
angular momentum equation in §2.2, we take the steady state
TURBULENT MIXING IN ACCRETING NEUTRON STARS 13
limit to derive
. (72)
In the limit where tacc ≫ tmix that we are interested in, the
term on the left-hand side can be dropped. Finally, making
the approximation that ρDd/dz ≈ y/tmix, we find
. (73)
This equation mimics the properties we expect from mixing.
When mixing is strong tmix/t3α ≪ 1, and dY/dy≈ 0, the com-
position does not change with depth. At the depth where
tmix ≈ t3α, dY/dy ≈ −Y/y and the helium is depleted expo-
nentially. All the helium burns into carbon, so that carbon has
a mass fraction X12 = 1 −Y .
The envelope profiles are found from simultaneously in-
tegrating three differential equations: (1) radiative transfer,
equation (52), (2) the entropy equation, dF/dy = −ǫ3α, and
(3) continuity of helium, equation (73). We integrate using
a shooting method, but first we must set three boundary con-
ditions 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 = Y0E3αṁ + Fc,
where Fc = 150 keV nuc
−1〈ṁ〉 (as discussed in §4.1). The
surface temperature is set from the radiative zero solution
(Schwarzschild 1958). The helium abundance in the mixed
region, Ymix, is an eigenvalue. It is varied until shooting gives
the correct base flux of Fc. This is easily found through itera-
tion since when Ymix is set too large, too much burning occurs,
and the base flux is too small (and vice versa for small Ymix).
In Figure 13 we plot a steady-state envelope using Y0 = 1,
ṁ = 0.3ṁEdd, Ω = 0.1ΩK, and αTS = 1. The shooting method
demonstrates that the initial helium abundance within the
mixed layer is Ymix = 0.097. The top panel shows the tem-
perature profile (thin solid line), and the critical curve for sta-
bility where dǫ3α/dT = dǫcool/dT (thick dashed line). The
bottom panel shows that the majority of triple-α burning oc-
curs at ≈ 3×108 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
tmix = t3α (as required by construction), which is deeper than
where tacc = t3α (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ǫ3α/d lnT
and dǫcool/d lnT for three different accretion rates. If the
cooling derivative is always larger, then the model is thermally
stable. We find stable accretion at ṁ’s considerably less than
the stable accretion rate of ṁ ≈ 10ṁEdd ≈ 2× 10
6 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 ṁ and small Ω are the models found to be stable.
If we where to increase αTS, a wider range of the models in
Figure 11 become stable.
5.4. The Critical 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 accre-
tion rates and spins are required for this to occur, and how
does it depend on parameters such as αTS.
FIG. 13.— An example steady-state mixing and burning envelope model
using Y0 = 1, ṁ = 1.0ṁEdd , Ω = 0.03ΩK , and αTS = 1. The material within the
mixed layer has Ymix = 0.097. The top panel shows the temperature profile
(solid line), as well as the critical curve ignition using a helium mass fraction
of 0.052 (thick dashed line, which is the helium mass fraction at the burning
depth). The middle panel compares the key timescales. The bottom panel
shows the helium (dashed line) and carbon abundances (dotted line), as well
as the energy generation rate for helium burning (thick solid line).
First we must derive the correct condition for stability in-
cluding 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
d lnT
dǫcool
d lnT
, (74)
where the derivatives are taken at constant pressure. For ǫ3α ∝
T ζρχ, where ζ = 44/T8 −3, stability requires (Bildsten 1998a)
ζ − 4 +
∂ lnκ
∂ lnT
∂ lnρ
∂ lnT
∂ lnκ
∂ lnρ
< 0. (75)
Substituting the scalings for a free-free opacity and ideal gas
equation of state we find T8 > 3.26 is required for stability.
By substituting the analytic form we found for the mixed
ignition depth (eq. [67]) into the temperature profile (eq. [61])
we can find the temperature at the base of the mixed layer,
Tign = 2.54× 10
8 K α0.09TS ṁ
0.1 Ω
−0.067
0.1 . (76)
By simply asking when Tign > 3.26× 10
8 K, we derive a sta-
bilizing ṁ of
ṁcrit,2 = 1.0ṁEddα
−0.83
0.1 . (77)
This is in reasonable agreement with Figure 11, which shows
stability occurs for 0.9ṁEdd . ṁ . 1.0ṁEdd for Ω = 0.1ΩK. It
is interesting that ṁcrit,2 is near (within an order of magnitude)
a value where bursts are observed to change. It is conceivable
14 PIRO & BILDSTEN
FIG. 14.— A comparison of dǫ3α/d lnT (solid lines) and dǫcool/d ln T
(dashed lines) as a function of depth for ṁ = 0.1, 0.3, and 1.0ṁEdd (top to
bottom panel). All the models use Ω = 0.03ΩK , αTS = 1, and Y0 = 1. This
demonstrates that only the ṁ = 1.0ṁEdd model is stable out of these three.
that this mechanism may act to stabilize X-ray bursts for the
hydrogen-rich accreting systems, and given the strong scaling
ṁcrit,2 has with αTS it is possible that more detailed calcula-
tions 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 trans-
port in the surface layers of accreting NSs. We found that the
hydrodynamic instabilities used by Fujimoto (1993) in a pre-
vious 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 heat-
ing. The turbulent mixing is sufficiently large to have impor-
tant 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 burn-
ing regimes we find in Figure 15. These can be summarized
as follows:
• Mixing is strongest at large ṁ (when angular momen-
tum is being added at greater rates) and small Ω (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
FIG. 15.— Summary of the three regimes of burning found for models in-
cluding mixing. The boundaries between each regime are shown as shaded
regions to emphasize possible uncertainty in the strength of mixing (we con-
sider αTS = 0.7 − 1.5). The heavy shaded region divides where the buoyancy
barrier is overcome (ṁcrit,1, eq. [59]), and the light shaded region divides
between stable and unstable mixed burning (ṁcrit,2, eq. [77]). An additional
uncertainty in ṁcrit,1 is the value of ∆ lnµ, which could vary depending on
the results of previous bursts.
through this, it extends down to a depth where tmix = tacc,
which is generally much deeper than the accreted col-
umn. This means that mixing should turn on abruptly
(as a function of either ṁ or Ω). It also means that the
importance of mixing depends on the particular ashes
left over from previous bursts. If, for example, incom-
plete 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 ∼ 5 − 30 minutes. Sec-
ond, 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 ≈ 0.1ṁEdd (Cornelisse et al. 2003).
The mixed ignition case can be studied easily using the cur-
rent numerical experiments (e.g., Woosley et al. 2004) by just
artificially accreting fuel with a mixed composition of Ymix ≈
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 ig-
nition of carbon in the X-ray burst ashes at columns
of ≈ 1011 − 1012 g cm−2 (Cumming & Bildsten 2001;
Strohmayer & Brown 2002). This problem could be alle-
viated by enhanced heating from the core on the order of
TURBULENT MIXING IN ACCRETING NEUTRON STARS 15
1 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 car-
bon fractions this high (Schatz et al. 2003). Observationally,
in’t Zand et al. (2003) showed that the α-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 sim-
plified compositions, since most accreting NSs are expected
to be accreting a fuel abundant in hydrogen. If the only cru-
cial burning is triple-α, such envelopes can be considered
within the framework of our models by using a solar value
of Y0 = 0.3. Our models fail, though, when sufficient carbon
is produced by triple-α to feedback into hydrogen burning
(which burns via the hot-CNO cycle, Hoyle & Fowler 1965).
Such a scenario is interesting thought because it could poten-
tially produce very hydrogen poor bursts.
In the standard theoretical framework for X-ray bursts,
flashes should be hydrogen poor at low ṁ (when there is suf-
ficient time for the hot-CNO cycle to act), and then become
mixed hydrogen-helium fuel at higher ṁ. Paradoxically, ob-
servations show a transition at a seemingly universal luminos-
ity of ≈ 2× 1037 ergs s−1 (approximately 0.1ṁEdd), but in the
opposite sense (Cornelisse et al. 2003). Cooper & Narayan
(2006b) argue that their models in fact show this transi-
tion (also see Narayan & Heyl 2003) because carbon cre-
ated during helium simmering increases the the hot-CNO
rate and the temperature, which decreases the temperature
sensitivity of triple-α reactions. To explain the discrep-
ancies 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(α,γ)19Ne
(Cooper & Narayan 2006a; Fisker et al. 2006). Unfortu-
nately, the most recent experimental results do not support
a decreased rate (Tan et al. 2007; Fisker et al. 2007). Further-
more, helium accreting systems show this same transition in
bursting properties, which is not explained within their frame-
work.
Another idea that may recreate this trend is that the frac-
tion of the star covered by the fuel increases with the global
accretion fast enough that the local accretion rate actually de-
creases (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 lu-
minosity of 2× 1037 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 spread-
ing 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 ig-
nition takes place.
Since mixing gets stronger with ṁ, it may be that the ob-
served 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 expla-
nation is that slowly spinning NSs are more likely to have
their hydrogen depleted resulting in helium-rich bursts. As-
suming that the burst oscillation frequencies are indicative
of the NS spin frequency (which is close to true or ex-
actly true for all current explanations Heyl 2004; Lee 2004;
Piro & Bildsten 2005; Payne & Melatos 2006), the systems
can be broken into slow spinning (∼ 300 Hz) and fast spin-
ning (∼ 600 Hz) classes. This slowly spinning class in-
cludes 4U 1916 − 053 (270 Hz, Galloway et al. 2001), 4U
1702 − 429 (330 Hz, Markwardt et al. 1999) and 4U 1728 − 34
(363 Hz, Strohmayer & Markwardt 1999). Do these systems
show helium-rich looking bursts, and if so, is this due to turbu-
lent 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 ≈ 50 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 α-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 dis-
cussion 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
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 ∆V . The equation for the energy per unit mass expresses the fact that
16 PIRO & BILDSTEN
energy from differential rotation can go into either viscous shearing or vertical mixing
(∆V )2
ν(qΩ)2 −
2∆ρgH
ρtmix
ν(qΩ)2 − 2N2D, (A1)
where D = H2/tmix is the mixing diffusivity. The ratio of the two right hand terms is called the flux Richardson number (Fujimoto
1988)
R f =
ν(qΩ)2
Ri, (A2)
where Ri ≡ N2/(q2Ω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 R f ∼ 0.1 − 1). From this result we can solve for the mixing diffusivity
ν. (A3)
This matches (up to a factor of order unity, R f/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.
REFERENCES
Acheson, D. J. 1978, Philos. Trans. R. Soc. Lond. A, 289, 459
Baker, N., & Kippenhahn, R. 1959, Z. Astrophys., 48, 140
Bhattacharya, D., & van den Heuvel, E.P.J., 1991, Physics Reports, 203, 1
Bildsten, L. 1995, ApJ, 438, 852
Bildsten, L. 1998a, in The Many Faces of Neutron Stars, ed. R. Buccheri, J.
van Paradijs & A. Alpar (Dordrecht: Kluwer), 419
Bildsten, L. 1998b, ApJ, 501, L89
Bildsten, L. 2000, in AIP Conf. Proc. 522, Cosmic Explosions, ed. S. S. Holt
& W. W. Zhang (New York: AIP), 359
Bildsten, L. & Brown, E. F. 1997, ApJ, 477, 897
Bildsten, L., & Cumming, A. 1998, ApJ, 506, 842
Boirin, L., Keek, L., Mendez, M., Cumming, A., in ’t Zand, J. J. M., Cottam,
J., Paerels, F., & Lewin, W. H. G. 2007, accepted for publication in
A&A(astro-ph/0701575)
Braithwaite, J. 2006, A&A, 449, 451
Brown, E. F. 2000, ApJ, 536, 915
Brown, E. F. 2004, ApJ, 614, L57
Brown, E. F. & Bildsten, L. 1998, ApJ, 496, 915
Brown, E. F., Bildsten, L., & Chang, P. 2002, ApJ, 574, 920
Chabrier, G., & Potekhin, A. Y. 1998, Phys. Rev. E, 58, 4941
Chakrabarty, D., Morgan, E. H., Muno, M. P., Galloway, D. K., Wijnands,
R., van der Klis, M., & Markwardt, C. B. 2003, Nature, 424, 42
Cooper, R. L. & Narayan, R. 2006a, ApJ, 648, L123
Cooper, R. L. & Narayan, R. 2006b, ApJ, 652, 584
Cornelisse, R., in’t Zand, J. J. M., Verbunt, F., Kuulkers, E., Heise, J., den
Hartog, P. R., Cocchi, M., Natalucci, L., Bazzano, A., & Ubertini, P. 2003,
A&A, 405, 1033
Cumming, A. 2003, ApJ, 595, 1077
Cumming, A. & Bildsten, L. 2000, ApJ, 544, 453
Cumming, A., & Bildsten, L. 2001, ApJ, 559, L127
Cumming, A., Zweibel, E., & Bildsten, L. 2001, ApJ, 557, 958
Cumming, A., Macbeth, J., in’t Zand, J. M. M., & Page, D. 2006, ApJ, 646,
Denissenkov, P. A. & Pinsonneault, M. 2007, ApJ, 655, 1157
Endal, A. S. & Sofia, S. 1978, ApJ, 220, 279
Fisker, J. L., Görres, J., Wiescher, M., & Davids, B. 2006, ApJ, 650, 332
Fisker, J. L., Tan, W., Görres, J., Wiescher, M., & Cooper, R. L. 2007,
submitted to ApJ(astro-ph/0702412)
Fricke, K. 1968, Z. Astrophys., 68, 317
Fujimoto, M. Y. 1987, A&A, 176, 53
Fujimoto, M. Y. 1988, A&A, 198, 163
Fujimoto, M. Y. 1993, ApJ, 419, 768
Fujimoto, M. Y., Hanawa, T., & Miyaji, S. 1981, ApJ, 247, 267
Fujimoto, M. Y., Sztajno, M., Lewin, W. G. H., & van Paradijs, J. 1987, ApJ,
319, 902
Fushiki, I., & Lamb, D. Q. 1987, ApJ, 317, 368
Galloway D., Chakrabarty D., Muno M., & Savov, P. 2001, ApJ, 549, L85
Galloway, D. K., Muno, M. P., Hartman, J. M., Savov, P., Psaltis, D., &
Chakrabarty, D. 2006, submitted to ApJS(astro-ph/0608259)
Goldreich, P. & Schubert, G. 1967, ApJ, 150, 571
Grindlay, J. E., Bailyn, C. D., Cohn, H., Lugger, P. M., Thorstensen, J. R., &
Wegner, G. 1988, ApJ, 334, L25
Gupta, S. Brown, E. F., Schatz, H., Möller, P, & Kratz, K.-L. 2006,
submitted to ApJ(astro-ph/0609826)
Haensel, P. & Zdunik, J. L. 1990, A&A, 227, 431
Haensel, P. & Zdunik, J. L. 2003, A&A, 404, L33
Heger, A., Cumming, A., & Woosley, S. E. 2005, submitted to
ApJ(astro-ph/0511292)
Heger, A., Langer, N., & Woosley, S. E. 2000, ApJ, 528, 368
Heyl, J. S. 2004, ApJ, 600, 939
Hoyle, R. & Fowler, W. A. 1965 in Quasi-Stellar Sources and Gravitational
Collapse, ed. I. Robinson, A. Schild, & E. L. Schucking (Chicago:
Univesity of Chicago Press)
in’t Zand, J. J. M., Kuulkers, E., Verbunt, F., Heise, J., & Cornelisse, R.
2003, A&A, 411, L487
Inogamov, N. A. & Sunyaev, R. A. 1999, Astron. Lett., 25, 269
Keek, L., in’t Zand, J. J. M., & Cumming, A. 2006, A&A, 455, 1031
Lee, U. 2004, ApJ, 600, 914
Lewin, W. H. G., van Paradijs, J., & Taam, R. E. 1995, in X-ray Binaries, ed.
W. H. G. Lewin, J. van Paradijs, E. P. J. van den Heuvel (Cambridge:
Cambridge Univ. Press), 175
Lovelace, R. V. E., Kulkarni, A. K., & Romanova, M. M. 2007, ApJ, 656,
Maeder, A. & Meynet, G. 2004, A&A, 422, 225
Markwardt, C. B., Strohmayer, T., & Swank, J. 1999, ApJ, 512, L125
Miller, M. C. 2000, ApJ, 537, 342
Miralda-Escudé, J., Paczyński, B., & Haensel, P. 1990, ApJ, 362, 572
Muno, M. P., Chakrabarty, D., Galloway, D. K., & Psaltis, D. 2002, ApJ,
580, 1048
Muno, M. P., Galloway, D. K., & Chakrabarty, D. 2004, ApJ, 608, 930
Narayan, R. & Heyl, J. S. 2003, ApJ, 599, 419
Paczyński, B. 1983, ApJ, 267, 315
Payne, D. J. B. & Melatos, A. 2006, 652, 597
Pedlosky, J. 1987, Geophysical Fluid Dynamics (New York: Springer)
Piro, A. L. & Bildsten, L. 2005, ApJ, 629, 438
Podsiadlowski, P., Rappaport, S., & Pfahl, E. D. 2002, ApJ, 565, 1107
Revnivtsev, M., Churazov, E., Gilfanov, M., & Sunyaev, R. 2001, A&A,
372, 138
Schatz, H., Bildsten, L., Cumming, A. & Ouelette, M. 2003, Nucl. Phys. A,
718, 247
Schatz, H., Bildsten, L., Cumming, A. & Wiescher, M. 1999, ApJ, 524, 1014
Schwarzschild, M. 1958, Structure and Evolution of the Stars (New York:
Dover Publications)
Spitzer, L. J. 1962, Physics of fully Ionized Gases (New York: Interscience
Publ.)
Spruit, H. C. 1999, A&A, 349, 189
Spruit, H. C. 2002, A&A, 381, 923
Spruit, H. C. 2006 (astro-ph/0607164)
Strohmayer, T. E. & Bildsten, L. 2006, in Compact Stellar X-Ray Sources,
eds. W. H. G. Lewin and M. van der Klis (Cambridge University Press)
(astro-ph/0301544)
Strohmayer, T. E., & Brown, E. F. 2002, ApJ, 566, 1045
Strohmayer, T. E., & Markwardt C. B. 1999, ApJ, 516, L81
Tan, W. P., Fisker, J. L., Görres, J., Couder, M., & Wiescher, M. 2007,
submitted to Phys. Rev. Lett.
Tassoul, J.-L. 1978, Theory of Rotating Stars (Princeton : Princeton Univ.
Press)
Tayler R.J., 1973, MNRAS, 161, 365
Townsend, A. A. 1958, J. Fluid Mech., 4, 361
van Paradijs, J., Penninx, W., & Lewin, W. H. G. 1988, MNRAS, 233, 437
http://arxiv.org/abs/astro-ph/0701575
http://arxiv.org/abs/astro-ph/0702412
http://arxiv.org/abs/astro-ph/0608259
http://arxiv.org/abs/astro-ph/0609826
http://arxiv.org/abs/astro-ph/0511292
http://arxiv.org/abs/astro-ph/0607164
http://arxiv.org/abs/astro-ph/0301544
TURBULENT MIXING IN ACCRETING NEUTRON STARS 17
von Zeipel, H. 1924a, MNRAS, 84, 684
von Zeipel, H.. 1924b, MNRAS, 84, 665
Wasiutyński, J. 1946, Astrophys. Norvegica, 4, 1
Wijnands, R. & van der Klis, M., 1998, Nature, 394, 344
Woosley, S. E., et al. 2004, ApJS, 151, 75
Yoon, S.-C. & Langer, N. 2004, A&A, 419, 623
Zahn, J.-P. 1992, A&A, 265, 115
Zdunik, J. L., Haensel, P., Paczyński, B., & Miralda-Escude, J. 1992, ApJ,
384, 129
