Anomalous diffusion across a tera-Gauss magnetic field in accreting neutron stars

Abstract

When mass falls on the polar regions of a neutron star in a binary X-ray source system, it tends to spread out over the entire surface. A long-standing question in research on this problem is: will the mass be anchored on the magnetic field lines and drag the field with it or is there a special mechanism that allows the mass to slip through the magnetic field lines, leading to much less distortion? As the amount of mass falling on the neutron star can actually be comparable with the neutron star mass, the question of which alternative holds is very important. We suggest an efficient mechanism that will allow the mass to slip through the lines. The mechanism is based on a strong ideal Schwarzschild (Structure and Evolution of the Stars. Princeton University Press, 1958) instability. As the instability itself is ideal, it cannot directly force the mass to slip though the lines. However, it can create a cascade of eddies whose scale extends down to a resistive scale, at the same time mixing the field lines up without breaking them. On this scale the mass can cross the lines. This instability is efficient enough that it can produce a mass flow in the plasma without growing to a large amplitude but saturates at a small one. The instability determines the mass per flux distribution of the accumulated material on different lines so that the equilibrium is marginal to the instability on every line. This makes the equilibrium unique. Thus, as the extra mass on the neutron star grows, the state of the outer shell proceeds through a sequence of unique critically unstable equilibria. In an appendix, an attempt is made to track the critical equilibria over long times.

1. Introduction

This paper investigates the role of anomalous diffusion, caused by a magnetohydrodynamic (MHD) instability, in spreading freshly accreted matter at the bottom of the massive accretion columns at the magnetic poles of ultra-luminous X-Ray (ULX) pulsars: accreting neutron stars with strong ($B \gtrsim 10^{13}\ \textrm {G}$) magnetic fields. It should help understand the data from several X-Ray spacecraft observing ULX sources in external nearby galaxies. These objects have extremely high X-ray luminosities of approximately $10^{40}\text {--}10^{41}\ \textrm {erg}\ \textrm {s}^{-1}$, exceeding, by orders of magnitude, the Eddington luminosity $10^{38}\,(M/M_\odot)\ \textrm {erg}\ \textrm {s}^{-1}$ at which the radiation pressure force on the electrons equals the gravitational attraction on protons by an object with mass $M$. The majority of observers have traditionally believed that these sources are ‘intermediate-mass’ black holes with masses of the order of hundreds or thousands solar masses. Recently, however, the NuSTAR and Swift spacecraft detected regular pulsations with periods of the order of a second from several of such objects, implying that they are neutron stars with huge accretion rates.

X-ray binary star emission is believed to arise from mass falling from an accretion disc onto a neutron star (see Basko & Sunyaev 1976). The binary consists of a neutron star and a normal star. The normal star evolves and fills its Roche lobe with plasma that flows into the accretion disc around the neutron star. This mass is driven radially inward through the disc until it comes into contact with those outermost magnetic field lines of the neutron star. Then magnetic reconnection allows the mass to transfer to these magnetic field lines.

As the contact point of this transfer is very far from the neutron star compared with its radius, the mass is transferred to the outermost field lines, which are anchored in a very small polar region about the neutron star's magnetic axis. Then, gravity pulls the mass onto the neutron star with great velocity along the field lines and onto this polar region. We assume the radius of this region is of order of several hundred metres.

The actual details of the falling material, how it is slowed down by a radiative shock, and finally reaches the neutron star surface, are discussed by Basko & Sunyaev (1976). In the present paper we assume that it arrives on the surface and accumulates at rest there. Our concern is what happens subsequently and how this added mass manages to merge with the ambient mass of the neutron star shell.

We assume that the neutron star has a very strong, solid, highly conducting crust below its surface. This crust forms a strong lower boundary to the outer shell of the neutron star. The plasma in this shell consists of protons and highly degenerate zero-temperature electrons. Its pressure is mainly due to the Fermi energy of the electrons. This plasma acts very much like an ordinary high-temperature plasma of high electrical conductivity. In the absence of instabilities it is strongly frozen onto the field lines of the strong $10^{12}$–$10^{ 13}$ Gauss magnetic field of the neutron star.

The weight of the accumulating plasma quickly becomes too great to be supported by the neutron star surface and a large portion of it sinks below its surface under the influence of the gravitational field of the neutron star, and adds a large amount of pressure to the original ambient pressure of the neutron star shell. This additional pressure, localized to the polar region, pushes the shell outward from the polar axis, and as a result distorts the magnetic field.

Over time this pressure forces the incoming plasma to spread over the entire surface, greatly distorting the magnetic field of the neutron star in a way that can be observed externally. (See figures 1 and 2.)

Figure 1. There are four regions involved in the infalling plasma. In region A, the plasma is slowed down by radiation that is being emitted by the neutron star. In region B, the plasma comes to rest or is slowed down and its accumulated weight presses down on the neutron star and it submerges below the surface. In region C, the ambient plasma is pressed down and its pressure is increased. It expands radially against the plasma on closed lines. Region D, which lies on closed lines, has no new plasma until flow across the magnetic field lines occurs.

Figure 2. The region, $\mathcal {R}$, where the instability occurs. This region is on the lower part of a line of force. The $\boldsymbol {s}$ direction perpendicular to the line of force is shown. Owing to the increase of force because of a finite gradient in $h(\psi)$, the pressure can increase in this direction and produce the instability, which always occurs on the lower part of the distorted line containing the instability.

In this paper, we explore this distorted magnetic field and the resulting redistribution of the additional plasma in the neutron star's outer shell. These depend on the flux freezing interaction between the plasma and the magnetic field. We show that when these distortions are large enough an ideal plasma instability will develop that will lead to breaking the flux-freezing constraint and allow the plasma to flow across the field lines. This cross-field flow plays an important role in the redistribution of the field lines and matter over the surface of the neutron star.

If there were no such cross-field flow, then, at any time, the amount of mass on each field line, between the crust and the surface would be equal to the amount of plasma that has fallen on it plus the ambient plasma initially on this line. For the lines on which plasma is falling, the extra mass on these lines would be equal to the amount of plasma that has flown along this line up to this time. In this paper, we refer to the inflowing plasma as simply the plasma and the ambient plasma as simply the ambient plasma.

On the other hand, in the case of no flux freezing, the mass of plasma in the outer neutron star shell will be redistributed and determined by a method we present in this paper.

Now, the distortion of the plasma and field at any time is determined by solving a certain partial differential equation called the Grad–Shafranov equation, after its relation to a similar equation in plasma physics. However, the unique solution of this equation depends on the total mass on each of the magnetic flux tubes. For this reason, it is important to show that, after the magnetic field is sufficiently stressed, a local instability develops on each field line. These instabilities cause the plasma to flow across the field lines to an extent we will determine. Thus, the equilibrium will determine the instabilities and the instabilities will determine the distribution of the plasma on the field lines, which will, in turn, determine the equilibrium.

Our theory of these instabilities is based on the similar situation that occurs in the convective region of the Sun and other stars. There, an analogous instability convects the luminosity across this region. The instability adjusts its strength to just the level necessary to provide the right convection. If the instability is too strong, too much heat flows, and if the instability is too weak, not enough heat flows. As the convective zone must carry the excess heat not carried by radiation, the instability is crucial to the structure of the Sun and stars. It turns out that the strength of this instability makes the entropy nearly adiabatic. A very similar situation holds in the neutron star's outer shell where the instability forces the magnetic type of entropy to be nearly constant, and allows the mass to flow across the magnetic field. It is this analogy which has guided our understanding of our problem.

The paper is divided into the following sections. Section 2 describes the condition for static equilibrium.

Section 3 describes the instability and the conditions for the equilibrium to be unstable to it on any line.

In § 4, we describe the nonlinear properties of the instability and describe the cascade of fluctuations of the magnetic field strength and plasma density which it produces. We present a simple physical picture as to how, when the eddies reach the resistive scale, the magnetic strength and plasma density fluctuations cause mass to cross the magnetic field lines.

In § 5, the analytic details of the evolution of the cascade is presented and the rate of cross-field flow is calculated. An example is given to show how close to marginality all the instabilities need to be for a steady state in an actual case.

