Pinch dynamics in a low-$\beta$ plasma

TL;DR

This paper investigates the relaxation dynamics of a low-$eta$ plasma with helical magnetic fields, exploring the effects of an $oldsymbol{ extalpha}$-effect and demonstrating that the field does not relax to a Taylor state under the studied conditions.

Contribution

It introduces an exploratory model of $oldsymbol{ extalpha}$-effect in plasma relaxation, showing how it influences magnetic field evolution and prevents relaxation to a Taylor state.

Findings

Weak $oldsymbol{ extalpha}$-effect causes oscillations and reversed fields.

Field does not become a Taylor state due to non-uniform $oldsymbol{ extgamma}$.

High-frequency oscillations occur when $oldsymbol{ extalpha}$ exceeds a critical value.

Abstract

The relaxation of a helical magnetic field ${\bf B}({\bf x}, t)$ in a high-conductivity plasma contained in the annulus between two perfectly conducting coaxial cylinders is considered. The plasma is of low density and its pressure is negligible compared with the magnetic pressure; the flow of the plasma is driven by the Lorentz force and and energy is dissipated primarily by the viscosity of the medium. The axial and toroidal fluxes of magnetic field are conserved in the perfect-conductivity limit, as is the mass per unit axial length. The magnetic field relaxes during a rapid initial stage to a force-free state, and then decays slowly, due to the effect of weak resistivity $\eta$, while constrained to remain approximately force-free. Interest centres on whether the relaxed field may attain a Taylor state; but under the assumed conditions with axial and toroidal flux conserved inside every cylindrical Lagrangian surface, this is not possible. The effect of an additional $\alpha$-effect associated with instabilities and turbulence in the plasma is therefore investigated in exploratory manner. An assumed pseudo-scalar form of $\alpha$ proportional to $q\,\eta\, ({\bf j}\cdot {\bf B})$ is adopted, where $ {\bf j}=\nabla\times {\bf B}$ and $q$ is an $\mathcal{O}(1)$ dimensionless parameter. It is shown that, when $q$ is less that a critical value $q_c$, the evolution remains smooth and similar to that for $q=0$; but that if $q>q_c$, negative-diffusivity effects act on the axial component of $\bf B$, generating high-frequency rapidly damped oscillations and an associated transitory appearance of reversed axial field. However, the scalar quantity $\gamma={\bf j}\cdot {\bf B}/B^2$ remains highly non-uniform, so that again the field shows no sign of relaxing to a Taylor state for which $\gamma$ would have to be constant.

Equations

$${\mathcal{H}}_V={\int }_V\mathbf{A}\cdot \mathbf{B}\text{d}V,$$

$$M(t)=\frac{1}{2}{\int }_{\mathcal{D}}{\mathbf{B}}^2\text{d}V$$

$$M(t)>\lambda |\mathcal{H}|,$$

$$\mathbf{B}=B_0\left(0,b_{\theta }(r,t),b_z(r,t)\right).$$

$$\mathbf{j}=\nabla \times \mathbf{B}=B_0\left(0,-\frac{\partial b_z}{\partial r},\frac{1}{r}\frac{\partial }{\partial r}(r{b}_{\theta })\right).$$

$$\mathbf{u}=\left(u(r,t),0,0\right)$$

$$\partial \mathbf{B}/\partial t=\nabla \times (\mathbf{u}\times \mathbf{B}),$$

$$0=\mathbf{j}\times \mathbf{B}+{\mu }_s{\nabla }^2\mathbf{u}+\left(\frac{1}{3}{\mu }_s+{\mu }_b\right)\nabla (\nabla \cdot \mathbf{u}),$$

$$dM/dt=-{\int }_V\left[{\mu }_s(\nabla \times \mathbf{u}{)}^2+\left(\frac{4}{3}{\mu }_s+{\mu }_b\right)(\nabla \cdot \mathbf{u}{)}^2\right]\text{d}V.$$

