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 (Reference Shapiro and Teukolsky1983). 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 (Reference Schwarzschild1958), 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 (Reference Newcomb1961), 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 (Reference Kulsrud and Anderson1992), which itself is equivalent to the direct-interaction formalism of Kraichnan & Nagarajan (Reference Kraichnan and Nagarajan1967). 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 (Reference Shapiro and Teukolsky1983) 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. (Reference Hameury, Bonazzola, Heyvaerts and Lasota1983) 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 (Reference Geppert and Urpin1994), and later Cumming, Zweibel & Bildsten (Reference Cumming, Zweibel and Bildsten2001), assumed that the incoming mass was not magnetized, but they included resistivity to see whether the mass penetrated the neutron star field. Brown & Bildsten (Reference Brown and Bildsten1998) assumed that the mass did not spread out radially, but simply sunk into the neutron star. Litwin, Brown & Rosner (Reference Litwin, Brown and Rosner2001) discovered and analysed the Newcomb instability of the equilibrium and determined critical conditions for instability. Melatos & Phinney (Reference Melatos and Phinney2001), 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 (Reference Payne and Melatos2004, Reference Payne and Melatos2007), Vigelius & Melatos (Reference Vigelius and Melatos2008) 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 (Reference Priymak, Melatos and Payne2011) 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. (Reference Litwin, Brown and Rosner2001) and those studied in our present paper.

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

Some attempts (e.g. Urpin & Geppert Reference Urpin and Geppert1995; Urpin Reference Urpin2005) 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 Reference Kulsrud and Anderson1992) (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.