Section 6 shows that the quasi-static equilibrium is a sequence of marginally stable states and how they are controlled by the plasma inflow onto the neutron star.

Section 7 describes other research on this problem.

Section 8 summarizes the results and draws the conclusion promised in the introduction as to the method of determining the distribution of the mass per flux and the distortion of the outer shell.

There are several appendices that supplement the calculations in the text. Appendix D presents a very approximate attempt to predict the long-time behaviour of the neutron star's shape distorted by the flow.

2. The equilibrium

The rate of the infalling plasma on the neutron star is slow enough that the equilibrium, under the combined pressure gradient, gravitational and magnetic forces, is essentially static. (We neglect the centrifugal force owing to the star's rotation.) At first, we consider the equilibrium only in the region near the north pole. The symmetric solution holds near the south pole.

This equilibrium is described by the so-called Grad–Shafranov equation.

Let the incoming matter and the neutron star's magnetic field be axisymmetric about its axis and have no toroidal component. At first, we restrict ourselves to times before the infalling matter has spread over a large region of the neutron star's surface, so that we are able to treat the ambient surface as planar. We employ cylindrical coordinates, $r, \theta , z$, and take the ambient outward magnetic field that existed before the infall of any matter, as uniform, of constant strength $B_0$ and parallel to the axis of symmetry. The $\boldsymbol {z}$ coordinate increases radially inward so that $z$ gives the depth below the ambient surface. Therefore, $\boldsymbol {B}_0$ is in the negative $\hat {z}$ direction (i.e. vertically upward). We assume that the gravitational field $\boldsymbol {g}$ is constant, of magnitude $g$ and is in the positive (downward) $z$ direction. The surface $z=0$ is the ambient neutron star surface, that exists before any mass has fallen. All quantities are in cgs units.

The components of the magnetic field can be expressed in terms of its flux function $\psi$ as

(2.1a,b)\begin{equation} B_r= -\frac{1}{r} \frac{\partial \psi }{\partial z} ; \quad B_z= \frac{1}{r} \frac{\partial \psi }{\partial r} . \end{equation}

Here $\psi=$ constant is a flux surface (as $\boldsymbol {B} \boldsymbol {\cdot } \boldsymbol {\nabla } \psi =0$). The function $\psi (r,z)$ is, thus, equal to $1/2 {\rm \pi}$ times the magnetic flux enclosed by the magnetic surface through $r$ and $z$.

The current density, $\boldsymbol {j}$, which is in the $\theta$ direction, is

(2.2)\begin{align} 4 {\rm \pi} j_\theta &= \frac{ \partial B_r}{\partial z}- \frac{ \partial B_z}{\partial r} \nonumber\\ &= -\frac{\varDelta^{*} \psi }{r}, \end{align}

where

(2.3)\begin{equation} \varDelta^{*} \equiv \left(\frac{\partial^2 }{\partial r^2} -\frac{1}{r} \frac{\partial }{\partial r} + \frac{\partial^2 }{\partial z^2} \right). \end{equation}

Now, for a magnetostatic equilibrium,

(2.4)\begin{equation} \boldsymbol{j} \times \boldsymbol{B} = \boldsymbol{\nabla} p -\rho \boldsymbol{g} . \end{equation}

The radial component of this equation is

(2.5)\begin{equation} - \frac{1}{4 {\rm \pi} r^2} (\varDelta^*\psi) \frac{\partial \psi }{\partial r}= \left( \frac{\partial p }{\partial r} \right)_z , \end{equation}

where $p$ is at first taken as a function of $z$ and $r$.

The neutron star pressure is mainly due to cold degenerate electrons for which $p= p(\rho)$. Non-relativistically, in cgs units $p =1.0 \times 10^{ 13} \rho ^{5/3}$, and ultrarelativistically, $p= 1.5 \times 10^{ 15} \rho ^{4/3}$. In general, we write $p=K \rho ^{\gamma _{A}}$ where $\gamma _{A}$ is the adiabatic index, 5/3 in the non-relativistic limit and 4/3 in the ultrarelativistic limit.

The component of (2.4) parallel to the magnetic field is

(2.6)\begin{equation} \frac{\textrm{d} p}{\textrm{d} \rho } \left( \frac{\partial \rho }{\partial z}\right)_{\psi} = \rho g , \end{equation}

where the subscript $\psi$ indicates the $z$ variation at constant $\psi$, i.e. taken along a fixed line $\psi$. Its non-relativistic solution is

(2.7)\begin{equation} \rho(z, \psi) = 22.6 (z+\textrm{const})^{3/2} . \end{equation}

We have taken $g= 2 \times 10^{14}\ \textrm {cm}\ \textrm {s}^{-2}$. The integration constant is to be determined on each line by the boundary condition at the surface. It should be zero when $z = -h(\psi)$ at the top of the accumulated plasma, which rises a height $h(\psi)$ above the ambient surface. Thus,

(2.8)\begin{equation} \rho(\psi,z) =22.6 [z+h(\psi)]^{3/2}\ \textrm{g}\ \textrm{cm}^{-3} , \end{equation}

and

(2.9)\begin{equation} p(z,h)=1.83 \times 10^{ 15}(z+h)^{5/2}\ \textrm{ergs}\ \textrm{cm}^{-3} . \end{equation}

Cancelling the factor $(\partial \psi /\partial r)$, the non-relativistic Grad–Shafranov becomes

(2.10)\begin{equation} - \frac{\varDelta^* \psi }{4 {\rm \pi} r^2} = 4.58 \times 10^{ 15} [z+h(\psi)]^{3/2} \frac{\textrm{d} h}{\textrm{d} \psi } . \end{equation}

In the relativistic region, the pressure is more complex because when the electron Fermi energy is greater than 2 MeV, the electrons combine with protons to make neutrons. The resulting equation of state for relativistic degenerate electrons in neutron stars is quite complicated and is discussed in chapter 3 of Shapiro & Teukolsky (1983). We do not get into the details and simply use $P_{\textrm {am}} (z)$ for the ambient non-relativistic and relativistic pressures.

For the general case, let $P_{\textrm {am}}(z)$ and $\rho _{\textrm {am}}(z)$ be the hydrostatic solution for the ambient pressure and density. Then, because the hydrostatic equation (2.6) does not involve $z$, $P_{\textrm {am}}[z+h(\psi)]$ and $\rho _{\textrm {am}}[z+h (\psi)]$ are also solutions along $B$ with the pressure and density still related by the same equation of state. This latter solution is the hydrostatic solution of a column of plasma lifted by a height $h(\psi)$ along the line of force $\psi$. Thus, we have

(2.11)\begin{equation} p(z, \psi) = P_{\textrm{am}}[z+h(\psi)] \end{equation}

is valid for a plasma that rises a height $h(\psi)$ above the neutron star surface. The general Grad–Shafranov equation (2.5) is then

(2.12)\begin{equation} -\frac{1}{4 {\rm \pi} r^2} \varDelta^* \psi = \frac{\partial P_{\textrm{am}}}{\partial z} \left[ z+h(\psi)\right] \frac{\textrm{d} h}{\textrm{d} \psi }. \end{equation}

However, the equations for the pressure and density are not valid within a distance $h(\psi)$ of the crust and the pressure and density have been extended by a Taylor expansion of $p$ and $\rho$. To lowest order, in this expansion the density is a constant equal to the ambient density at the crust, $\rho _{\textrm {crust}}$. In addition, the conductivity of the crust is so large that the magnetic lines of force are frozen in it and do not move during the distortion. The primeval mass per flux in the extension is equal to this density times the volume which equals the height times the cross-sectional area of the flux tube; therefore, $m(\psi)$, the extra mass of the fallen matter per $\textrm {d} \psi$, is

(2.13)\begin{equation} m(\psi) = \frac{h(\psi) \rho_{\textrm{crust}}}{B_0} . \end{equation}

Here $2 {\rm \pi} m(\psi)$ is the mass per flux. We see that the function $h(\psi)$ has three different interpretations. It gives the height of the plasma above the neutron star surface, the extra mass per unit flux and the right-hand side of the Grad–Shafranov equation. As it gives the distribution of mass, it should lead to a unique static solution as can be seen by an energy argument. The equilibrium state with a given distribution of mass on flux tubes is the lowest energy state with this constraint.

Thus, the determination of $h(\psi)$ is a fundamental problem. Its determination is related to our main question of whether there are any cross-field flows. The answer to this question is related to the existence of an instability on each line which has the capacity to drive these flows. A quasi-steady state depends on there being some cross-field flows at every line of force and, therefore, on there being an instability driving turbulent transport.

3. The instability

In general, every line of force in the static equilibrium may be subject to our instability. We present a simplified derivation of this instability here. A more rigorous derivation, based on the MHD energy principle, is given in appendix A.

In figure 2, consider a particular small region on the lower part of a line of force, and let us visualize it in local coordinates (see figure 3).

Figure 3. The diagram of the elementary physics of the buoyant instability. The orientation of the instability is indicated by the direction $\boldsymbol {s}$ in which the buoyancy occurs. This direction is at an angle to $\boldsymbol {g}$. The coordinates of the instability are indicated, and $\alpha$ is the angle between $-\boldsymbol {s}$ and $\boldsymbol {g}$.

The instability occupies a cylindrical region $\mathcal {R}$ on the lower part of a line of force as shown in figure 2. It extends a distance $\ell$ along this line of force. We describe the instability in the localized Cartesian coordinates, $s,y,\zeta$, where $\zeta$ is taken along the line on which the plasma is unstable, $y$ is in the $\theta$ direction and $s$ is perpendicular to $\zeta$ in the $r z$ plane and increases towards the axis. Here $g_s = g \cos \alpha$ is the component of $\boldsymbol {g}$ in the negative $\boldsymbol {s}$ direction where $\alpha$ is the angle between $-\boldsymbol {s}$ and $\boldsymbol {z}$.

We assume that the region $\mathcal {R}$ extends a sufficient distance along the magnetic line that locally the motions are two-dimensional and independent of $\zeta$. The lines that are straight and rigid will roll over each other. At first, we ignore any tension forces associated with the ends of the region.

In these coordinates, the initial ambient pressure, $p$, decreases along $\boldsymbol {s}$ as $\psi$ decreases, but after being upshifted by $h(\psi)$ the total distorted pressure may actually increase. The magnetic field strength $B$ decreases along $\boldsymbol {s}$ because the magnetic field is pushed outward by the added infalling mass.

As in the solar case, described in Schwarzschild (1958), consider a two-dimensional (cylindrical) bubble rising a small distance $\delta s$ away from the original line. Let the bubble expand its volume slightly by a factor $1 + \epsilon$. Inside the bubble, the density will change by the factor $(1 -\epsilon)$ and the plasma and the magnetic pressures will change by factors of $(1 -\epsilon) ^{\gamma _A}$ and $(1 -\epsilon)^{2}$.

The total pressure inside the bubble will be

(3.1)\begin{align} P_{\textrm{in}} =\left( p +\frac{B^2}{8 {\rm \pi} }\right) &= p_0 (1 -\epsilon)^{\gamma_A} +\frac{B_0^2}{8 {\rm \pi} } (1-\epsilon)^2 \nonumber\\ &= p_0 + \frac{B_0^2}{8 {\rm \pi} } - \left( \gamma_A p_0 + 2 \frac{B_0^2}{8 {\rm \pi} }\right) \epsilon , \end{align}

where the zero index refers to the bubble's initial pressure and magnetic field strength. The total pressure outside the bubble in the undisturbed plasma and magnetic field is

(3.2)\begin{equation} P_{\textrm{out}} =\left( p_0 +\frac{B_0^2}{8 {\rm \pi} }\right) + \left( p'_0 + \frac{B_0 B_0'}{4 {\rm \pi}} \right) \delta s , \end{equation}