$$\mathbf{j}=\gamma \mathbf{B},$$

$$\eta \mathbf{n}\times \mathbf{j}=\mathbf{n}\times \mathbf{E}=0\text{on}r=\delta \text{and on}r=1,$$

$$\eta \partial b_z/\partial r=0,\eta \partial (r{b}_{\theta }/)\partial r=0,\text{on}r=\delta \text{and on}r=1.$$

$$\frac{\partial \mathbf{B}}{\partial t}=-\nabla \times \mathbf{E}=\nabla \times (\mathbf{u}\times \mathbf{B})+\eta {\nabla }^2\mathbf{B}.$$

$$\frac{\partial \rho }{\partial t}=-\nabla \cdot (\rho \mathbf{u}).$$

$$\frac{D}{Dt}\left(\frac{\mathbf{B}}{\rho }\right)=\left(\frac{\mathbf{B}}{\rho }\cdot \nabla \right)\mathbf{u}+\frac{\eta }{\rho }{\nabla }^2\mathbf{B},$$

$$(\mathbf{B}\cdot \nabla )\mathbf{u}=B_0(0,b_{\theta }u/r,0),$$

$$\frac{D}{Dt}\left(\frac{b_z}{\rho }\right)=0\text{and}\frac{D}{Dt}\left(\frac{b_{\theta }}{r\rho }\right)=0,$$

$$\frac{\partial b_{\theta }}{\partial t}=-\frac{\partial }{\partial r}(u{b}_{\theta })+\eta \frac{\partial }{\partial r}\frac{1}{r}\frac{\partial }{\partial r}(r{b}_{\theta }),$$

$$\frac{\partial b_z}{\partial t}=-\frac{1}{r}\frac{\partial }{\partial r}(ru{b}_z)+\eta \frac{1}{r}\frac{\partial }{\partial r}r\frac{\partial b_z}{\partial r}.$$

$${\Phi }_z={\int }_{\delta }^1{b}_z(r,t)2\pi r\mathrm{d}r,$$

$$\frac{1}{2\pi }\frac{\mathrm{d}{\Phi }_z}{\mathrm{d}t}={\int }_{\delta }^1\frac{\partial b_z}{\partial t}r\mathrm{d}r={\int }_{\delta }^1\frac{\partial }{\partial r}\left(-ru{b}_z+\eta r\frac{\partial b_z}{\partial r}\right)\mathrm{d}r={\left[-ru{b}_z+\eta r\frac{\partial b_z}{\partial r}\right]}_{\delta }^1=0.$$

$${\Phi }_{\theta }={\int }_{\delta }^1{b}_{\theta }(r,t)\mathrm{d}r,$$

$${\Phi }_z=\text{cst.},{\Phi }_{\theta }=\text{cst.},$$

$$\mathbb{M}={\int }_{\delta }^1\rho (r,t)2\pi r\mathrm{d}r=\text{cst.}$$

$$b_z(r,0)=\frac{1}{\pi (1-{\delta }^2)}(\text{so}{\Phi }_z=1),$$

$$r{b}_{\theta }(r,0)=\frac{c}{3(1-\delta {)}^3}\left[3(r-\delta {)}^2(1-\delta )-2(r-\delta {)}^3\right]{\text{e}}^{-k(1-r{)}^2}.$$

$$\delta =0.5,c=0.5,k=10,$$

$$u(r,0)=0,\rho (r,0)={\rho }_0,$$

$$\mathbf{j}\times \mathbf{B}=B_0^2\left(-\frac{1}{2}\frac{\partial }{\partial r}(b_{\theta }^2+b_z^2)-\frac{b_{\theta }^2}{r},0,0\right).$$

$$\rho \left(\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial r}\right)=-B_0^2\left(\frac{1}{2}\frac{\partial }{\partial r}(b_{\theta }^2+b_z^2)+\frac{b_{\theta }^2}{r}\right)+\mu \left(\frac{{\partial }^{2}u}{\partial r^2}+\frac{1}{r}\frac{\partial u}{\partial r}-\frac{u}{r^2}\right),$$

Full text

Keywords: Magnetic relaxation; Taylor conjecture; $\alpha$-effect; reversed-field pinch

1 Introduction

Magnetic relaxation is the process by which a magnetic field in a highly conducting fluid seeks a minimum energy state subject to pertinent topological constraints (Moffatt 1985). In the perfect conductivity limit, the lines of force (‘$\bf B$-lines’) are frozen in the fluid, a topological constraint represented by the family of magnetic helicity invariants

[TABLE]

where ${\bf B}=\nabla\times{\bf A}$, and $V$ is any Lagrangian volume on whose surface $\partial V$ (with unit normal $\bf n$)  ${\bf n}\cdot{\bf B}=0$. These invariants represent the conserved degree of linkage of $\bf B$-lines within $V$. Of particular importance is the global magnetic helicity $\cal H$, integrated over the whole domain $\cal D$ of fluid. The magnetic energy

[TABLE]

then has a lower bound (Arnold 1974)

[TABLE]

where $\lambda>0$ is a constant that depends on the scale and geometry of $\cal D$.

The particular problem addressed in this paper concerns relaxation in a plasma of extremely low-density $\rho$, in which the fluid pressure $p$ is negligible compared with the magnetic pressure $p_{M}={\textstyle\frac{1}{2}}\mu_{0}{\bf B}^{2}$, i.e. $\beta=p/p_{M}\ll 1$. Flow of the plasma is then driven solely by the Lorentz force ${\bf j}\times{\bf B}$. The low density implies further that inertia is negligible compared with the viscous force in the Navier-Stokes equation, the viscosity $\mu$ being essentially independent of $\rho$ in the limit $\rho\rightarrow 0$. These approximations have been adopted in a cartesian model by Bajer & Moffatt (2013) who treated relaxation of a single-component field, and by Moffatt (2015) who considered the case of a two-component field with non-zero helicity. Here we shall consider the situation in a cylindrical geometry, for which the Lorentz force includes the ‘hoop stress’ associated with curvature of the $\bf B$-lines.

2 Relaxation in a low-$\beta$ plasma

We consider a two-component helical field in cylindrical polar coordinates $\{r,\theta,z\}$ of the form

[TABLE]

The associated current distribution is given by333For simplicity of notation, we absorb the conventional constant $\mu_{0}$ in the definition of $\bf j$.

[TABLE]

The first objective is to determine how the relaxation of this field is constrained by the initial magnetic helicity distribution. The component $b_{\theta}$ of the field is responsible for the classic ‘pinch effect’ (Bennet 1934). In this scenario, it is natural to suppose that the initial $z$-component of field is uniform, and that the initial $\theta$-component is concentrated near the outer cylindrical boundary, such a field then providing a radial Lorentz force that tends to drive the plasma inwards. We shall suppose that the resulting radial motion

[TABLE]

is controlled by viscosity, which, as indicated above, dominates over inertia when the plasma density is sufficiently small. We further suppose that the fluid pressure is negligible compared with the magnetic pressure, i.e. this is a ‘low-$\beta$’ plasma. Of course, the fluid pressure increases in the inner region where the density increases; the effect of this increase can be included without difficulty in the numerical treatment.

In the perfect conductivity limit $\eta=0$, the magnetic field evolves according to the ‘frozen-field’ equation

[TABLE]

and, with the neglect of inertia and pressure gradient, the Navier-Stokes equation degenerates to

[TABLE]

where $\mu_{s}$ and $\mu_{b}$ are the shear and bulk viscosities. From (7) and (8), an equation may easily be derived for the magnetic energy:

[TABLE]