where throughout the rest of this paper (except in one case), a prime will always denote an $s$ derivative.

If the bubble rises slowly enough, the total pressure inside the bubble must equal that outside it. Thus, $P_{\textrm {in}}= P_{\textrm {out}}$ and, hence,

(3.3)\begin{equation} \epsilon = \frac{-P_0'}{\varGamma P_0} \delta s . \end{equation}

where $\varGamma P_0 = \gamma _A p_0 +2 B_0^2/8 {\rm \pi}$. Thus, $\varGamma$ is the effective adiabatic index for the total pressure $P_0 = p_0+ B_0^2/8 {\rm \pi}$.

The gravitational (buoyancy) force per unit volume on the bubble in the $\boldsymbol {s}$ direction, $F_s$, is

(3.4)\begin{equation} F_s =-g_s( \rho_{\textrm{in}} -\rho_{\textrm{out}}) = -g_s [ -\epsilon \rho_0 -(\delta s) \rho_0')] = - g_s \rho \left( \frac{ P_0'}{\varGamma P_0 } -\frac{\rho'_0}{\rho_0} \right) \delta s . \end{equation}

Now, from our equation of state $p'/p = \gamma _A \rho ' /\rho$, so expanding the expression in the parentheses and dropping the zero subscripts, we obtain

(3.5)\begin{align} \frac{P'}{\varGamma P}-\frac{\rho' }{\rho } &=\frac{ P'}{\varGamma P }-\frac{p' }{\gamma_A p } \nonumber\\ &= \frac{1}{4{\rm \pi} \varGamma P \gamma_A p}[\gamma_A p ( B B' + 4{\rm \pi} p') - p' ( B^2+ 4 {\rm \pi} \gamma_{A} p)] \nonumber\\ &= \frac{1}{4 {\rm \pi} \varGamma P \gamma_A p} ( \gamma_A p B' B - p' B^2) \nonumber\\ &= \frac{B^2}{4 {\rm \pi} \varGamma P \gamma_A } \left( \frac{\gamma_A B' }{B } - \frac{ p' }{p } \right)= \frac{1}{\gamma_A (1+\gamma_A \beta/2)} \left( \frac{ \gamma_{A} B'}{B}-\frac{p'}{ p}\right) . \end{align}

As a result the force in the $\boldsymbol {s}$ direction is

(3.6)\begin{equation} F_s =\frac{g_s \rho}{\gamma_A} \frac{\varDelta }{ (1 + \gamma_A \beta/2) } \delta s =g_s \rho \frac{\varDelta }{C} \delta s , \end{equation}

where

(3.7)\begin{align} \varDelta &= \frac{p'}{p} -\gamma_A \frac{B'}{B} \nonumber\\ =& \frac{\textrm{d}}{\textrm{d}s} \ln{ \left( \frac{p}{B^{\gamma_A}}\right) } , \end{align}

and

(3.8)\begin{equation} C=\gamma_{A} (1+\gamma_{A} \beta /2). \end{equation}

As usual, $\beta \equiv 8 {\rm \pi} p/B^2$. In the ambient medium $\beta = 4.6 \times 10^{ -8} z^{2.5}$ for $B=10^{ 12} \ \textrm {G}$. The factor $C$ depends on $\beta$. When $\beta$ is small, $C$ is $\approx \gamma _{A}$, and when $\beta$ is large, $C$ is of the order of $\gamma _{A}^2 \beta /2$. We treat $C$ as of order one in our estimates.

The product $g_s \varDelta$ is the critical stability quantity. Assuming $g_s >0$, as is the case on the lower part of a line of force (see figure 2), $\varDelta$ must be positive for instability. If $g_s < 0$, as is the case on the upper part of the line, then $\varDelta <0$ is the condition. However, on this upper part of the line $\varDelta$ is almost always positive and there is no instability.

This instability was discovered by Newcomb (1961), and he derived its force in the form given in the parentheses of (3.4). We have reduced it to the form given in (3.6), which is more convenient for our purposes.

It is useful to express $\varDelta$ approximately in terms of the density gradient. In the local approximation, we assume $p' + B B'/4 {\rm \pi} \approx 0$ because the gravitational and curvature forces are small compared with the gradient forces in a thin region. In addition, $p'/p =\gamma _{A} \rho '/\rho$ for degenerate electrons.

Then it is easy to show that

(3.9)\begin{equation} \varDelta \approx C \frac{\rho'}{\rho} , \end{equation}

thus we have an instability if $\rho '$ (or $p'$) is greater than zero on the lower parts of the lines.

Further, using (2.11), we obtain

(3.10)\begin{align} p' &= \frac{\textrm{d} P_{\textrm{am}}[z + h] }{\textrm{d} s} = \frac{\textrm{d} P_{\textrm{am}}}{\textrm{d} z} \left( \frac{\textrm{d} z}{\textrm{d} s} +\frac{\textrm{d} h}{\textrm{d} s} \right) \nonumber\\ &= \frac{\textrm{d} P_{\textrm{am}}}{\textrm{d} z}(-\cos \alpha +\textrm{d}h/\textrm{d}s) , \end{align}

where $\alpha$ is the angle between the $- \boldsymbol {s}$ and $\boldsymbol {z}$. Therefore, (because $\textrm {d} P_{\textrm {am}}/\textrm {d} z$ is positive), on the lower part of the line where $g_s>0$ the plasma is unstable if and only if

(3.11)\begin{equation} \frac{\textrm{d} h}{\textrm{d}s} > \cos \alpha . \end{equation}

Thus, if $h$ changes fast enough in the $s$ direction, then one has instability. In other words, if too much mass accumulates inside the flux surface of the line, $\textrm {d}h/\textrm {d}s$ will be large enough to produce the instability as was hinted at in the introduction.

In our derivation of the stability condition, we have neglected the fact that our cylindrical bubble has ends and we have assumed that the magnetic field is only locally perturbed along a length $\ell$ and is unperturbed outside. However, there is a tension force ${\rm \pi} B_0^2(\delta s/\ell ^2)$ per unit volume. Thus, $\varDelta$ must be greater than zero for instability. An estimate shows that the critical value of $\varDelta$ for instability is

(3.12)\begin{equation} \varDelta_{\textrm{crit}} =\frac{{\rm \pi}^2 v_A^2 C}{g_s \ell^2 }. \end{equation}

Thus, the critical condition for instability is $\varDelta '= \varDelta -\varDelta _{\textrm {crit}} >0$. (This prime does not refer to a derivative.) When necessary, we can take care of the tension forces by replacing $\varDelta$ by $\varDelta '$. It turns out that in our applications $\varDelta '$ is very close to $\varDelta$ so, where appropriate in our subsequent discussions, we can replace $\varDelta '$ by $\varDelta$.

The quantity $\varDelta$ is a function of $s$ through the $s$ dependence of $p'/p$ and $B'/B$, so that the force $F_s$ will be positive over a limited region $- \xi < s < \xi$, where $\xi$ is a mixing length, whose magnitude we estimate in § 5. The force $F_s$ is the gradient of a local potential. If, for simplicity, we take the potential as parabolic, its depth is

(3.13)\begin{equation} \phi = \rho g_s \frac{\varDelta }{2 C} \xi^2 . \end{equation}

As an approximation to the nonlinear orbit of the unstable fluid element, we can take its corresponding velocity as a harmonic motion with amplitude

(3.14)\begin{equation} v_0 = \xi \sqrt{\frac{g_s \varDelta }{C}}\ \textrm{cm}\ \textrm{s}^{-1}. \end{equation}

It is useful to have some idea as to the magnitude of these quantities. With $g = 2 \times 10^{14} \ \textrm {cm}\ \textrm {s}^{-2}$ and the mixing length $\xi =10^{3}\ \textrm {cm}$, the velocity is $v_0 =1.4 \times 10^{10} \sqrt { \varDelta /C}\ \textrm {cm}\ \textrm {s}^{-1}$. We expect $\varDelta$ to be at most of order $10^{-4}\ \textrm {cm}^{-1}$ because it is the difference of the gradient of two equilibrium quantities and we take their scale heights as of order $10^{4}\ \textrm {cm}$. Thus, (in this case) $v_0$ should be at most $1.4 \times 10^{8}\ \textrm {cm}\ \textrm {s}^{-1}$. The nonlinear growth rate $\gamma _{\textrm {nl}} \equiv k_0 v_0/2$, where $k_0 =2 {\rm \pi} / \xi$, is from (3.14)

(3.15)\begin{equation} \gamma_{\textrm{nl}} =\frac{k_0 v_o}{2} ={\rm \pi} \sqrt{\frac{g_s \varDelta }{C}}, \end{equation}

thus $\gamma _{\textrm {nl}} =4.4 \times 10^{ 7} \sqrt {\varDelta /C}\ \textrm {s}^{-1}$, or at most $4.4 \times 10^{ 5}\ \textrm {s}^{-1}$.

When the instability becomes nonlinear, we expect there to be many such independent unstable modes in the unstable $k$ region. Treating the modes as independent, and adding all of these modes, we can treat the resulting velocity field statistically.

Expand $\boldsymbol {v}$ in a two-dimensional Fourier series in $\boldsymbol {k}$ and a Fourier integral over $\omega$,

(3.16)\begin{equation} \boldsymbol{v} = \varSigma_{\boldsymbol{k}} \int\ \textrm{d} \omega \boldsymbol{v}_{\boldsymbol{k}\omega } \exp({\textrm{i} \boldsymbol{k} \boldsymbol{\cdot} \boldsymbol{r} - \textrm{i}\omega t}). \end{equation}

We assume that $\boldsymbol {v}$ is perpendicular to $\boldsymbol {B}_0$ and independent of $\zeta$, the coordinate along $\boldsymbol {B}_0$.

Then, we use the random phase approximation to write

(3.17)\begin{equation} \langle \boldsymbol{v}_{\boldsymbol{k'}\omega'}^{*} \boldsymbol{v}_{\boldsymbol{k}\omega }\rangle = J(k,\omega) \delta_{\boldsymbol{k'} \boldsymbol{k}} \delta (\omega' - \omega) (\boldsymbol{I} - \hat{\boldsymbol{k}} \hat{\boldsymbol{k} } ) . \end{equation}

Here, $J(k, \omega)$ is the turbulent velocity spectrum that we assume axisymmetric about the magnetic field direction.

Although one expects nonlinear effects to spread the unstable velocities over many different wave length scales, we make the simplifying assumption that all the velocity modes remain on the scale of the inverse mixing length, $k_o= 2 {\rm \pi} / \xi = 6.7 \times 10^{-3}\ \textrm {cm}^{-1}$. However, the frequency of the velocity modes is not monotonic.

4. The cascade of the fluctuations of the magnetic field strength

Using these facts about the instability, we are now in a position to discuss the nonlinear cascade developed by the instability, when it saturates and how it can move plasma across the very strong neutron star magnetic field lines.

We assume that the turbulent velocities are localized to the cylindrical region $\mathcal {R}$ of radius $\xi$ and length $\ell$ along a given field line. Outside this region in the longitudinal direction along $\boldsymbol {B}_0$, or in the radial direction across $\boldsymbol {B}_0$, the plasma will be stable and there will be no turbulent velocity (see figure 2). Owing to the great strength of the magnetic field, its lines of force will be straight and rigid during the perturbations. However, as they are displaced, the magnetic strength at a fixed point varies.

Inside $\mathcal {R}$, the turbulent velocities move the straight magnetic field lines around, mixing them up. These lines will remain tied to the static outside lines so this line tying will lead to tension that will try to keep the lines in $\mathcal {R}$ from moving. As long as the instability persists owing to a non-zero (positive) $\varDelta ^{*}$, the unstable forces will be able to overcome the tension forces and mixing will continue.

Even though we have assumed that the transverse turbulent velocities remain on the mixing length scale $k_0$, these velocities will have a shear that can generate magnetic and density perturbations at much smaller scales. As long as the plasma is ideal, these perturbations will extend down to the smallest conceivable scales. However, at sufficiently small scales, resistive effects become important and truncate the magnetic and density spectra. On these scales, the plasma will no longer be ideally tied to individual lines, and resistivity can lead to transfer of mass between different lines.

The combination of these two effects, namely the propagation of the perturbations down to resistive scales and the transfer of mass to different lines on these scales, leads to systematic motions of the plasma across the lines. With time, these motions flatten the density gradient, which was driving the instability and it will terminate. Finally, the lines in $\mathcal {R}$ will have returned to their original position under the influence of the tension force tying the lines to the external region, but the plasma will not have returned to its original position. On the contrary, some of the plasma will have moved downstream in the negative $s$ direction. In this way, flux freezing will be broken, there will be a bulk motion of the plasma across the lines and the accumulated infalling mass will be able to spread over the neutron star without completely dragging the field lines with it.

The rest of the paper is devoted to deriving the rate at which this all happens.

Before proceeding to an analytical treatment of the evolution of the fluctuation spectrum let us pause and understand in a general way what is going on. We do this by making use of a simple model (see figures 4 and 5).

Figure 4. Side view of the model. The region $\mathcal {R}$ is indicated in which the rods exist. Their size is the resistive size. Two particular rods $\phi _1$ and $\phi _2$ are shown. They are being mixed up by the eddies produced by the turbulent velocity shear.

Figure 5. End-on view of the model. $(a)$ The two rods $\phi _1$ and $\phi _2$ with different densities are initially some distance apart. $(b)$ After some time they are brought close together and resistive effects will cause the plasma to flow between the rods and equalize their densities. A repetition of these collisions will equalize the densities of all the rods. $(c)$ Plot of how the density would vanish along a cross cut of the region.

Imagine the magnetic field strength and density having fluctuations at only a single scale, say the scale $k_{\eta }$. Ignoring the region along $\boldsymbol {B}_0$ outside of $\mathcal {R}$, which is stable, let us replace the magnetic field by a large collection of rigid rods oriented in the $\boldsymbol {B}_0$ direction with radii of the order of this scale (see figure 4). Let each of the rods, which we label by $\phi$, internally contain a uniform plasma and magnetic field, whose density and strength can be different for different rods. Let the rods be moved around randomly by the shear of large-scale turbulent velocities, which are perpendicular to the rods and uniform in $\ell$ (i.e. rolls).

As the rods move, let them remain rigid. Initially, let the different plasmas in the rods vary in space in a way consistent with the unstable equilibrium. Let the mass in each individual rod remain constant as the rods move.

In figure 5, we display an end view of the rods. At first, neighbouring rods will have nearly the same densities. After the rods are mixed up by transverse motions, two rods $\phi _1$ and $\phi _2$ in figure 5(a) with quite different densities and field strengths may come close together as in figure 5(b). When this happens a surface current $\boldsymbol {j}$ develops between the two rods. From Ohm's law in the frame of motion $\boldsymbol {v} \times \boldsymbol {B} = \eta \boldsymbol {j}$, this current will give rise to diffusional velocity of plasma between the two rods. It is easy to see that if the rods both have the same total pressure $p + B^2/8 {\rm \pi}$ the diffusion is from the high-density to the low-density rods and this equalizes their two densities.

If the mixing continues long enough, all the neighbouring pairs in the same neighbourhood will approach each other and all the densities in the entire cylindrical region will equalize. As a consequence the unstable density gradient $\rho ' = \rho \varDelta /C$ driving the instability will be reduced to zero, and the turbulent velocity driving the mixing will cease.

However, when one compares the original density, which had a finite gradient, with the final density, which has a stable gradient, we see that the plasma will have been shifted relative to the magnetic field lines, see figure 6. Thus, when the turbulent cascade reaches the resistive scale, it is truncated by the resistive effects described previously, which flatten the magnetic and density gradients and damp the instability.

Figure 6. An explanation of why some mass can cross magnetic field lines. A plot of the density before and after the instability. Before the instability develops, there is an unstable gradient of the density in the $s$ direction, which produces a cascade (here, $s$ increases to the left). After the cascade saturates, the density gradient vanishes and the mass distribution becomes different, so that there is less mass upstream (positive $s$) and more mass downstream, violating flux freezing. The density before is indicated by lines slanting to the left and after lines slanted to the right.

We see that the cascade proceeds in two stages. In the first stage, the fluctuations of $B$ and $\rho$ cascade down the scales to a size where resistive effects come into play. We show in the next section that this happens at a rate of order $\gamma _{\textrm {nl}}/16$, see also appendix B. In the second stage, the resistive effects will flatten the density spectrum at the resistive diffusion rate $k^2 \eta$ (see appendix C). At large scales, this diffusion rate is small compared with the cascading rate, but at a certain scale $k_{\textrm {eq}}$, it will balance the rate of flow of the first stage so as the fluctuations propagate down the cascade they will be damped by the second stage effects. This scale will be roughly equal to $k_{\eta }$ where

(4.1)\begin{equation} k_{\eta}^2 \eta = \gamma_{\textrm{nl}} . \end{equation}

At still shorter wave lengths there will be no fluctuations.

We will separate these two stages. In the first stage, treated in the next section, we will derive the evolution of the magnetic strength fluctuations as though there were no resistivity. In the second stage treated in appendix C we will derive the resistive damping for $k \approx k_{\eta }$. Then we will match the two rates to find the scale $k_{\textrm {eq}}$ where the two rates balance. The first scale is characterized by a the time $t_{\gamma , D}$ for the cascade to reach $k_{\eta }$. and the second is characterized by the resistive time, $t_{\eta ,D}$, at $k_{\eta }$. Once the two times have been found a simple iteration leads us to $k_{\textrm {eq}}$. We will find that the two times do not balance. We will then alter $k$ to get balance, and find $k_{\textrm {eq}} = 4.5 k_{\eta }$.

The mathematical details of the first stage are given in the next section and the second stage in appendix C.

The total life of the instability from its initial time to the time when it is terminated is

(4.2)\begin{equation} t_D = t_{D \gamma } + t_{D \eta} \end{equation}

with the times first evaluated for $k = k_{\eta }$ and then for $k_{\textrm {eq}}$.

4.1. The rate of mass flow across magnetic lines

We can express the flow rate in terms of the instability lifetime $t_D$ as follows. After the instability is terminated, the initial density gradient has been reduced by

(4.3)\begin{equation} \delta \rho' = \rho_0 \varDelta / C , \end{equation}

where $\rho '$ is the initial gradient and the final gradient is zero, see (3.9).

After the cascade has terminated an amount of mass per unit area

(4.4)\begin{equation} \delta M = \frac{\xi^2 \delta \rho'}{ 2} = \frac{\varDelta}{C} \frac{\rho_0 \xi^2}{2 } \ \textrm{g}\ \textrm{cm}^2\ \textrm{s}^{-1} \end{equation}

will have been transferred across the plane $s =0$ by the cascade, see figure 6. Dividing this mass by $t_D$, the life of the instability, we obtain the average rate of the local mass flow across the magnetic lines during the life of this instability

(4.5)\begin{equation} \frac{ \delta M}{t_D} = \frac{\rho_0 \xi^2}{2 t_D} \frac{\varDelta}{C} . \end{equation}

So far, we have concentrated on the flow of local mass in region $\mathcal {R}$ crossing the magnetic field the field line. However, in a steady state all the mass in the entire flux tube must cross the field (see figure 7). However, the flow along the field line is easy so that as the local mass is removed from the unstable region, a pressure gradient will exist along the line which will pull the plasma to the unstable region. This happens fast enough that one sees a continual flow of the total mass across the field. If $\epsilon$ is the fraction of the line's length in the unstable region, then the amount of mass needed to be transferred will be larger by a factor approximately $1/\epsilon$. However, the time for complete transfer will be longer by the same factor. Thus, (4.5) still gives the mean flow during the life of the instability.

Figure 7. The correction for the flow along field lines. The mass along the entire fluid slows along the tube until it reaches the instability and crosses the line at that point, increasing the length of the flow time, but increasing the amount of mass during the longer time.

Assuming that the instabilities are continually excited with the same $\varDelta$, this expression yields the relation between the mean cross-flow rate and $\varDelta$, once the life of the instability $t_D$ is known. Our goal is thus reduced to evaluating $t_D$.

5. The evolution of the cascade

In this section, we derive the ideal equation for the time-dependent spectrum of the intensity of the magnetic fluctuations at a general point and time of the quasi-static equilibrium.

Let the field strength be

(5.1)\begin{equation} B=B_0 + b(\boldsymbol{x} ,t), \end{equation}

where $B_0$ is the unperturbed quasi-static magnetic field. The initial unstable gradient of the magnetic field strength is included in $b$.

We assume $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {v} =0$ because the magnetic field is so strong. Then the magnetic field strength satisfies the scalar equation

(5.2)\begin{equation} \frac{\partial b}{\partial t} = -\boldsymbol{\nabla} \boldsymbol{\cdot} (\boldsymbol{v} b) . \end{equation}

If we perform a Fourier analysis of $b$,

(5.3)\begin{equation} b= \varSigma_{\boldsymbol{k}} b_{\boldsymbol{k}}(t) \exp({\textrm{i} \boldsymbol{k} \boldsymbol{\cdot} \boldsymbol{x} }), \end{equation}

then this equation can be written as

(5.4)\begin{equation} \frac{\partial b_{\boldsymbol{k}} }{\partial t} = - \textrm{i} \boldsymbol{k} \boldsymbol{\cdot} \int \varSigma_{\boldsymbol{k''}} \left( \boldsymbol{v}_{\boldsymbol{k''}} b_{\boldsymbol{k'} } \right) \,\textrm{d} \omega , \end{equation}

where we use the notation

(5.5)\begin{equation} \boldsymbol{k} = \boldsymbol{k'} + \boldsymbol{k^{\prime\prime}} . \end{equation}

This equation can be interpreted as a $\boldsymbol {k}$ magnetic mode being changed by the interaction of a $\boldsymbol {k''}$ velocity mode with another $\boldsymbol {k'}$ magnetic mode. The integration over $\textrm {d}\omega$ allows for the time dependence of the turbulent velocity.

This equation is of the same type as the dynamo equation treated by Kulsrud & Anderson (1992), which itself is equivalent to the direct-interaction formalism of Kraichnan & Nagarajan (1967). It can be solved in a similar way with two changes. Here the perturbed magnetic field strength $b$ is the passive scalar. Further, the velocities are perpendicular to the unperturbed $\boldsymbol {B}_0$ and independent of the position along the unperturbed field line. The details of the solution, with the appropriate modifications, are given in appendix B.

We have used discrete Fourier modes for the analysis because they make the direct interaction between the turbulent velocity and the magnetic field clearer. The Fourier series are most easily interpreted in terms of continuous functions, so with box normalization, we replace an axisymmetric sum of the square of the Fourier modes by

(5.6)\begin{equation} \varSigma_k |b_{\boldsymbol{k}}|^2 = \int M(k)\,\textrm{d}k. \end{equation}

Thus, $M(k)$, is the spectrum of the intensity of the perturbations in the magnetic field strength. In the box normalization formalism $M(k)$ is defined by

(5.7)\begin{equation} M(k) = \left(\frac{L}{2 {\rm \pi} } \right)^2 \int |b_{\boldsymbol{k}}|^2 k\,\textrm{d} \hat{\boldsymbol{k}}, \end{equation}

where $L^2$ is the two-dimensional volume of our box, $k$ is a scalar, $\boldsymbol {\hat {k}}$ is a unit vector and the integral is over a cylindrical surface of radius $k$.

As shown in appendix B, for $k \gg k_0$ the time evolution of $M(k)$ is given by the differential equation

(5.8)\begin{equation} \frac{\partial M}{\partial t} = \frac{\gamma }{16} \left( k^2 \frac{\partial^2 M}{\partial k^2} +k \frac{\partial M}{\partial k} - M \right), \end{equation}

where from (B 28)

(5.9)\begin{equation} \gamma = \left( \frac{L}{2 {\rm \pi} } \right)^2 \int k''^2 J(k'',0) d \boldsymbol{k}'' . \end{equation}

(This ignores the fact that (5.8) for $\partial M(k)/\partial t$ is actually the large $k$ limit of an integral equation. Assuming the differential equation is valid down to $k_0$ does make a small error. However, because the range in $k$ is very large this is probably not serious.)

We have assumed that the velocity spectrum is concentrated at $k_0$. Then from this equation for $\gamma$ and (3.17)

(5.10)\begin{equation} v_0^2 = 2 \left(\frac{L}{2 {\rm \pi} } \right)^2 \int J((k^{\prime\prime},\omega)\,\textrm{d}^2 \boldsymbol{k}^{\prime\prime}\,\textrm{d} \omega , \end{equation}

it follows that

(5.11)\begin{equation} \frac{\gamma}{ v_0^2} = \frac{k_0^2}{2} \frac{ J(k,0)}{\int J(k,\omega)\,\textrm{d} \omega } = \frac{k_0^2 }{2 (\varDelta \omega) } \sim \frac{k_0^2}{2 k_0 v_0} , \end{equation}

where $v_0$ is the root-mean-square (rms) velocity. In this equation, $(\varDelta \,\omega)$, is the decorrelation rate of the turbulence, and is roughly equal to $k_0 v_0/2$, one half the turnover rate of the eddy. Thus, a reasonable approximation for $\gamma$ is the nonlinear growth rate $\gamma _{\textrm {nl}}$ given in (3.15). We see that $\gamma = \gamma _{\textrm {nl}}$ is also the turnover rate of the turbulent velocity, and is given in our numerical example in § 3 as $4.4 \times 10^{ 7 } \sqrt { \varDelta /C}\ \textrm {s}^{-1}$.

To solve (5.8), we transform the independent variable $k$ to

(5.12)\begin{equation} x= \ln (k/k_0) , \end{equation}

Then, the equation for $M$ becomes

(5.13)\begin{equation} \frac{\partial M}{\partial t} = \frac{\gamma }{16} \left( \frac{\partial^2 M}{\partial x^2} - M \right) . \end{equation}

This is a damped diffusion equation. Its solution, normalized so that the initial level of fluctuations $\left [\int M(k) \,\textrm {d}k =k_0 \int M(x) \textrm {e}^{x}\,\textrm {d}x\right]$ is unity, at $t=0$ is

(5.14)\begin{equation} M(x,t) =\frac{4 }{ k_0 \sqrt { {\rm \pi} \gamma t}} \exp({- \gamma t/16 -4 x^2/\gamma t}) . \end{equation}

(In fact, $M$ should be normalized by the initial condition $(\delta B)^2/B_0^2= (\gamma \beta /2 C) \varDelta$, but we do not make any use of this normalization.)

Although it turns out that $k_{\textrm {eq}}$ is not exactly equal to $k_{\eta }$ it is of the same order of magnitude. Thus, we are first interested in the integrated energy about $k_{\eta }$,

(5.15)\begin{equation} k_{ \eta } M(k_{ \eta },t) = \frac{4 }{\sqrt{ {\rm \pi} \tau }} \exp({- \tau /16 -4 x_{\eta}^2/\tau + x_{\eta} }), \end{equation}

where $k_{\eta } = \sqrt { \gamma /\eta }$, $\tau = \gamma t$ and $x_{\eta }= \ln (k_{\eta } /k_0$. If we take $\eta =1.0\ \textrm {cm}^2\ \textrm {s}^{-1}$ and $\varDelta = 10^{-4}\ \textrm {cm}^{-1}$, we find $x_{\eta } =9.26$.

The value of the resistivity in the conditions of the neutron star shell is somewhat controversial. (As the resistivity enters into our results only logarithmically, we may choose the value of unity given by Shapiro & Teukolsky (1983) for it.)

For a fixed $x_{\eta }$, (5.15) is a function of $\tau$ alone. Its asymptotic expression for large $x_{\eta }$ can be found by expanding its exponent about its maximum value $\tau = \tau _0 =8 x_{\eta }$. The result is

(5.16)\begin{equation} k_{ \eta } M(k_{ \eta } ,t)= \frac{4 k_0 } {\sqrt{ {\rm \pi} \tau_0 }} \exp{ \frac{-(\tau-\tau_0)^2}{16 \tau_0}} . \end{equation}

Thus, the time for the eddies to cascade down to $k_{\eta }$ is $t_{D \gamma } = \tau _0/\gamma = 8 x_{\eta } /\gamma =74 /\gamma \ \textrm {s}$.

A similar result holds for the density whose spectrum is given by $N(k,t)$.

These results give the initial condition for the second stage: the resistive evolution of $b$ and $\delta$, the perturbed magnetic strength and density, respectively, at $k_{\eta }$. In appendix C, it is shown that $\delta$ and $b$ satisfy a third-order system of differential equations. Their solution can be expressed in terms of three normal modes. One mode damps at the rate $k_{\eta }^2 \eta /2 = \gamma /2$, and the other two at the rate $\gamma /4$. Therefore, all these modes should disappear in a time of order $4/\gamma$.

Hence, it would seem that the total time during which the instability lasts, $t_D$, is

(5.17)\begin{equation} t_D = t_{D \gamma } + t_{D \eta } = \frac{74 }{ \gamma} + \frac{4}{\gamma }, \end{equation}

we see that the first term dominates.

However, the two times, or rates, are different for stage one and two. This is because we chose $k_{\eta }$ only as a first guess for $k_{\textrm {eq}}$. Thus, our guess for $k_{\textrm {eq}}$ must be adjusted. Note that the first stage time $t_{\gamma ,D}$ depends only logarithmically on the terminal $k$, but the second stage time depends on its square. Thus, as a first iteration we may keep $t_{\gamma ,D}$ fixed and vary $k$ to make the second time equal to the first. This gives $k_{\textrm {eq}} = 4.3 k_{\eta }$ and the total time will be just twice the unchanged first time.

Thus, the total decay time is

(5.18)\begin{equation} t_D = \frac{16 x_{\eta} }{\gamma }= \frac{148}{k_0 \xi \sqrt{ g_s \varDelta/C}} = \frac{24}{\sqrt{g_s \varDelta/C} } . \end{equation}

In this formula, we use $x_{\eta } = 9.26$ rather than the more correct value of $x_{\eta } =13.87-0.5 \times \ln (\varDelta /C)$. This is slightly inaccurate, but it keeps the results more compact.

5.1. The mean rate of flow across the field lines

Let $u_F$ be the velocity of mass transfer across the field lines. Then $\rho _0 u_F$ is given by the total mass divided by $t_D$, which from (4.5) is

(5.19)\begin{equation} \rho_0 u_F = \frac{\delta M}{ t_D}= \frac{1 }{48 } \rho_0 \sqrt{ g_s} \xi^2 \left( \frac{\varDelta}{C} \right)^{3/2} , \end{equation}

thus

(5.20)\begin{equation} u_F= \frac{1}{48} \sqrt{ g_s} \xi^2 \left(\frac{\varDelta}{C} \right)^{3/2} \ \textrm{cm}\ \textrm{s}^{-1}. \end{equation}

To determine the mixing length $\xi$ we use the fact that the critical quantity $\varDelta$ is only positive over a limited range $-\xi < s < +\xi$ as in the figure 3. Now, when $\varDelta$ is considered to be a function of $s$, then $2 \xi$ is the difference between its zeros. The maximum value of $\varDelta (s)$ is the critical value $\varDelta$ so $\varDelta = -(\xi ^2/2) (\textrm {d}^2 \varDelta)(\textrm {d} s^2)$. However, $\varDelta$ is itself a logarithmic derivative, thus

(5.21)\begin{equation} (\textrm{d}^2 \varDelta/\textrm{d} s^2) =- \frac{\textrm{d}^3 \ln (p/B^{\gamma})}{\textrm{d} s^3} , \end{equation}

is of order $-1/\lambda ^3$, where the length, $\lambda$, depends on the actual equilibrium. Without knowing the detailed equilibrium obtainable from a numerical calculation, we can only guess at the value of $\lambda$, but presumably it is of the order of magnitude of the two scale heights of $B$ and $p$. Estimating $\xi$ in this way, we have

(5.22)\begin{equation} \xi^2 \approx 2 \lambda^3 (\varDelta/C) . \end{equation}

Thus, with this estimate for $\xi$,

(5.23)\begin{equation} u_F =\frac{1}{48} \sqrt{ g_s} \left( \frac{\varDelta}{C}\right)^{3/2} \xi^2 =\frac{1 }{24 } \lambda^3 \sqrt{ g_s} \left( \frac{\varDelta }{C} \right)^{5/2}\ \textrm{cm}\ \textrm{s}^{-1}. \end{equation}

This result gives the relation between the instability parameter $\varDelta$ and the flow velocity of matter across the strong neutron star field lines. It can be considered indicative of the rate of loss of flux freezing.

5.2. The magnitude of $\varDelta$

We have stated that positive values of $\varDelta$ in our normal solutions can be treated small enough as to be considered negligible. To indicate that this is true, we estimate mass flow across a general surface. Then use (5.23) to calculate $\varDelta$ and verify that it is indeed small.

Let $\dot {M_{\textrm {total}}}$ be the time derivative of $M_{\textrm {total}}$ equal to $6 \times 10^{ 17} L$ grams per second, which is equivalent to $10^{ -9 } L$ solar masses per year. In a quasi-steady state $\dot {M_{\textrm {total}}}$ is the same for all surfaces except those close to the poloidal axis. The surface area of any given flux tube is its mean radius $r$ times its vertical curved length, which we approximate as $D_c$, the depth of the crust. Then the mean flow velocity across the surface, $U_F$, is given by

(5.24)\begin{equation} \dot{M_{\textrm{total}}} = 2 {\rm \pi} \rho_0 r D_c U_F . \end{equation}

For $D_c = 1\ \textrm {km}$ and $\rho _0 = 10^{ 11}\ \textrm {g}\ \textrm {cm}^{-3}$, we find

(5.25)\begin{equation} U_F = \frac{\dot{M_{\textrm{total}}}}{2 {\rm \pi} \rho_0 r D_c } \approx \frac{10 L}{r}\ \textrm{cm}\ \textrm{s}^{-1} . \end{equation}

This velocity is slow because the huge mass near the crust flows along the field to the instability and across the flux line.

Taking $\lambda \approx 10^{4}\ \textrm {cm}$ and $g_s =2 \times 10^{ 14 }\ \textrm {cm} \ \textrm {s}^{-2}$, we find from (5.23) that

(5.26)\begin{equation} u_F = 5 \times 10^{17} \left( \frac{\varDelta }{C } \right)^{5/2}\ \textrm{cm}\ \textrm{s}^{-1} , \end{equation}

so equating this to $U_F$ ,

(5.27)\begin{equation} \frac{\varDelta }{C} = (5 \times 10^{ -17} L/r)^{2/5} \approx 2.1 \times 10^{ -7} \left( \frac{L}{r} \right)^{2/5}\ \textrm{cm} . \end{equation}

(Here $r$ is in centimetres. For $r$ in metres, we have $\varDelta = 3.3 \times 10^{ -8 } (L/r_m)^{2/5}$.) Remembering that an arbitrary $\varDelta$ would be of order $10^{ -4}$, its small size here is due to the cancellation of $p'/p$ and $\gamma B'/B$. We see that the effective value of $\varDelta$ is smaller than expected by a factor of $2.1 \times 10^{ -3 }/r$.

Thus, we may consider the effective value $\varDelta$ to be very small. This estimate justifies our claim that $\varDelta$ is negligibly small.

6. The evolution of the quasi-steady equilibrium as a sequence of self-organized criticality states

We can make use of the smallness of $\varDelta$ to determine the evolution of the quasi-static states. We have treated the instabilities as starting from a marginal state and considering the amplitude growing to a magnitude that can balance the damping of the instability owing to the flow of mass across the neutron star surface, and then decreasing to zero. It would be more realistic to imagine the instabilities are always at the amplitude corresponding to equation (5.27).

We do not expect the amplitudes of each instability to be constant, but rather the rms of their amplitudes to be such as to lead to a balance between the dissipation rates of all the amplitudes and their growth rates.

That is, the mean value of the square of the amplitudes gives the correct flow according to (5.23), and the size of the $\varDelta$ is such as to continually balance the dissipation rate.

We do not need to know the actual value of $\varDelta$ because their magnitude is small enough that they do not affect the quasi-static equilibrium so we may simply take the size of the $\varDelta$ at the instability sites large enough to balance the damping of the unstable modes. We call a quasi-static state normal if there is a point on every line such that $\varDelta$ is zero and is non-positive everywhere else.

This situation is duplicated in stellar convection zones where, in practice, the radial entropy gradient $\boldsymbol {\nabla } \ln (p/\rho ^\gamma)$ can be taken as zero everywhere in the convection zone (except near the surface). There the gradient is small and just positive enough for the instability amplitude to drive the necessary amount of convective heat flow, but not large enough to be noticeable, so in the convective zone the gradient of the entropy can be taken to be zero.

The significant difference between our treatment of the distribution of mass on the neutron star surface and the distribution of temperature in the convection zone is that our $\varDelta$ need to be zero only on a limited region of each magnetic surface and negative everywhere else.

It should be noted that our critical quantity for stability $\varDelta = \boldsymbol {\nabla } [\log (p/B^\gamma)]$ is the gradient of the magnetic form of entropy.

Thus, to find the quasi-steady state at $t$ we find the normal state in which the total excess mass $M_{\textrm {total}}$ is equal to the total mass that has fallen up to this time. We believe such normal states only depend on $M_{\textrm {total}}$ and are unique solutions of the Grad–Shafranov equation. As $M_{\textrm {total}}$ increases from zero they give the correct evolution of the problem of fallen mass on the surface of a neutron star.

These normal states are marginally stable states that are often referred to in the literature as ‘self-organized criticality states’.

From the fact that in our steady state, the flow across the field would seem to occur so freely, one might conclude that the field itself would not be changing. This cannot be the case because the flow across the lines can only occur if there is an instability and this instability will only occur under certain conditions. These conditions change as the amount of accumulated mass changes and to accommodate this the magnetic field itself must change. To find these changes, one must rely on there always being instabilities and the system must in a marginal state and satisfy the Grad–Shafranov equation.

In appendix D, we speculate as to the long-time behaviour of the evolution of the neutron star's external shell.

7. Comparison with other research

There has been a large number of papers presented on this problem. We found the following to be noteworthy.

Hameury et al. (1983) were the first to derive the Grad–Shafranov as the equation for equilibrium. They solved this equation for special cases of what we refer to as $h(\psi)$ equilibrium equations. Then Geppert & Urpin (1994), and later Cumming, Zweibel & Bildsten (2001), assumed that the incoming mass was not magnetized, but they included resistivity to see whether the mass penetrated the neutron star field. Brown & Bildsten (1998) assumed that the mass did not spread out radially, but simply sunk into the neutron star. Litwin, Brown & Rosner (2001) discovered and analysed the Newcomb instability of the equilibrium and determined critical conditions for instability. Melatos & Phinney (2001), in a seminal paper, investigated conditions when a very large amount of infalling neutral mass, comparable with $10^{ -5}$–$10^{ -2 }$ solar masses, accumulated on the neutron star. They could not be definite about their results because of their uncertainty as to whether flux freezing or resistive diffusion dominated. However, they did discuss the possibility that this amount of mass could reduce the magnetic moment of the neutron star. Then Payne & Melatos (2004, 2007), Vigelius & Melatos (2008) and others made extensive numerical simulations of the time evolution of the mass. However, in these papers it was assumed that the neutron star surface was rigid so that none of the new mass penetrated the surface. They further explored the depression of the magnetic moment and investigated the size of the resulting quadrupole moment of the neutron star, to see whether there could be significant emission of gravitational radiation. However, because it is now appreciated that the mass does penetrate the surface, the significance of their results is uncertain. In most of these papers, the emphasis was on nuclear reactions that would be produced by the infalling mass rather than on the equilibrium itself.

Priymak, Melatos & Payne (2011) returned to the stability question. They approached it by time-integrating the dynamic equations after a strong perturbation was produced by adding random extra mass. This led to some results on the gross instabilities, but they could not reach the small-scale instabilities discussed by Litwin et al. (2001) and those studied in our present paper.

Mukherjee & Bhattacharya (2012) and Mukherjee, Bhattacharya & Mignone (2013) carried out a similar approach to the stability question with similar results.

Some attempts (e.g. Urpin & Geppert 1995; Urpin 2005) were made to find mechanisms to enhance the resistivity by the effects of the infalling mass, but apparently they were not successful.

8. Summary and conclusions

The main goal of this paper has been to show that when external mass flows onto a neutron star, it may not be inhibited by flux freezing, from flowing across the neutron star's magnetic field from the pole towards the equator. We have shown that when the flow of mass from the other star sets in, stresses lead to an instability that allows mass to cross field lines. Although the instability is ideal on its initial scale, it nonlinearly generates a cascade of magnetic field strength fluctuations that extends down to resistivity scales where flux freezing is broken and matter is allowed to cross the field lines.

By an analytical analysis, we show that the instability is so strong that even a small increase in the parameter for instability, $\varDelta$, above its critical condition, will lead to the necessary amount of flow to balance the inflow.

This flow leads to a redistribution of the excess mass per flux over the neutron star surface. (This distribution of mass uniquely determines the magnetostatic equilibrium.) However, the latter is very difficult to derive from the properties of the instability, because the amount of flow is so sensitive to its amplitude.

Instead the flow is determined indirectly at any time $t$, by finding a distribution of the mass that allows a normal solution of the equations for a quasi-static equilibrium. We call a solution normal if it leads to an instability on every flux surface. (An additional requirement is that the equilibrium contain the amount of excess mass that has accumulated up to this time.) The necessary flow between the states can be found from the different equilibria at neighbouring times. These flows can be used to calculate the amplitude of the instability that we have shown to be very small.

In this paper, we first derive the equation for magnetostatic equilibrium and show that its solution depends on the distribution of mass per flux. We next derive the ideal instability and show that it indeed leads to a cascade of fluctuations.

Then we carry out a two-dimensional ‘Kulsrud–Anderson’ analysis (Kulsrud & Anderson 1992) (first developed for dynamo theory) to derive the equation for the time-dependent evolution of the cascade. We apply resistive MHD theory to calculate the amount of cross-field flow that results from these fluctuations and relate it to the marginal value $\varDelta$. Then we demonstrate that, for the rate of flow of mass onto the neutron star, very small values of $\varDelta$ are needed.

Because these values of $\varDelta$ are so small, this is all the information we need from the instability. The whole problem reduces to looking for the sequence of normal states. This result was our goal.

Based on these results and making use of a very crude marginal condition for instability, we speculate in appendix D about what might happen on longer time scales. This approximate condition can be expressed in terms of $h(\psi)$, the variation of the height of the mass above the ambient neutron star surface. It is $\textrm {d} h/\textrm {d}r < \cos \alpha$, where $r$ is the distance measured along the neutron star surface and $\alpha$ is a certain angle, which we cannot ascertain without a detailed numerical calculation.

Assuming that mass flow on the star is $10^{ -9 }$ solar masses per year, we find that it will take $\approx 12/a$ thousand years for the mass to finally reach the equator, where $a = \cos \alpha$. The height of the mass will be ($1.5/a$) kilometres at the poles and will be a linear decreasing function of $r$.

We found that our normal states may not stretch all the way to the equator with free cross-field flows. Instead, if the neutron star field is less than $10^{13}$ Gauss, there will be an unmagnetized cylindrical region near the neutron star pole where the plasma may be field free and allow the flowing mass to fill this unmagnetized region with no trouble. This region would be surrounded by a normal type of solution. This, in turn, could be surrounded by an undisturbed region where the remaining flux is compressed.

Needless to say these speculations are unlikely to be even qualitatively valid, because they are based on a highly approximate treatment of the instability. However, they do yield interesting results. They predict the form of the altered surface magnetic field over the entire surface and indicate a time-dependent macroscopic change in the neutron star dipole field. After a few tens of thousand of years, a steady state for the magnetic field will be reached. The neutron star will still increase its mass owing to more mass falling on it and its radius will increase at about one-tenth of a centimetre per year, but the distribution of mass and magnetic field will be in a steady state and the neutron star will slowly increase its radius under all this new mass to accommodate its inflow.