Thus $M(t)$ is monotonic decreasing and bounded below by (3). Equilibrium is attained only when ${\bf u}\equiv 0$, and then from (8),  ${\bf j}\times{\bf B}=0$, i.e. the field is ‘force-free’. Hence

[TABLE]

for some pseudo-scalar function $\gamma(\bf x)$ satisfying $({\bf B}\cdot\nabla)\,\gamma=0$.

For the particular one-dimensional geometry considered in this paper, it is therefore to be expected that, when the magnetic diffusivity $\eta$ is sufficiently weak, the field will relax rapidly during an initial stage to a force-free state (with here $\gamma=\gamma(r)$) that has minimum energy compatible with its initial (conserved) topology. Minimising energy subject to the single topological constraint of conserved global helicity yields a force-free field structure with $\gamma=$cst., a condition that provides reversed axial field near the outer boundary (Taylor 1974). However, the dynamical process through which such a reversed field may spontaneously appear is not revealed by the simple process of seeking a minimum-energy state. Such a reversal cannot in fact appear for so long as the pinching motion is purely radial. However, it seems possible that instabilities of the basic relaxing field may lead to an $\alpha$-effect, which could conceivably achieve reversal. We shall explore this possibility in §8; first however, we treat simple radial relaxation neglecting any instabilities that may be present.

We suppose that the plasma is contained in the cylindrical annulus $\delta\,<r/a<1$, where $0<\delta<1$; we shall scale all lengths so that, in effect, $a=1$. The boundaries $r=\delta$ and $r=1$ are assumed to be thin perfectly conducting cylinders, separating the plasma from the internal and external regions; these boundaries can therefore support current sheets with both $z$- and $\theta$-components. The electric field is given by ${\bf E}=\eta\,{\bf j}-{\bf u}\times{\bf B}=(0,E_{\theta},E_{z})$, and we suppose that ${\bf E}=0$ in the internal and external regions, assumed insulating, i.e. for $r<\delta$ and $r>1$. Since $u=0$ on both boundaries and both tangential components of ${\bf E}$ are continuous across them, it follows that, with ${\bf n}=(1,0,0)$,

[TABLE]

i.e. that

[TABLE]

Some type of boundary-layer behaviour is to be expected in the limit $\eta\rightarrow 0$.

3 Field evolution and flux conservation

When the magnetic diffusivity $\eta$ is nonzero, the field evolution is described by the induction equation,

[TABLE]

and the fluid density $\rho$ satisfies the mass conservation equation

[TABLE]

These equations may be combined to give

[TABLE]

where $D/Dt\equiv(\partial/\partial t+{\bf u}\cdot\nabla)$, the Lagrangian (or ‘material’) derivative. Noting that, from (4) and (6),

[TABLE]

when $\eta=0,\,$eqn.(15) gives

[TABLE]

i.e. when following any material element of fluid, $b_{z}/\rho$ and $(b_{\theta}/r)/\rho$ are constant.

In Eulerian form, eqn.(13) has components

[TABLE]

and

[TABLE]

The flux of $b_{z}$ between the two cylinders is

[TABLE]

and we note, using (19) and the conditions $u=0,\,\eta\,\partial b_{z}/\partial r=0$ on both boundaries, that

[TABLE]

It follows that $\Phi_{z}=$ cst.

Similarly, the flux of $b_{\theta}$ in the $\theta$-direction, per unit axial length between the cylinders, is

[TABLE]

and it follows in the same way from (18), and the conditions $u=0,\,\eta\,\partial(rb_{\theta})/\partial r=0$ on both boundaries, that $\Phi_{\theta}=$ cst. also. The constraints

[TABLE]

provide an important check on the numerical computations that follow (see Fig. 2(f)). The mass per unit axial length $\mathbb{M}$ is of course similarly constant:

[TABLE]

The results (23) are clearly compatible with (17) when $\eta=0$.

4 Initial conditions

We adopt as initial conditions for the magnetic field

[TABLE]

and

[TABLE]

These are chosen to be compatible with the boundary conditions (12) for any values of the parameters $\{\delta,c,k\}$, and to satisfy $b_{\theta}(1,0)=c$. Fig. 1 shows the field $rb_{\theta}(r,0)$ for $\delta=c=0.5$ and for three values $k=1,\,10,\,30$, together with corresponding values of the flux $\Phi_{\theta}$. Increasing $k$ leads to increasing concentration of $rb_{\theta}(r,0)$ near the outer boundary. Rather arbitrarily we choose

[TABLE]

in the computations that follow.

We further adopt initial conditions for the velocity and density fields,

[TABLE]

i.e. the plasma is initially at rest with uniform density $\rho_{0}$.

5 Dynamics of the relaxation process

The Lorentz force in the cylindrical geometry considered here takes the form

[TABLE]

The Navier-Stokes equation, including this term, has only a radial component:

[TABLE]

where $\mu=4\mu_{s}/3+\mu_{b}$ is an effective viscosity. As in Moffatt (2015), it is now convenient to introduce dimensionless variables

[TABLE]

and dimensionless parameters

[TABLE]

we assume that both these parameters are small: $\kappa\!\ll\!1$ (i.e. small diffusivity); and $\epsilon\ll 1$ (i.e. low density). With the variables (31), and immediately dropping the hats, eqn. (14) is unchanged, while in eqns. (18) and (19), $\eta$ is simply replaced by $\kappa$:

[TABLE]

and

[TABLE]

The momentum equation (30) becomes

[TABLE]

Here, the term ${\textstyle\frac{1}{2}}\partial(b_{\theta}^{2}+b_{z}^{2})/\partial r$ is the gradient of magnetic pressure, and the term $-b_{\theta}^{2}/r$ is the additional ‘hoop stress’ that arises due to curvature of the ${\bf B}$-lines.

6 Numerical integration

We can now proceed to numerical integration of these equations, with boundary conditions as already stated. When $\kappa\neq 0$, these are

[TABLE]

and initial conditions as specified in the previous section. Results obtained with Mathematica are summarised in Fig. 2(a -h), for the particular choice of parameters (27), together with $\epsilon=0.01,\kappa=0.0001$ (the behaviour for this choice is quite typical). The panels of the figure show (a) the inward movement of the density field $\rho(r,t)$ in response to the negative radial Lorentz force (b), which collapses rapidly to near zero; (c,d) the corresponding evolution of the magnetic field components; (e) the rise of both field components at the inner boundary; (f) the decay of magnetic energy; (g) the pseudo-scalar coefficient $\gamma(r)$ when the nearly force-free state has been established; and finally (h) the fluxes $\Phi_{z}$ and $\Phi_{\theta}$ which, as expected, remain constant to within numerical error throughout the whole computational period $0<t<500$.

The following points are particularly worth noting. First, there are clearly two phases to the evolution: an initial phase (here $0<t\lesssim 7$) during which magnetic diffusion is negligible and the magnetic energy decreases relatively rapidly on the (dimensional) time-scale $\mu/B_{0}^{2}$; and a slow diffusive stage $t\gtrsim 7$, during which $b_{z}$ slowly relaxes back to its initial uniform value (its flux remaining constant); during this phase, $rb_{\theta}$ also slowly decays to a constant $C\,(=\Phi_{\theta}/\log(1/\delta))$, implying ultimate concentration of axial current as a current sheet on the inner boundary $r=\delta$.

Second, although not shown in the figure, the energy ${\textstyle\frac{1}{2}}\int b_{z}^{2}\,{\textnormal{d}}V$ of the $b_{z}$-field actually increases during the initial phase (reaching a maximum at $t\approx 18$), but this increase is more than compensated by the decrease of energy of the $b_{\theta}$-field. Later, during the slow diffusive stage, both contributions to energy slowly decrease. Fig. 2(f) shows the decay of energy for $\kappa=0.00001$ as well as for $\kappa=0.0001$; as expected, the diffusive effect is less apparent in the former case, but the initial (non-diffusive) stage is quite similar.

Third, as previously noted, when the Lorentz force is effectively zero (i.e. during the diffusive phase), ${\bf j}=\gamma(r,t){\bf B}$; the coefficient $\gamma$ is then given by

[TABLE]

This coefficient, shown in Fig. 2(g), is far from uniform in $r$, so this is certainly not a Taylor state (for which $\gamma$ would necessarily be uniform). The weak time-dependence of $\gamma$ results from slow continuing evolution during the diffusive stage. This stage is interesting because the field components continuously adjust themselves in such a way that the force-free condition is maintained; in other words, this is not a pure diffusive process, but one that is still constrained to remain nearly force-free through the dynamics encapsulated in eqn. (35).

7 Limiting behaviour as $\eta\rightarrow 0$

In the limit $\eta=0$, i.e. $\kappa=0$, we drop the diffusion terms in (33) and (34), and only the boundary condition $u=0$ of (36) survives, (12) being then automatically satisfied. The two fluxes $\Phi_{z}$ and $\Phi_{\theta}$ are still conserved, within numerical error, as in Fig. 2. Comparison with the situation when $\kappa=0.00001$ is interesting. The difference is in fact very slight up to about $t=100$, but becomes visible, particularly near the boundary $r=1$ by $t=200$, as shown in Fig. 3(a,b). The expanded views in Figs. 3(b,d) show that indeed when $\kappa=0$ the boundary condition $\partial(rb_{\theta})/\partial r=0$ is not maintained at $r=1$; when $\kappa=0.00001$ a weak boundary layer is required, within which the solution adapts to this boundary condition.

As time advances, the curve for $\kappa=0.00001$ continues to evolve due to weak diffusion towards the situation $rb_{\theta}=C$, whereas the curve for $\kappa=0$ remains static. In this case of $\kappa=0$, there is only the initial phase of non-diffusive relaxation to a force-free minimum-energy field, which then remains essentially static.

8 Inclusion of an $\alpha$-effect

It is obvious that purely radial flow cannot lead to local reversal of an axial field that is initially uniform, even in conjunction with diffusion. However, the situation considered by Taylor (1974) in the context of the reversed-field pinch was that of a fully turbulent plasma in which more complex processes must be present. Such turbulence presumably results from persistent instability of the evolving mean field $\bf B$, which has non-zero helicity represented both by the integrals (1) and also by the non-zero pseudo-scalar ${\bf j}\cdot{\bf b}$. Such instabilities must inherit the helicity of this mean field, and may be expected to provide an $\alpha$-effect, such that the (non-dimensionalised) mean electric field becomes

[TABLE]

where $\bf u$ still represents the mean radial flow, and $\kappa$ now includes turbulent as well as molecular diffusivity. A stability analysis following the approach of Furth, Killeen & Rosenbluth (1963) has been carried out by Mizerski (2017), and a resulting anisotropic $\alpha$-effect deduced. Here, we adopt the simpler isotropic prescription

[TABLE]

where $q$ is a pure scalar constant (positive or negative); the factor $\kappa$ is included here in recognition of the diffusive origin of the $\alpha$-effect. The choice (39) of course ensures that $\alpha$ has the same pseudo-scalar character as ${\bf j}\cdot{\bf b}$.

With this prescription for $\alpha$, the boundary condition ${\bf n}\times{\bf E}=0$ becomes

[TABLE]

and, in the cylindrical geometry, this is still obviously satisfied by $j_{\theta}=j_{z}=0$ on $r=\delta,1$. The only modification required is therefore the inclusion of a term

[TABLE]

in the induction equation. In exploratory vein, we include only the z-component444Justification is suggested by study of the interaction of tearing modes proportional to $\exp i(k_{z}z+m\theta)$ (Mizerski 2017, particularly eqn. (44)). It is found that the $\theta$- and $z$- components of the curl of the resulting mean emf differ by a factor proportional to $k_{z}/m$, and that, summing over all Fourier modes, one of these components must vanish., thus replacing (34) by

[TABLE]

while leaving (33) unaltered. From (4) and (5), we then have

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where, for convenience, we include the remaining unaltered evolution equations.

8.1 Field reversal and negative diffusivity

Consider first whether and under what circumstances reversals of the axial field component $b_{z}$ may occur. Suppose that this does happen, and that the reversal is initiated at a critical time $t^{*}$; then $b_{z}(r,t^{*})$ must equal zero for some $r=r^{*}$ where $r^{*}$ may be an internal point or an end-point of the closed interval $[{\delta,1}]$. Moreover, for $t<t^{*}$, $b_{z}(r,t)>0$ for all $r\in[\delta,1]$. Then, at $t=t^{*}$, we have

[TABLE]

and

[TABLE]

since the curvature is necessarily positive at this point. From eqn.(43), we then have

[TABLE]

all other terms vanishing by virtue of (47). Hence it would appear that $b_{z}(r,t)$ will indeed become negative in a neighbourhood of $r^{*}$, provided

[TABLE]

This can therefore occur only if $q>0$ and $b_{\theta}$ is sufficiently strong; and it is likely to occur first in the region where $|b_{\theta}|$ is maximal. More generally, (43) may be written

[TABLE]

The second term on the right-hand side has a diffusive character, but with negative diffusivity in any region where $q\,b_{\theta}^{2}>1$. (The remaining terms of the right-hand side involve only $b_{z}$ and $\partial b_{z}/\partial r$.) It follows that $b_{z}$ can become negative only if this is ‘triggered’ by a period of negative diffusivity in some $r$-interval.

We continue to use the parameter values (27). With these values, $b_{\theta}(r,0)$ is maximal at $r\approx 0.9765$, with maximum value $0.5059$; an interval of negative diffusivity therefore occurs for $q\gtrsim 3.907$.

8.2 Results for $q\leq 4$

Numerical integration for $q\lesssim 3.5$ were quite regular, and not greatly different from the situation when $q=0$. However, as might be expected from the above discussion, numerical instabilities that are difficult to control appear when $q\gtrsim 4$. This however is the regime that must be investigated in seeking possible reversal of $b_{z}(r,t)$.

This led us to adopt a controllable numerical procedure, specifically 4th-order finite-differences in the radial direction and 2nd-order Adams-Bashforth time-stepping with Crank-Nicolson treatment of the diffusive terms. The results for $q=0$ were as expected in complete agreement with those obtained using Mathematica. We focus first on the near-critical situation when $q=4$. Let $\rho(r,t;\,q)$ denote the density computed for any particular value of $q$, and let

[TABLE]

similarly for $\Delta b_{z}(r,t;\,q),\Delta rb_{\theta}(r,t;\,q),\Delta u(r,t;\,q),\,\Delta(\mathbf{j}\times\mathbf{b})_{r}(r,t;\,q)$ and $\Delta\gamma(r,t;\,q)$. Figs. 4 and 5 show numerical results for $\kappa=0.01$; the left-hand columns show curves for $q=0$, while the right-hand columns show the differences $\Delta b_{z}(r,t;\,q)$, etc.

Note first from the last row of Fig. 4 that the negative (pinching) velocity is initially decreased in magnitude by the $\alpha$-effect when $q=4$; it actually becomes weakly positive near $r=1$ for $t\gtrsim 2$ and is positive over the whole range $(\delta,1)$ for $t\gtrsim 2.5$. This implies a corresponding net decrease in the inward transport of mass; however, the decrease in transport for the magnetic field component $b_{z}(r,t)$ is more than compensated by the direct action of the $\alpha$-effect. When $q=4$, $b_{z}(r,t)$ decreases much more rapidly than when $q=0$ (by a factor of about 3) near $r=1$ in the early stage of relaxation; $rb_{\theta}(r,t)$ increases more rapidly near $r=1$, but by a more modest amount ($\sim 6\%$). The function $\gamma(r,t)$ defined by eqn. (37) is also changed by $\sim 25\%$ when $q=4$, but there is no apparent tendency for $\gamma(r,t)$ to become more uniform.

8.3 Results for $q=5.5$

As indicated above, we anticipated numerical problems for $q\gtrsim 4$, and we did indeed run into these. Typically, in the range of $q$ between 4 and 6, a packet of oscillations in $b_{z}$ of very short wavelength appears at $t=0+$ in the region of negative diffusivity near $r=1$. These oscillations move inwards, in tandem with $rb_{\theta}$ which decreases till the local diffusivity $\kappa(1-qb_{\theta}^{2})$ becomes positive, at which stage the oscillations in $b_{z}$ are damped out, the subsequent evolution being quite smooth. Fig. 6 shows this subsequent evolution for $q=5.5$. Note that already at the early time $t=0.05$, $b_{z}$ is negative near $r=1$, presumably a consequence of the early negative diffusivity in this region. The field then relaxes back, becoming positive at $r=1$ by time $t=0.15$.

The short period of reversed $b_{z}$ near $r\!=\!1$ is interesting in the context of the reversed-field pinch. However, we can’t be certain that this is a genuine physical effect rather than just a consequence of adopting a unphysical model for the $\alpha$-effect yielding a period of negative diffusion. The behaviour for $t=0+$ is evidently non-analytic; an asymptotic treatment of the behaviour as $t\downarrow 0$ is presented in Appendix A.

9 Conclusions

We have investigated the relaxation of an axisymmetric magnetic field having both axial and toroidal components in a cylindrical geometry with perfectly conducting boundaries. The density is assumed very small and fluid pressure is neglected compared with magnetic pressure. A purely radial flow is driven by the Lorentz force and energy is dissipated by viscosity. In the zero-resistivity limit, the field rapidly relaxes to a force-free state. When weak resistivity is taken into account, the initial rapid relaxation is followed by slow decay of the field which is constrained to remain nearly force-free with ${\bf j}=\gamma(r,t){\bf B}$. However $\gamma(r,t)$ is quite strongly non-uniform, so this is not a Taylor state.

In §8, we have explored the possibility that an $\alpha$-effect, with $\alpha$ proportional to ${\bf j}\cdot{\bf b}$, might be capable of causing axial field reversal near the outer boundary where the toroidal field component $b_{\theta}$ is initially strong. We have found that if this $\alpha$-effect is sufficiently strong, it can produce a region of negative effective diffusivity of the $b_{z}$-field near the outer boundary; this can instantaneously generate high-frequency short-wavelength oscillations which are rapidly damped as they move into the interior region of positive diffusivity. A transitory reversal of $B_{z}$ occurs near the outer boundary during this process.

There are serious numerical difficulties in handling such a situation; nevertheless, this work points to one possibility whereby a reversed axial field, as observed in the reversed field pinch (Taylor 1974), can be dynamically generated from an initially uniform axial field, through the combined action of pinching by the $b_{\theta}$-field and a suitably contrived $\alpha$-effect. Of course, it would be desirable to derive a correct form of this $\alpha$-effect, through investigation of the turbulence that results from instabilities of the relaxing field. Work is ongoing on this aspect of the problem.

Acknowledgments. The partial funding of the Ministry of Science and Higher Education of Poland within the grant no IP 2014 031373 and statutory activities No 3841/E-41/S/2015 is gratefully acknowledged. This work was also partially funded from the Leading National Research Centre (KNOW) received by the Centre for Polar Studies in Poland for the period 2014-2018.