2.1. Magnetic surfaces

We can evade the constraint (2.4) if we prevent the lines of force closing after one circuit of the torus. A simple way to do this – which is the basis for the tokamak and the toroidal pinch – is to induce a toroidal current (see figure 4). This produces a poloidal field so that, after one circuit around the torus, a line of force is displaced poloidally from its starting point. Then (hopefully) if the line of force is followed indefinitely around the torus it will generate a closed ‘toroidal magnetic surface’ (or ‘flux surface’), as shown in figure 5. The dots indicate successive intersections of a single field line with a plane at a fixed toroidal position, known as a Poincaré plot. The average angle, ${\it\iota}$ , between successive intersections is the ‘rotational transform’ and the inverse of ${\it\iota}$ is the tokamak ‘safety’ factor, $q=2{\rm\pi}/{\it\iota}$ . (This is also the ratio of the number of circuits a field line makes in the toroidal direction to the number it makes in the poloidal direction.)

Figure 5. Poincaré plot generated by following a magnetic field line through many circuits around the torus.

Since $\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}p=0$ , we see that $p$ must be constant over a flux surface. Therefore, what we really require for confinement is a set of nested flux surfaces (see figure 6). Each magnetic surface is defined by ${\it\Psi}(x,y,z)=\text{constant}$ , where $\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\Psi}=0$ and $p=p({\it\Psi})$ . (The degenerate surface at the centre, where we want maximum plasma pressure, is the ‘magnetic axis’.)

Figure 6. Nested magnetic flux surfaces.

At this point we should consider what has become of the constraint (2.4) when we have flux surfaces rather than closed field lines. There is a very useful trick for dealing with a switch from closed lines of force to a closed toroidal flux surface. Consider an integral along a line of force such as

Now introduce a thin flux tube of cross-section ${\it\delta}A$ , with flux ${\it\delta}F=B\,{\it\delta}A$ , enclosing this field line. Then

Here the line integral has been converted to a volume integral. If the field line now sweeps out a surface, then the thin flux tube enclosing it sweeps out a thin shell between two adjacent toroidal flux surfaces, as shown in figure 7. So we effectively replace

where $V({\it\Psi})$ is now the volume within a toroidal flux surface.

Figure 7. Shell swept out by thin flux tube surrounding flux surface.

Applying this result to the constraint (2.4), it becomes

which is automatically satisfied since $\boldsymbol{{\rm\nabla}}p$ is normal to the toroidal flux surface (i.e. parallel to $\text{d}\boldsymbol{S}$ ).

We should also see what becomes of the constraint (2.4) when the magnetic field is non-axisymmetric. For this we write the vector potential (as we can any vector) as $\boldsymbol{A}={\it\alpha}\times \boldsymbol{{\rm\nabla}}{\it\beta}+\boldsymbol{{\rm\nabla}}{\it\phi}$ , so $\boldsymbol{B}=\boldsymbol{{\rm\nabla}}{\it\alpha}\times \boldsymbol{{\rm\nabla}}{\it\beta}$ . The functions ${\it\alpha}(x,y,z)$ and ${\it\beta}(x,y,z)$ , known as Clebsch variables, are constant along a field line and the pressure is a function $p=p({\it\alpha},{\it\beta})$ . Then we have

and (2.4) becomes

Remembering that ${\it\alpha}$ and ${\it\beta}$ are constant along $\text{d}l$ , this implies

where

Hence, in a toroidal system with closed, but not necessarily purely toroidal, field lines, equation (2.4) imposes the constraint $p\equiv p(U)$ .

In conclusion, therefore, the consequences of the equilibrium force balance $\boldsymbol{j}\times \boldsymbol{B}=\boldsymbol{{\rm\nabla}}p$ in a torus are:

1. (i) If field lines generate flux surfaces ${\it\Psi}$ , then $p=p({\it\Psi})$ .

2. (ii) If field lines close, $p$ is constant along $\boldsymbol{B}$ and $p=p(U)$ , where $U=\oint \text{d}l/B$ .

This means that, in a non-axisymmetric torus, there is a potential conflict between the constraints $p=p({\it\Psi})$ and $p=p(U)$ at the so-called ‘rational surfaces’ where the rotational transform ${\it\iota}=2{\rm\pi}m/n$ with $m$ and $n$ integers and the lines of force close after $n$ circuits. Mathematically, these rational surfaces have zero measure compared to the irrational surfaces. Obviously, the simple fluid model is inadequate to discuss such fine details. Nevertheless, low-order rational surfaces (where $m$ and $n$ are small) play an important role in the stability of toroidal plasmas.

2.2. Magnetic surfaces and the Grad–Shafranov equation

We saw that for good containment we need an equilibrium with closed toroidal magnetic surfaces. But do solutions of $\boldsymbol{j}\times \boldsymbol{B}=\boldsymbol{{\rm\nabla}}p$ with closed magnetic surfaces really exist? In general this is a difficult question, related to the Kolmogorov–Arnold–Moser (KAM) theorem and mathematical chaos, and further complicated by the conflict mentioned above. However, solutions do exist in an axisymmetric torus and in other configurations with an ignorable coordinate (Grad & Rubin Reference Grad and Rubin1958).

2.3. Rotational transform

An important property of a toroidal magnetic surface is its rotational transform ${\it\iota}({\it\Psi})$ , or safety factor $q=2{\rm\pi}/{\it\iota}$ . Indeed, a rotational transform is essential for toroidal confinement (see § 2.1). In general,

where the limit is as the number of circuits $\rightarrow \infty$ . The simplest configuration in which a rotational transform appears is a circular cylinder, with periodicity $2{\rm\pi}R$ , to simulate a torus. In this case each field line is a uniform spiral about the $z$ axis, $r\,\text{d}{\it\theta}/\text{d}z=B_{{\it\theta}}/B_{z}$ . The magnetic surfaces are circular cylinders and

where $I_{z}$ and ${\it\psi}_{z}$ are the current and flux within the surface. Note that this does not depend on the presence of plasma – the current $I_{z}$ may as well flow in a thin wire along the $z$ axis.

Of course, this is trivial; everyone knows an axial current generates an encircling field! What is not trivial is that if we add an externally generated helical field (to resemble a straight stellarator) then a rotational transform exists without any axial current!

One can be forgiven if initially this appears to contradict Ampere’s law, for if we take a circular loop integral around the $z$ axis, $\oint rB_{{\it\theta}}\,\text{d}{\it\theta}$ must vanish when there is no axial current. This paradox is resolved if we note that the helical field modulates $B_{{\it\theta}}$ , so that a loop round a flux surface has fluctuating $r$ and $B_{{\it\theta}}$ . To first order these fluctuations cancel, but a second-order effect remains.

A detailed calculation (Johnson et al. Reference Johnson1958; Spitzer Reference Spitzer1958) shows that a helical field ${\sim}B_{hel}\cos (l{\it\theta}+kz)$ generates a rotational transform ${\it\iota}\sim r^{2l-4}(B_{hel}^{2}/B_{tor}^{2})[2(l-1)+k^{2}r^{2}]$ near the magnetic axis. The transform therefore differs according to whether $l=1$ (when it vanishes on the axis and increases slowly with minor radius $r$ ), $l=2$ (when the transform is finite near the axis and increases slowly with $r$ ) or $l\geqslant 3$ (when the transform vanishes on the axis and increases rapidly with $r$ ). These differences are important for plasma stability and $l=3$ is usually a preferred value – because shear in the rotational transform generally enhances stability.

Surprisingly, neither axial current nor helical windings are necessary for a rotational transform! Indeed, the original suggestion for generating a rotational transform (Spitzer Reference Spitzer1951) was to take a simple toroidal solenoid, with no additional currents or fields, and to ‘twist’ it into a figure-of-eight (see figure 10). In fact, any solenoidal field with a twisted magnetic axis generates a rotational transform. This is purely a geometrical effect and ${\it\iota}$ is the integral of the torsion around the (closed) magnetic axis, independent of minor radius.

Figure 10. Sketch of figure-eight stellarator. (Reproduced with permission from Spitzer (Reference Spitzer1958). Copyright 1958, AIP Publishing LLC.)

3.1. Chaotic regions

We have seen that a single helical perturbation of circular surfaces produces islands at a resonant surface, but these are still magnetic surfaces, albeit more complex ones. This is because the perturbed system retains an ignorable coordinate $(m{\it\theta}+kz)$ . Note that this is true only of the straight helical system; it is not true if we add a helical field to a torus (as in a toroidal stellarator!).

When we impose two or more perturbations of different helical pitch there is no longer an ignorable coordinate (Rosenbluth et al. Reference Rosenbluth1966). If the amplitudes of the two perturbations are small, numerical computations show that each perturbation produces a separate chain of islands at its own resonant surface – but the interaction between them produces small chaotic regions near the X-points of each chain. As the perturbations increase, these chaotic regions expand to fill the space between the two chains.

This behaviour is illustrated in a plane slab model with two incommensurable Fourier perturbations (Wesson Reference Wesson2004). The fields in this model are

and the field lines are determined by integrating

The resulting Poincaré plots at fixed $z$ are shown in figure 12. As each chaotic zone is generated by following a single field line, plasma can readily flow throughout a chaotic zone. An interesting question is therefore: What is the magnitude of the perturbations that cause the chaotic zones to first extend from one island chain to the next? This can be estimated by calculating the width of the islands in each chain, according to the formula given above, and assuming that the last magnetic surface between island chains disappears when their individual widths overlap. However, as we remarked earlier, an accurate calculation is a fundamental problem in mathematical chaos theory.

Figure 12. Illustration of the break-up of magnetic surfaces in the presence of multiple perturbations with increasing strength: (a)  ${\it\delta}B_{y}=0.1$ ; (b)  ${\it\delta}B_{y}=0.2$ ; (c)  ${\it\delta}B_{y}=0.3$ ;(d)  ${\it\delta}B_{y}=1.0$ .

3.2. Mappings

One way to study chaotic break-up of surfaces is by numerical calculation of field lines and construction of the Poincaré plot as in the figures above. However, as field lines must be followed very accurately for many circuits of the torus, such calculations are lengthy. Therefore one often uses a simple ‘Iterative Mapping’ to illustrate the break-up of magnetic surfaces (Lichtenberg & Lieberman Reference Lichtenberg and Lieberman2010).

The mapping

simulates successive intersections of a field line with the Poincaré surface of section without consideration of the path between intersections. Such mappings can easily be rapidly and accurately iterated many times.

In order to represent $\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{B}=0$ , the mapping must be area-preserving. That is, any area must map into an equal area. A general area-preserving map can be written implicitly in the form

where $G\equiv G(x_{n+1},{\it\theta}_{n})$ .

The simplest map that simulates the formation of magnetic islands and chaotic regions is the standard or Chirikov–Taylor map (Lichtenberg & Lieberman Reference Lichtenberg and Lieberman2010),

This corresponds to

The important features of this mapping are:

1. (i) it is periodic in ${\it\theta}$ , representing the poloidal angle;

2. (ii) when $k=0$ , $x_{n}=$ constant, representing the behaviour of a line of force lying on a magnetic surface;

3. (iii) the difference $({\it\theta}_{n+1}-{\it\theta}_{n})$ represents the rotational transform ${\it\iota}$ on a flux surface (note that this varies from one surface to another); and

4. (iv) the term $k\sin {\it\theta}$ represents a perturbation of the magnetic field.

In a Poincaré plot generated by a continuous field line between intersections of the Poincaré plane, $k\sin {\it\theta}$ would represent perturbation by a single helical mode. As we have seen, a single helical perturbation does not create chaotic regions, but because the mapping is discontinuous, the term $k\sin {\it\theta}$ is equivalent to a discontinuous perturbation of the field containing many helical modes.

The results of iterating the standard map are shown in figure 13 for several values of the perturbation amplitude $k$ . As expected, when $k$ is small, islands appear with small chaotic regions near X-points. These chaotic regions expand, at the expense of continuous surfaces, as $k$ increases until the last continuous surface between island chains disappears. By the overlap criterion above, this occurs at $k_{c}\sim {\rm\pi}^{2}/4$ ; the true value is 0.971635…(Greene Reference Greene1968; Lichtenberg & Lieberman Reference Lichtenberg and Lieberman2010).

Figure 13. Result of the standard map (3.15)–(3.16) for varying $k$ : (a)  $k=0.5$ ; (b)  $k=1.0$ ; (c)  $k=2.5$ ; (d)  $k=4.0$ .

4.1. Fluid model with resistivity

The resistive fluid model is described by

and the evolution of the field is given by

The two important time scales in a plasma can be seen in the behaviour of small-amplitude waves, transverse to $\boldsymbol{B}_{0}$ in an incompressible plasma (Alfvén waves). Then $b\sim \exp (-\text{i}{\it\omega}t+\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{r})$ and (4.1) reduce to

giving

Equation (4.5) defines the characteristic speed $v_{A}^{2}=B_{0}^{2}/{\it\mu}_{0}{\it\rho}_{0}$ of Alfvén waves, the Alfvén time scale ${\it\tau}_{A}^{2}\sim a^{2}{\it\mu}_{0}{\it\rho}_{0}/B_{0}^{2}$ and the resistive time scale of the wave damping ${\it\tau}_{R}\sim {\it\mu}_{0}a^{2}/{\it\eta}$ .

The ratio of the resistive to Alfvén time scales is the Lundquist number, $S$ ,

In high-temperature laboratory plasmas $S\sim 10^{6}{-}10^{8}$ (in astrophysical plasmas it is ${\sim}10^{14}$ ). These large values of $S$ create serious difficulties but can also lead to great simplification. The difficulties lie in computing over such an enormous range of time scales. The simplification is that field diffusion can often be ignored and the resistivity set to zero. This leads to some crucial properties of an ideal plasma, or of a real plasma on a time scale much less than ${\it\tau}_{R}$ .

(i) Flux conservation. The magnetic flux ${\it\phi}=\int _{S}\boldsymbol{B}\boldsymbol{\cdot }\text{d}\boldsymbol{A}$ through an area $S$ bounded by a closed curve $C$ that moves with the fluid is constant.

The rate of change of the flux is

where the first term is due to the change in the field $\boldsymbol{B}$ and the second is due to motion of the boundary moving with the fluid. The total change can be written as

which, by virtue of the field evolution equation (4.2), is always zero when ${\it\eta}=0$ .

(ii) Field line preservation. A line moving with the fluid that is initially a line of force always remains a line of force. (Alternatively, if fluid particles initially lie on a common line of force then they continue to do so.)

Figure 14. Magnetic island formation.

To show this, it is convenient to reintroduce Clebsch variables as in § 2.1. Then

where ${\it\alpha}(x,y,t)$ and ${\it\beta}(x,y,t)$ are constant along a line of force (crucially they label each line of force). Therefore

However, we also have

where $\boldsymbol{v}$ is the fluid velocity. Although (4.10) and (4.11) do not completely determine $\partial _{t}{\it\alpha}$ and $\partial _{t}{\it\beta}$ , a consistent solution is

That is to say, there is a labelling of the lines of force in a moving ideal plasma such that each line moves with the fluid velocity $\boldsymbol{v}_{F}=\boldsymbol{E}\times \boldsymbol{B}/B^{2}$ . Because the field lines are locked into the fluid, they cannot break or intersect one another, since this would require a discontinuous fluid velocity.

Note that the choice of field line labels, and of their velocity, is not unique. In fact, instead of saying that field lines move with the velocity of the fluid, $\boldsymbol{v}$ , we could equally say that they move with velocity $(\boldsymbol{v}+\boldsymbol{u})$ , where $\boldsymbol{u}\times \boldsymbol{B}=\boldsymbol{{\rm\nabla}}G$ , and $G$ is any function satisfying $\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}G=0$ . Of course, no physical quantity would be affected by this. It would simply be a different labelling of the imaginary lines of force.

4.2. Effect of field line preservation

The effect of field line preservation in an ideal plasma can be seen if one tries to calculate the growth of magnetic islands in the so-called ‘sheared slab’ configuration $\boldsymbol{B}_{0}=(0,sxB,B)$ , shown in figure 14. In this configuration $B_{y}$ varies only in the $x$ direction and changes sign at $x=0$ . The question is: Can the field lines in this equilibrium reconnect across the $x=0$ plane to form magnetic islands? This would occur if there were an unstable perturbation with $B_{x1}\neq 0$ at $x=0$ . For example, a single Fourier mode ${\sim}\exp (\text{i}ky)$ could create a chain of islands like those in figure 14. (For such a mode, the ‘resonant surface’ $\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{B}_{0}=0$ is located at $x=0$ .)

A perturbation ( $\boldsymbol{B}_{1}$ , $\boldsymbol{v}_{1}$ ) of the sheared slab equilibrium satisfies

and for the single Fourier mode $B_{1}\sim B_{1}(x)\exp (\text{i}ky)$ ,

We see from this that growth of a reconnecting field across the $x=0$ plane would require an infinite fluid velocity at $x=0$ if the plasma resistivity were precisely zero!

5.1. Resistive instabilities: the tearing mode

The basic ‘tearing mode’ instability involves the growth of magnetic islands in the sheared slab that we considered in the previous section (Furth et al. Reference Furth, Killeen and Rosenbluth1963). As we found, if ${\it\eta}=0$ any such growth is prevented by the field line preservation constraint, but now we wish to calculate the rate of growth when ${\it\eta}\neq 0$ . In particular, we are interested in instabilities that grow rapidly compared to resistive diffusion, but slowly compared to MHD Alfvén time scales.

We again start from a sheared slab equilibrium $\boldsymbol{v}=0$ and $\boldsymbol{B}_{0}=(0,B_{0y}(sx),B_{0z})$ . We assume that $B_{0y}(sx)$ is symmetric about $x=0$ and that the plasma is incompressible. Then the growth of a linear perturbation ( $\boldsymbol{B}_{1}$ , $\boldsymbol{v}$ ) is given by

and

It is convenient to write $\boldsymbol{B}_{1}=\hat{z}\times \boldsymbol{{\rm\nabla}}A_{1}$ . Then for a Fourier mode $A\sim A_{1}(x)\exp ({\it\gamma}t+\text{i}ky)$ , $B_{x1}=-\text{i}kA_{1}$ , and we can extract equations for ( $A_{1}$ , $v_{x}$ ) alone, decoupled from all other components,

and

We saw that, when ${\it\eta}=0$ , island growth requires $v_{x}\rightarrow \infty$ at $x=0$ . Hence we expect the effect of resistivity to be most important near $x=0$ . We therefore divide the calculation into a thin inner layer around $x=0$ , of thickness ${\it\delta}$ , where resistivity is important, and two outer regions $|x|>{\it\delta}$ , where resistivity can be neglected. A solution is then constructed in each outer region and matched to the solution in the inner layer.

Because we expect the tearing mode to grow slowly compared to the Alfvén time, we can neglect inertia in the outer regions. Then in these regions

Note that this is exactly the equation for ideal plasma equilibrium in the outer regions – this equilibrium merely adjusts itself to the slowly changing inner solution.

We can integrate (5.5) inwards from $\pm$ infinity, or from the plasma boundaries, and adjust the (linear) amplitude of the left and right solutions to be equal at $x=0$ . Then in general we will find that there is a discontinuity in the slope $(\text{d}A/\text{d}x)$ between the left and right solutions at $x=0$ . We define this discontinuity as the ‘stability index’,

This ${\it\Delta}^{\prime }$ is all we require from the outer solutions and we will see later that ${\it\Delta}^{\prime }>0$ implies instability.

We now need to find a solution in the inner region. In this region $x$ is small but $v_{x}$ and $A_{1}$ vary very rapidly with $x$ . We can therefore neglect the last term in (5.4) as well as derivatives in the $y$ direction. Then, in the thin resistive inner layer, equations (5.3) and (5.4) reduce to

We now ‘stretch’ the thickness of the inner layer by introducing ${\it\tau}=x/{\it\delta}$ , and select

so that the layer has unit width in ${\it\tau}$ . We also rescale the velocity by introducing $w=(sB_{0}{\it\delta}/{\it\gamma})v_{x}$ . Then the resistive layer equations become

This essentially completes the solution! The change in $\text{d}A/\text{d}{\it\tau}$ across the inner layer will be of order ${\it\gamma}{\it\mu}_{0}{\it\delta}^{2}/{\it\eta}$ . Hence, when we match the change in $\text{d}A/\text{d}x$ to the discontinuity ${\it\Delta}^{\prime }$ in the outer solutions, we will have

Therefore ${\it\Delta}^{\prime }=0$ is the stability threshold and, when ${\it\Delta}^{\prime }>0$ ,

In terms of Alfvén and resistive times, ${\it\tau}_{A}=a\sqrt{{\it\mu}_{0}{\it\rho}}/B_{0}$ and ${\it\tau}_{R}=a^{2}{\it\mu}_{0}/{\it\eta}$ , defined relative to a characteristic dimension $a$ of the outer plasma, equation (5.13) becomes

Consequently ${\it\gamma}{\it\tau}_{R}\sim S^{2/5}\gg 1$ and ${\it\gamma}{\it\tau}_{A}\sim S^{-3/5}\ll 1$ . So our assumptions about the tearing mode are confirmed: when resistivity is small, the growth rate is indeed rapid compared to resistive diffusion but slow compared to the Alfvén speed, and the resistive layer is thin, with the width ${\it\delta}$ given by (5.9). Note that in the limit ${\it\eta}\rightarrow 0$ we recover an ideal plasma, but with a discontinuity in the form of a current sheet at the resonant layer $\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{B}_{0}=0$ – a possibility we did not consider in our discussion of ideal plasmas.

It is important to recognise the significance of tearing modes. Of course, we expect that introducing a small resistivity will perturb the motion of a previously ideal plasma. But this is not what tearing modes do. Tearing modes are new modes involving critical layers in which resistivity is dominant no matter how small it may be. Mathematically they arise because the differential equations governing resistive plasma are of higher order than those governing ideal plasma. This means that there are solutions for the resistive plasma that do not exist when resistivity is strictly zero. (This is very similar to the effect of viscosity in fluid dynamics, where it also introduces additional solutions and is dominant in critical boundary layers.)

In addition to its intrinsic interest, there are several reasons why the calculation of tearing modes is important: (a) it describes an instability that has been seen in experiments, (b) it illustrates the importance of small departures from ideal Ohm’s law, (c) it confirms that reconnection can occur much faster than resistive diffusion and (d) the matching procedure described in the calculation is used in many other calculations when there is a small departure from the ideal plasma. In these calculations more complicated models are used in the resonant layer, but the outer region calculation is essentially unchanged.

In one respect, however, our calculation of the tearing mode growth rate is atypical. Recall that in the outer region

We assumed that $j_{0}$ was symmetric about $x=0$ , so the right-hand side of this equation tends to a constant times $A$ as $x\rightarrow 0$ and the general solution near $x=0$ is

In this case $\text{d}A/\text{d}x$ and $A$ are well behaved at $x=0$ and there is no difficulty in calculating ${\it\Delta}^{\prime }$ . But if $j_{0}$ is not symmetric about $x=0$ , the right-hand side of the equation ${\sim}A/x$ as $x\rightarrow 0$ . Then the general solution near $x=0$ is

so that $A$ and $\text{d}A/\text{d}x$ are divergent at $x=0$ . In this case, rather than matching the derivatives of inner and outer solutions, we match the large ( ${\it\alpha}$ ) and small ( ${\it\beta}$ ) solutions in the inner and outer regions. (If we had included plasma pressure in our calculation, the general solution would have been of the form ${\it\alpha}|x|^{q}+{\it\beta}|x|^{1-q}$ near $x=0$ , but we again match the large and small solutions.)

One might reasonably question whether at this level of detail the plasma can really be represented by a fluid model. However, the underlying concept of reconnection occurring in a thin layer, but dependent on a parameter ${\it\Delta}^{\prime }$ calculated from the plasma outside this layer, appears to be robust. This view is reinforced by another interpretation of ${\it\Delta}^{\prime }$ . Suppose we remove the resistive layer and replace it by a thin flexible insulating membrane. The magnetic pressure on this membrane is $B^{2}/2{\it\mu}_{0}$ , so if we deform it by a tearing-mode-like displacement ${\it\xi}_{x}$ , the work done by the plasma (per unit area) on one side of the membrane is

As there is no change of plasma volume, the contribution from $B_{0}^{2}$ is zero, and from the earlier calculation of tearing modes, in the ideal plasma

where ${\it\gamma}{\it\xi}_{x}=v_{x}$ . Therefore the net change of energy is

Hence ${\it\Delta}^{\prime }$ is a measure of the energy available in the ideal plasma if the field line conservation constraint is relaxed at the resonant surface $\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{B}_{0}=0$ .

5.2. Forced reconnection

As we saw, tearing mode reconnection is driven by the energy released from an unstable plasma. However, reconnection may occur even more rapidly when it is ‘forced’, that is, driven by plasma flows or imposed on the reconnecting region. Conceptually we might view the forcing as being applied externally, as in figure 15, effectively increasing the value of ${\it\Delta}^{\prime }$ . In practice we will be more interested in reconnection driven by turbulence.

Figure 15. Faster reconnection: $P$ represents the forcing perturbation at the wall.

An often quoted idealised form of forced reconnection is the Sweet–Parker model. This represents plasma flow through an X-point of a two-dimensional magnetic field, as illustrated in figure 16. The antiparallel reconnecting fields are forced towards the X-point by flows in the $\pm x$ direction, reconnect, and then flow outwards along the $\pm y$ directions. We can readily imagine that the forced inflows cause ‘flattening’ of the X-point to eventually form a current sheet of thickness ${\it\delta}$ and length $l$ , as indicated in figure 16, across which the field $B$ changes direction. We can then envisage a steady state, with plasma forced into the current sheet across the field, at a velocity $v_{x,in}$ controlled by resistivity, and flowing out along the field at a much higher velocity,  $v_{y,out}$ .

Figure 16. Model of forced reconnection.

From Ohm’s law we have for the inflow, $v_{in}\sim {\it\eta}/{\it\delta}{\it\mu}_{0}$ . Pressure balance, ${\it\rho}v_{y}^{2}/2=B_{0}^{2}/2{\it\mu}_{0}$ , gives for the outflow, $v_{out}\sim B/\sqrt{{\it\mu}_{0}{\it\rho}}$ , and mass conservation requires that $v_{out}{\it\delta}\sim v_{in}l$ . Consequently, $v_{in}^{2}\sim l^{2}/({\it\tau}_{R}{\it\tau}_{A})$ so that flux is swept in, and reconnects at a rate ${\sim}({\it\tau}_{A}{\it\tau}_{R})^{-1/2}$ , where ${\it\tau}_{A}$ and ${\it\tau}_{R}$ are again the Alfvén and resistive diffusion times relative to the outer dimensions of the reconnecting region (with scale length ${\sim}l$ ).

Of course, both the tearing mode and the Sweet–Parker model are over-idealised and, except in carefully controlled experiments, can scarcely represent real situations. In particular, the strictly two-dimensional nature of these models is likely to be unrepresentative, and the Sweet–Parker outflow is also likely to be unstable. (Numerical simulations generally show that reconnection is faster than it is in the idealised models.)

In a turbulent plasma, a more realistic picture is one in which elements of plasma, carrying differently oriented frozen-in fields, are forced together. In the absence of reconnection, they could not penetrate each other and a current layer would form between them. However, resistivity allows some penetration and reconnection. If this is not fast enough to keep up with the inflow, magnetic field will pile up against the current sheet, making the Alfvén speed larger and the current sheet thinner, so that reconnection speeds up. Consequently, in a turbulent resistive plasma, we can expect reconnection to proceed freely.

7.1. Relaxation constraints

We have seen that, for a perfectly conducting plasma, important constraints follow from

which we now write as

where $\boldsymbol{A}$ is the vector potential and ${\it\chi}$ is an arbitrary scalar. Consequently, changes in $\boldsymbol{A}$ must satisfy a similar equation,

where ${\bf\xi}$ is an arbitrary plasma displacement (and $\overline{{\it\chi}}$ follows from ${\it\chi}$ ). It is clear that this does not restrict ${\it\delta}\boldsymbol{A}_{\bot }$ , but ${\it\delta}A_{\Vert }$ is constrained by

for any closed field line or magnetic surface (see § 2.1).

A more convenient, and more fundamental, form of this constraint is that, for every infinitesimal flux tube surrounding a closed field line labelled by Clebsch variables ${\it\alpha}$ , ${\it\beta}$ , and for every magnetic flux surface, the volume integral

is invariant. This provides an infinity of constraints – which replace the single constraint of the flexible wire.

We can confirm the invariance of $K_{{\it\alpha}{\it\beta}}$ directly:

where $V$ and $S$ are respectively the volume and surface area of the infinitesimal flux tube. The last term arises from the motion of the flux tube and the others from the local change in $(\boldsymbol{A}\boldsymbol{\cdot }\boldsymbol{B})$ . Then, from (7.2),

and from (7.1),

(recall that $\boldsymbol{v}$ is tangential to the boundary $S$ ).

Hence $\text{d}K/\text{d}t=0$ , and to find the relaxed state we need to minimise the magnetic energy $\int _{V}\text{d}{\it\tau}\,B^{2}/2{\it\mu}_{0}$ over all variations ${\it\delta}\boldsymbol{A}$ subject to the infinity of constraints $K_{{\it\alpha}{\it\beta}}=\text{constant}$ . This leads to the variational problem

This is very inconvenient because of the restriction ( $\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\lambda}=0$ ) on ${\it\lambda}$ that is implicit in the use of Clebsch variables ${\it\alpha}$ and ${\it\beta}$ . We can avoid this by introducing the restriction explicitly through another Lagrange multiplier. Then the variational problem becomes

We also need to consider the appropriate boundary condition. When the magnetic field is entirely contained within a closed conducting shell, the tangential component of ${\it\delta}\boldsymbol{A}$ on the boundary is zero.

Returning to the variational problem, the minimisation requires

From (7.12), we see that

therefore

where the first term is zero. Then (7.11) becomes

Redefining ${\it\lambda}$ , we obtain finally

Now, note well, this cannot be the correct description of the quiescent state, because ${\it\lambda}_{{\it\alpha}{\it\beta}}$ must be determined from the infinity of $K_{{\it\alpha}{\it\beta}}$ , which in turn will reflect every detail of the initial conditions – entirely contrary to the observations.

The resolution of this dilemma is to recognise, as we did earlier, that even small resistivity leads to field line breaking and reconnection. Then it is impossible to identify individual field lines and the individual $K_{{\it\alpha}{\it\beta}}$ are no longer of any significance. But the sum of all the $K_{{\it\alpha}{\it\beta}}$ is still invariant. This is just the total magnetic helicity

(Apart from the fact that they are rendered irrelevant by reconnection, the individual $K_{{\it\alpha}{\it\beta}}$ differ from the total helicity $K_{0}$ in another fundamental way. For the $K_{{\it\alpha}{\it\beta}}$ to be relevant, we would need to calculate and identify all the field lines – a knowledge of the field $\boldsymbol{B}$ itself is not enough – whereas $K_{0}$ is well defined by the field alone. In a loose analogy with particle dynamics, we could say that $K_{0}$ is an isolating integral, which constrains the available configuration space, whereas the $K_{{\it\alpha}{\it\beta}}$ are non-isolating integrals, which do not constrain the evolution of the field.)

We conclude, therefore, that the relaxed state of a slightly resistive turbulent plasma is obtained by minimising the magnetic energy subject to a single constraint, here denoted by the Lagrange multiplier ${\it\mu}$ ,

That is

and upon redefining ${\it\mu}$ , finally

This is no longer any force-free state, depending on the initial conditions; it is a unique state dependent only on a single parameter ${\it\mu}$ . Thus we have already recovered the uniqueness of the quiescent phase.

7.2. Helicity and field line linkage

Before discussing properties of the relaxed states given by (7.21), we should note some features of the magnetic helicity $K_{0}$ .

The helicity $K_{0}=\int \boldsymbol{A}\boldsymbol{\cdot }\boldsymbol{B}\;\text{d}{\it\tau}$ is defined only if the volume of integration is bounded by a closed flux surface. Otherwise it is not gauge-invariant. This reflects the fact that $\boldsymbol{A}\boldsymbol{\cdot }\boldsymbol{B}$ is not a local quantity. One cannot specify where the helicity resides, or a helicity density; only the total helicity within a closed flux surface can be defined.

In fact, helicity is a topological quantity related to the linkage of one flux tube with another. For example, for two flux tubes as in figure 17, the helicity is $K_{0}=2{\it\psi}_{1}{\it\psi}_{2}$ if they are linked as in the figure, and zero if they are not linked (Moffat Reference Moffat1978). Where the helicity is located when they are linked is a meaningless question.

Figure 17. Linked flux tubes.

If the volume concerned is multiply connected, as in a torus, then even though the bounding surface is a flux tube, the simple formula (7.18) for $K_{0}$ is not fully gauge-invariant. The full gauge-invariant form should be

where $C_{1}$ and $C_{2}$ are circuits the long and short way around the torus boundary. However, in most applications the circuit integrals are themselves constant and the simpler form can be used.

7.3. Back to reality

If all this sounds too much like exotic mathematics, there is a down-to-earth engineering interpretation of $K_{0}$ – as the volt-seconds from the driving transformer. (Volt-seconds are a standard measure of transformer performance and something you would certainly know was real if you touched the terminals!)

To see this connection, we introduce a small gap in the toroidal shell, as in figure 18, in order to apply the loop voltage $V_{loop}$ . Then, using (7.22)

But $\partial \boldsymbol{A}/\partial t\neq 0$ only in the gap. There we may write $\text{d}\boldsymbol{S}=\text{d}\boldsymbol{l}\times \text{d}\boldsymbol{s}$ and

where $\oint (\boldsymbol{A}\boldsymbol{\cdot }\text{d}\boldsymbol{s})={\it\Psi}_{tor}$ and $\oint (\partial _{t}\boldsymbol{A}\boldsymbol{\cdot }\text{d}\boldsymbol{l})=-V_{loop}$ so

Therefore, if loop volts are applied at a constant toroidal flux,

Figure 18. Interpreting magnetic helicity in terms of stored volt-seconds.

In principle (7.25) and (7.26) describe the start-up of the toroidal plasma, but as this is ill-controlled and volt-seconds are lost in ionisation and heating, we cannot usefully apply them to this initial stage. (The start-up loop voltage is much larger than that required to maintain the relaxed state.)

7.4. Uniqueness of relaxed states

In the next few lectures we will discuss specific solutions of equation (7.21) for relaxed states. But first we show that each relaxed state is uniquely determined – with no arbitrary or fitted parameters – by the externally controlled quantities ${\it\Psi}_{tor}$ and $K$ . To see this, scale all dimensions in the problem by ${\it\mu}$ . It is then clear, if it was not already, that the solution of (7.21) must be of the form $\boldsymbol{B}=c\,\hat{\boldsymbol{B}}({\it\mu}\boldsymbol{r};{\it\mu}a,{\it\mu}R)$ , $\boldsymbol{A}=(c/{\it\mu})\hat{\boldsymbol{A}}({\it\mu}\boldsymbol{r};{\it\mu}a,{\it\mu}R)$ , where $a$ and $R$ are minor and major radii and $c$ is the arbitrary scale factor (Taylor Reference Taylor1974b ). Then

and

giving

which fixes ${\it\mu}$ . The scale factor $c$ can then be set from $K$ or ${\it\Psi}_{tor}$ .

7.5. Summary

Now let us summarise what we have seen so far. Observations of turbulent plasmas in a toroidal pinch, with initial toroidal field $B_{0}$ and current $I$ , have shown that quiescent states form that depend only on the parameter ${\it\theta}\equiv 2I/aB_{0}$ , and if ${\it\theta}>{\it\theta}_{c}$ the toroidal field spontaneously reverses.

We interpret these as ‘relaxed’ states, that is, states of minimum energy subject to constraints. In an ideal plasma there would be an infinity of constraints, but this leads to a state in strong disagreement with observations. In a slightly resistive plasma, the infinity of individual constraints become irrelevant and the only surviving constraint is that the total magnetic helicity $K_{{\it\alpha}{\it\beta}}=\int _{{\it\alpha}{\it\beta}}\boldsymbol{A}\boldsymbol{\cdot }\boldsymbol{B}\;\text{d}{\it\tau}$ is invariant. This leads to relaxed states in agreement with observations, as we shall see in detail in the following sections.

8.1. Spontaneous field reversal

Since the ‘toroidal’ field is $B_{z}\sim \text{J}_{0}({\it\mu}r)$ , it is clear that the onset of spontaneous field reversal in the relaxed state occurs at ${\it\mu}a>2.4$ . This corresponds to an experimental pinch parameter ${\it\theta}>1.2$ – which again agrees well with the experimental observations.

The experimental data from reverse field pinches (RFPs) is usually presented in the form of an $F{-}{\it\theta}$ diagram, where

and the angled brackets denote the average over the cross-section. The value of $F$ in the relaxed state is

Note that this is calculated from ‘first principles’, and contains no empirical or fitted factors.

Equation (8.7) is in good agreement with experiment, as can be seen in figure 20. Note that the experimental points from several machines lie on a universal curve, close to the theoretical one but at somewhat higher ${\it\theta}$ . (This may be due to the lower current density near the wall, where the plasma is contaminated by impurities.)

Figure 20. An $F{-}{\it\theta}$ diagram showing data from the pinch experiments HBTX1, ALPHA and ZETA and the theoretical curve. (Reproduced with permission from Bodin & Newton (Reference Bodin and Newton1980).)

Even more striking evidence for relaxation can be seen in the time-dependent behaviour shown in figure 21. We see that if the operating point $(F,{\it\theta})$ is forced away from the locus of relaxed states (e.g. by a temporary increase in loop voltage), it very quickly ( $0.8~{\rm\mu}\text{s}$ ! – comparable to the Alfvén time ${\sim}1~{\rm\mu}\text{s}$ ) reverts back to it. However, as long as the operating point is changed slowly, the plasma stays close to the theoretical relaxed-state curve (figures 21 b and 22).

Figure 21. Time-dependent $F{-}{\it\theta}$ curve for HBTX1 in (a) fast mode and (b) slow mode. Time intervals are given in microseconds. (Reproduced with permission from Bodin & Newton (Reference Bodin and Newton1980).)

Figure 22. Time-dependent $F{-}{\it\theta}$ curves from pinches operated in slow mode: (a) ZT-40 experiment (from DiMarco Reference DiMarco, Post, Baldwin and Ryutov1983), (b) TPRE experiment (reproduced with permission from Tamaru et al. Reference Tamaru, Sugisaki, Hayase, Shimada, Hirano, Maejima, Ogawa, Hirano, Kitagawa and Sato1979) and (c) REPUTE experiment (reproduced with permission from Toyama et al. Reference Toyama, Asakura, Hattori, Inoue, Ishida, Matsuzuka, Miyamoto, Morikawa, Nagayama, Nihei, Shinohara, Ueda, Yamagishi, Yoshida, Pocs and Montvai1985).

9.1. Primitive and mixed states

The task of finding which is the lowest-energy solution is made easier by two comparison theorems (Taylor Reference Taylor1968):

(i) The first theorem is that if two solutions have the same value of ${\it\Psi}$ and $K$ , then the one with the smaller ${\it\mu}$ has the lower energy. The proof of this is as follows.

If $\int \boldsymbol{A}_{1}\boldsymbol{\cdot }\boldsymbol{B}_{1}\;\text{d}{\it\tau}=\int \boldsymbol{A}_{2}\boldsymbol{\cdot }\boldsymbol{B}_{2}\;\text{d}{\it\tau}$ with $\boldsymbol{A}_{1}=\boldsymbol{A}_{2}$ on the boundary, and both $\boldsymbol{{\rm\nabla}}\times \boldsymbol{B}_{1}={\it\mu}_{1}\boldsymbol{B}_{1}$ and $\boldsymbol{{\rm\nabla}}\times \boldsymbol{B}_{2}={\it\mu}_{2}\boldsymbol{B}_{2}$ , then

and

Therefore, if $W_{1}$ and $W_{2}$ are the energies of two relaxed states, then

Hence, if ${\it\mu}_{2}^{2}>{\it\mu}_{1}^{2}$ , $W_{2}>W_{1}$ .

(ii) The second theorem is that the eigenvalues of (9.12) for the same $k$ increase with increasing $m$ . The proof of this is similar to the argument above and to the Sturm–Liouville comparison theorems in the theory of differential equations (Ince Reference Ince1956).

By theorem (i), to find the solution with the lowest energy, we must select that with the smallest ${\it\mu}$ . For a Type I solution this means taking the smallest root of (8.5). By theorem (ii), for a Type II solution it means taking $m=1$ and then searching through $k$ at $m=1$ to find the smallest eigenvalue of (9.12). This occurs at $ka\approx 1.25$ and has the value ${\it\mu}a=3.11$ .

This essentially completes the solution of our problem. The primitive solution $\boldsymbol{B}={\it\alpha}\boldsymbol{B}^{00}$ is the correct relaxed state so long as ${\it\mu}a<3.11$ , i.e.  ${\it\theta}\leqslant 1.6$ . Otherwise, the correct state is a mixed, helical state, $\boldsymbol{B}={\it\alpha}_{0}\boldsymbol{B}^{00}+{\it\alpha}_{mk}\boldsymbol{B}^{mk}$ , with $m=1$ , $k=1.25$ and ${\it\mu}a=3.11$ .

Recalling that in the primitive state ${\it\mu}a=F(K/{\it\Psi}^{2})$ , this is equivalent to saying that the primitive state is the correct lowest-energy state for $K/{\it\Psi}^{2}$ less than some critical value, and the mixed state is the correct one for $K/{\it\Psi}^{2}$ greater than this value. This situation is illustrated schematically in figure 23. Equally, of course, we could say that the primitive state occurs when the volt-seconds are below a critical value and the mixed state occurs when the volt-seconds exceed that value.

Figure 23. Illustration of the energy of the relaxed states in a cylindrical plasma as a function of $K/{\it\Psi}^{2}$ .

9.2. Current saturation and spontaneous symmetry breaking

The Type II helical states with ${\it\mu}a=3.11~({\it\theta}\sim 1.6)$ have a remarkable property. In these states, since ${\it\mu}$ is fixed, an increase in volt-seconds (at fixed ${\it\Psi}$ ) produces no increase in the toroidal current $I$ ! In other words, a helical relaxed state acts as a current-limiting circuit element! When the voltage is increased, instead of producing an increase in current, the current remains the same but flows in a tighter helix with higher inductance. (This is reminiscent of the flexible wire with which we began the discussion of relaxed states!)

Evidence for this second type of relaxed state, and for current saturation, is shown in figure 24. This shows discharges in which ${\it\theta}$ is initially driven to a large value. However, it then falls quickly to ${\sim}1.6$ and remains there for the rest of the discharge. (Some discharges can be prevented from relaxing and held at ${\it\theta}>1.6$ , but at the cost of increased fluctuation levels and power losses; see figure 25.)

Figure 24. Evidence for second critical ${\it\theta}$ in HBTX1. (Reproduced with permission from Bodin & Newton (Reference Bodin and Newton1980).)

Figure 25. Driving a toroidal pinch current above the critical value. Fluctuation level versus ${\it\theta}$ in ZT-40. (Reproduced with permission from Watt & Nebel (Reference Watt and Nebel1983). Copyright 1983, AIP Publishing LLC.)

The transition from a primitive axisymmetric state to a helical state means that relaxation theory predicts not one, but two critical values of the pinch parameter ${\it\theta}$ :

1. (i) at ${\it\theta}\approx 1.2\rightarrow$ spontaneous field reversal occurs, but the plasma remains axisymmetric (see § 8.1);

2. (ii) at ${\it\theta}\approx 1.6\rightarrow$ current saturation occurs and a helical deformation sets in (see § 9.1).

Note that the second transition, from an axisymmetric state to the helical mixed state, is a neat example of ‘spontaneous symmetry breaking’.

9.3. Stability

Like all relaxed states, the relaxed states of the toroidal pinch are ‘minimum energy’ and therefore stable against any perturbations that do not change the helicity. This, of course, includes all ideal MHD perturbations.

The transition from Type I to Type II (helical) relaxed states at ${\it\mu}a=3.11$ coincides with a linear stability boundary for resistive tearing modes. However, the present analysis is not restricted to the linear regime. The nonlinear amplitude of the helical deformation is given precisely in terms of $K/{\it\Psi}^{2}$ as (Martin & Taylor Reference Martin and Taylor1974; Taylor Reference Taylor1986)

10.1. Nature of the eigenfunctions in multipinch

In the same way that we had to determine which of many eigenfunctions in a cylinder had the lowest energy, we still have to decide whether the lowest eigenvalue in the multipinch is in fact associated with an axisymmetric eigenfunction. In a cylinder the lowest eigenfunction was helical, but in a re-entrant cross-section the reverse may be the case.

Figure 28. (a) Idealised multipinch and (b) toroidal wavenumber corresponding to the lowest eigenvalue as a function of ${\it\Delta}$ . (Reproduced with permission from Gimblett et al. (Reference Gimblett, Hall, Taylor and Turner1987). Copyright 1987, AIP Publishing LLC.)

To see why, consider the extreme case of a multipinch made up of two circular tori, one above the other, connected only by a narrow slit where they make contact. In this case an axisymmetric eigenfunction, carrying no net flux, could be constructed from a circular cylinder solution in each torus, one carrying positive toroidal flux and the other carrying equal negative flux. To be continuous, the toroidal field would then have to vanish at the slit – so the two cylinder solutions must have ${\it\mu}a\sim 2.4$ , much lower than the helical eigenvalue ${\it\mu}a\sim 3.11$ . This suggests that the lowest eigenfunction in a multipinch is axisymmetric, at least when the ‘waist’ is narrow. (It also explains why the field reversal point is so close to the zero-flux eigenvalue.)

Numerical computations (Gimblett et al. Reference Gimblett, Hall, Taylor and Turner1987) confirm this picture. Figure 28 shows the toroidal wavenumber of the lowest eigenfunction for an idealised linear multipinch, with cross-section made up of circular arcs. When the ‘waist’ ${\it\Delta}$ is comparable with the diameter of the lobes ( ${\it\Delta}/a\sim 1$ so there is essentially no indentation at the waist), the helical wavelength $2{\rm\pi}/k_{c}$ is ${\sim}10$ times the lobe width. As the waist is narrowed, the wavelength increases and becomes infinite, that is, the eigenfunction becomes axisymmetric, when ${\it\Delta}=0.66$ . (This is somewhat larger than the actual experiment, but its cross-section is not made up of circular arcs.)

We can now describe the properties we would expect of relaxed states in the multipinch.

1. (i) For small $K/{\it\Psi}^{2}$ (low volt-seconds), the relaxed state will be axisymmetric and symmetric about the equator/mid-plane, as in figure 27(a). In this state, plasma current is the same in both lobes and increases with increasing volt-seconds.

2. (ii) When $K/{\it\Psi}^{2}$ , or the volt-seconds, reach a critical value, the relaxed state becomes a superposition of the axisymmetric state and the eigenfunction shown in figure 27(b). Therefore it continues to be axisymmetric but is no longer symmetric about the mid-plane. As the volt-seconds (at fixed toroidal flux) increase further, more current flows in one lobe and less in the other, but the total current does not change. Also, this saturated total current is proportional to the toroidal flux ${\it\Psi}$ .

10.2. Experimental results

The features described above are clearly demonstrated in the experiment. Figure 29 shows the peak plasma current as a function of the driving transformer voltage $V_{CB}$ , at fixed ${\it\Psi}$ . (Plasma current is a measure of ${\it\mu}a$ and $V_{CB}$ is a rough measure of $K/{\it\Psi}^{2}$ .) Note the following:

Figure 29. (a) Plasma current versus driving voltage and (b) variation of saturation current with toroidal flux in the multipinch. The value ${\it\theta}=1.56$ corresponds to ${\it\mu}a=2.42$ here. (Reproduced with permission from La Haye et al. (Reference La Haye, Jensen, Lee and Ohkawa1986).)

1. (i) At low $V_{CB}$ , the current is axisymmetric and symmetric about the equatorial mid-plane, i.e. equal in the upper and lower lobes, and it increases with $V_{CB}$ .

2. (ii) Above a critical $V_{CB}$ , the current remains axisymmetric but ceases to be symmetric about the mid-plane and it no longer increases with $V_{CB}$ . The total current saturates but more of it flows in one lobe and less in the other, and this asymmetry increases as $V_{CB}$ increases.

3. (iii) The saturation current, $I_{sat}$ , is proportional to the toroidal flux – in agreement with the relaxed-state prediction, $I_{sat}/{\it\Psi}\propto {\it\mu}a$ (see figure 29 b).

4. (iv) The observed proportionality between $I_{sat}$ and ${\it\Psi}$ corresponds to ${\it\mu}a\sim 2.42$ , in reasonable agreement with the theoretical value ${\it\mu}a=2.21$ for the zero-flux eigenvalue.

10.3. Flux generation

In fact, all the features of the experiment agree well with the properties calculated for the relaxed states – except one! According to theory, the total current should remain at its saturated value indefinitely as $V_{CB}$ increases beyond the critical value corresponding to ${\it\mu}=2.21$ . But this is not what happens in the experiment: at some point the current starts to increase again. The reason for this is that, as $V_{CB}$ increases, less and less current flows in one lobe and more in the other. This imbalance increases until at some point essentially all the current flows in one lobe and very little in the other. If current saturation were to continue beyond this point, as in theory it should, the current would have to reverse in one lobe. Instead, the discharge dies out in that lobe and becomes confined to the other, where it acts as a simple toroidal pinch whose current increases with $V_{CB}$ .

This feature of the multipinch brings out an important aspect of relaxation: it is not a passive decay process. It involves the generation of flux and current by plasma turbulence – as in a turbulent dynamo. When the multipinch reaches the state in which essentially all the current is in one lobe, and any further advance would require reverse current, the ‘dynamo’ is unable to act through the narrow gap from one lobe to the other – so the current remains confined to one lobe.

Finally we should note that current saturation in the multipinch is another example of spontaneous symmetry breaking, as it was in the cylinder. However, the symmetry involved is different in the two systems. In the cylinder it is the axial symmetry that is broken, whereas in the multipinch it is the symmetry about the equatorial plane.

11.1. Relaxed states of a spheromak

The helicity $K=\int \boldsymbol{A}\boldsymbol{\cdot }\boldsymbol{B}\;\text{d}{\it\tau}$ is conserved in a spheromak during relaxation (see § 7.1), and the relaxed state satisfies

Now, we recall that in a torus there were two invariants, $K$ and ${\it\Psi}$ , which were just sufficient to determine the relaxed state. However, in a spheromak, we have only one invariant, the helicity $K$ . The toroidal flux is no longer invariant because annihilation of flux can occur across the axis of symmetry. How, then, can the relaxed state be unique in a spheromak?

The key to this question is that in a torus there were two independent loop integrals on the surface of the torus: $\oint \boldsymbol{A}\boldsymbol{\cdot }\text{d}\boldsymbol{l}$ , the flux through the ‘hole’ in the torus; and $\oint \boldsymbol{A}\boldsymbol{\cdot }\text{d}\boldsymbol{s}={\it\Psi}$ , the toroidal flux inside the torus. But in a spheromak, being simply connected, all loop integrals on the surface are identical. In fact, since there is no flux through the boundary, they are all zero. Apart from an irrelevant gauge transformation this is equivalent to $\boldsymbol{A}\equiv 0$ on the boundary of a spheromak. But this is just the spheromak eigenvalue condition!

Hence for a spheromak, the only relaxed state is the lowest eigenfunction: the eigenvalue determines ${\it\mu}$ and the field amplitude is determined by the sole invariant $K$ . This means that, apart from a scale factor, the relaxed state in a spheromak is unique and entirely determined by the shape of the container.

11.2. Spheromak eigenfunctions

Axisymmetric eigenfunctions in a spheromak are easily found. The field is again written as (10.1) and the eigenvalue problem reduces to

with ${\it\chi}=0$ on the boundary, as in (10.6).

In a sphere of radius $a$ , the eigenfunction, in spherical coordinates ( ${\it\rho},{\it\theta},{\it\phi}$ ), is

where $j_{1}(x)\equiv \text{J}_{3/2}(x)/x^{1/2}$ , and the eigenvalue is ${\it\mu}a=4.49$ .

In a cylindrical container of height $h$ and radius $a$ (which, of course, is topologically spherical), the eigenfunction, in cylindrical coordinates $(r,{\it\theta},z)$ , is

where $kh={\rm\pi}$ , $la=3.83$ and the eigenvalue is ${\it\mu}^{2}=(3.83/a)^{2}+({\rm\pi}/h)^{2}$ .

11.3. Experimental data

In figure 32 the poloidal and toroidal field profiles given by (11.8) are compared with those measured in the BETA II experiment (Turner et al. Reference Turner, Goldenbaum, Granneman, Hammer, Hartman, Prono and Taska1983). This has a cylindrical flux conserver with $h/a\sim 1$ and a plasma created by a coaxial gun (see figure 31 a). The agreement is reasonable, particularly in view of the complex way the plasma is created, and it is particularly significant that the measured profiles remain essentially unchanged as the magnetic energy decays to ${\sim}1/8$ of its initial value. (Recall that the relaxed state of a spheromak is determined solely by the shape of the container.)

Figure 32. Magnetic field profiles in the BETA II spheromak. (a) Poloidal field and (b) toroidal field. Experimental measurements are shown as squares; the solid line is the theoretical prediction. (Reproduced with permission from Turner et al. (Reference Turner, Goldenbaum, Granneman, Hammer, Hartman, Prono and Taska1983). Copyright 1983, AIP Publishing LLC.)

The near-uniformity of ${\it\mu}$ in a relaxed state is supported by the measurements shown in figure 33, taken on the S-1 spheromak, which has an ellipsoidal flux conserver and a plasma produced by an inductive flux core (Hart et al. Reference Hart, Janos, Meyerhofer and Yamada1986). Figure 33(a) shows the ratio of poloidal current to poloidal flux across the mid-plane. The observations lie on a straight line whose slope agrees remarkably well with the computed value of ${\it\mu}\sim 5.5$ . Figure 33(b) shows how the ratio $J/B={\it\mu}$ on the mid-plane changes from an initial highly non-uniform profile (dashed line) to a much more uniform profile (solid line) after relaxation.

Figure 33. (a) Poloidal current versus poloidal flux and (b) measured profile of ${\it\mu}(r)$ in S-1 spheromak. (Reproduced with permission from Hart et al. (Reference Hart, Janos, Meyerhofer and Yamada1986). Copyright 1986, AIP Publishing LLC.)

Finally, figure 34 provides the most striking evidence for relaxation. This shows the evolution of the poloidal and toroidal fluxes, and of $q$ on the magnetic axis, during a plasma discharge in the S-1 spheromak (Janos et al. Reference Janos, Hart and Yamada1985a ,Reference Janos, Hart, Nam and Yamada b ). Figure 34(c) shows that $q$ rises rapidly from a small initial value and becomes close to that predicted for the relaxed state. It is then essentially constant for the remainder of the experiment. This rapid rise in $q$ is accompanied by an increase in toroidal flux (see figure 34 b) and a decrease in poloidal flux (see figure 34 a), during a period when there is strong plasma activity. This again shows that relaxation is not simply a passive decay but involves the creation and destruction of magnetic flux by plasma turbulence – or in this case the conversion of poloidal flux into toroidal flux.

Figure 34. Time dependence of magnetic fields in S-1 spheromak – see text for details. (Reproduced with permission from Janos et al. (Reference Janos, Hart and Yamada1985a ).)

12.1. Relative helicity

In the flux core spheromak, the plasma boundary is not a completely closed flux surface, and we noted earlier that in this case our definition of helicity is not gauge-invariant. This is more of a conceptual problem than a practical one, and can be overcome by introducing a modified definition of helicity.

The most intuitive modification is to imagine that the flux leaving or entering the boundary is extended outside as a vacuum field, where $\boldsymbol{{\rm\nabla}}\times \boldsymbol{B}=0$ . Then the integral $\int \boldsymbol{A}\boldsymbol{\cdot }\boldsymbol{B}\;\text{d}{\it\tau}$ , inside plus outside the sphere, is gauge-invariant. If the boundary of the sphere is a perfect conductor, the normal component of $\boldsymbol{B}$ is ‘frozen in’ and changes inside the sphere do not affect the hypothetical field outside. Therefore, the difference in $\int \boldsymbol{A}\boldsymbol{\cdot }\boldsymbol{B}\;\text{d}{\it\tau}$ for two configurations that differ inside the sphere but have the same $\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{B}$ on the surface (and therefore the same hypothetical extension field) is also well-defined and gauge-invariant. This is sometimes called the relative helicity, $K_{R}$ .

Note that $K_{R}$ must include the integral both inside and outside the sphere in order to be gauge-invariant, even though physically the latter does not change. This requirement reflects the fact that helicity is not a local quantity: it can be transferred from inside to outside by a gauge transformation!

The relative helicity satisfies

The exterior field does not change, therefore $\boldsymbol{E}_{exterior}=-\boldsymbol{{\rm\nabla}}{\it\phi}$ and

In the interior $\boldsymbol{E}+\boldsymbol{v}\times \boldsymbol{B}=0$ and the first integral vanishes. If the potential ${\it\phi}$ is constant over the surface, the last integral is also zero.

The purpose of introducing the flux core spheromak here is that it demonstrates another method of injecting helicity, in this case by a static voltage. If, instead of a uniform potential on the surface, one polar cap is held at a different potential from the other, helicity will be injected at a rate

where $V_{p}$ is the voltage between the polar caps and ${\it\Psi}_{p}$ is the flux through them (Taylor Reference Taylor1986).

12.2. Relaxed states of flux core spheromak

The relaxed state of the flux core spheromak satisfies the usual equation (11.1) and the relative helicity $K_{R}$ is invariant, but the appropriate solution is selected in a different way. In fact, a flux core spheromak can be operated in several ways. If the plasma relaxes from an initial state with given helicity $K_{R}$ , and with fixed flux ${\it\Psi}_{p}$ through the polar caps, then ${\it\mu}$ is specified by $K_{R}/{\it\Psi}_{p}^{2}$ somewhat in the same way as in a toroidal pinch (see figure 36 and § 9). However, this is only correct if the polar caps can provide whatever current is required by the resulting relaxed state. If this is not so, a voltage sheath would appear at the polar caps and the resulting voltage drop would reduce the helicity according to (12.3). This leads to another way of operating a flux core spheromak. If one polar cap is insulated from the other and connected to an external circuit that maintains a fixed current $I_{p}$ , then in the relaxed state ${\it\mu}=I_{p}/{\it\Psi}_{p}$ . This provides a method for external control of a flux core spheromak (see figure 37).

Figure 36. Analogy between tokamak and flux core spheromak.

Figure 37. Set-up to control relaxed state of flux core spheromak externally.

The relaxed states of a flux core spheromak are calculated in a similar way to those in a simple spheromak. Only the boundary conditions are different. Instead of $\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{n}=0$ on the boundary, it now has a fixed value on the polar caps and is zero elsewhere. A typical computed field profile is shown in figure 38. Analytic solutions for an idealised flux core system in which the polar caps are shrunk to a point while retaining finite flux through them are shown in figure 39. These show the changes that occur as the current through the polar caps, and therefore ${\it\mu}$ , is increased. When ${\it\mu}$ is much less than the lowest spheromak eigenvalue ( ${\it\mu}a=4.49$ ), the externally linked flux and current form a large part of the total flux and current. But as ${\it\mu}$ approaches the eigenvalue, the ratio of self-generated flux to externally linked flux increases indefinitely, and the externally linked flux becomes confined to a slim pencil along the axis of symmetry. This situation would constitute an extreme example of the generation of fields by the relaxation process. At the critical point all the fields would, in theory, be generated by the plasma itself.

Figure 38. Calculated field in a flux core spheromak, ${\it\mu}a=4.09$ (from Taylor Reference Taylor1986).

Figure 39. Evolution of analytic field profiles in an idealised flux core spheromak at various values of ${\it\mu}a$ : (a) ${\it\mu}a=0.001$ , (b) ${\it\mu}a=2.25$ , (c) ${\it\mu}a=2.70$ , (d) ${\it\mu}a=3.00$ , (e)  ${\it\mu}a=4.00$ , ( f) ${\it\mu}a=4.14$ and (g)  ${\it\mu}a=5.00$ . (Reproduced with permission from Turner (Reference Turner1984). Copyright 1984, AIP Publishing LLC.)

12.3. Experimental control of relaxed states

In a simple spheromak the field in the relaxed state will normally decay at the resistive diffusion rate (while retaining its profile). However, as we noted above, if the polar flux ${\it\Psi}_{p}$ and the polar current $I_{p}$ in a flux core spheromak are maintained, it should be possible to maintain its relaxed state at any amplitude and any ${\it\mu}=K_{R}/{\it\Psi}_{p}^{2}$ below the lowest eigenvalue. This idea is the basis for the experiments described below.

A somewhat distorted form of flux core spheromak, as shown in figure 40, has been extensively investigated at Los Alamos (CTX) (Jarboe et al. Reference Jarboe, Barnes, Platts and Wright1985; Barnes et al. Reference Barnes, Fernández, Henins, Hoida, Jarboe, Knox, Marklin and McKenna1986) and at UMIST (SPHEX) (Kitson & Browning Reference Kitson and Browning1990). The plasma in CTX is created by a plasma gun (on the left of the diagram). This injects plasma into a container, or ‘flux conserver’, at the right of the figure. If we concentrate attention on the flux conserver, we can see that there is a core of flux (shaded), which passes from the inner electrode of the gun, along the axis of the plasma in the flux conserver and then round the outside of the plasma back to the outer electrode of the gun. From the plasma viewpoint, the gun voltage appears across the flux core and acts as if it were between polar caps. This voltage can sustain the configuration against resistive decay for many times the usual resistive diffusion time.

Figure 40. Sustained configuration in CTX experiment (from Jarboe et al. Reference Jarboe, Barnes, Platts and Wright1985).

At this point we should note the time scales involved in this discussion:

1. (i) Relaxation time. This is much faster than the resistive diffusion time $a^{2}/{\it\eta}$ , and often close to the Alfvén time. This leads to the relaxed state. The helicity $K_{R}$ is invariant on this time scale.

2. (ii) The resistive diffusion time ${\sim}a^{2}/{\it\eta}$ . This is the time scale on which $K_{R}$ and the relaxed state itself decays, if it is not maintained by some external means.

3. (iii) The sustainment time. This is the time, longer than $a^{2}/{\it\eta}$ , on which the discharge can be maintained in its relaxed state by an external voltage.

12.4. The kinked Z-pinch

Once it was realised that a flux core spheromak could be maintained in the ways described above, it was found that they could also be created and maintained in more exotic configurations. One of these was the ‘kinked Z-pinch’ helicity source (Jarboe et al. Reference Jarboe, Barnes, Platts and Wright1985; Fernández et al. Reference Fernández, Wright, Marklin, Platts and Jarboe1989), shown in figure 41. In a highly idealised way, we can regard this device as creating a sequence of relaxed states in each of the three regions marked 1, 2 and 3. In region 1, we have a simple linear pinch maintained by the electrodes. Insofar as it is relaxed, ${\it\mu}$ would be set by the ratio of the current $I_{s}$ and flux ${\it\Psi}_{s}$ at the electrodes, ${\it\mu}_{1}=I_{s}/{\it\Psi}_{s}$ . In region 2, we have a helical relaxed state of a circular cylinder, with no net current or flux. If the cylinder were long, this would be the $m=1$ , $ka=1.23$ eigenfunction, which we discussed in connection with the cylindrical pinch. This would have ${\it\mu}_{2}a_{2}\sim 3.11$ as calculated in § 9.1. Finally, in region 3, we would have a spheromak with ${\it\mu}_{3}a_{3}$ , defined solely by the size and shape of the flux conserver (see § 11.2).

Figure 41. Illustration of kinked Z-pinch source.

Now, recall from § 9.1 that the energy for a given helicity in a relaxed state increases with ${\it\mu}$ . Consequently, it will be energetically favourable for the discharge to pass from region 1 to 2 and then into 3 if, and only if,

In fact, ${\it\mu}_{2}>{\it\mu}_{3}$ is set by the dimensions of the apparatus, but ${\it\mu}_{1}$ can be controlled externally and we can expect there to be a threshold for the successful injection of the spheromak into the flux conserver, set by

Such a threshold does seem to occur, as can be seen in figure 42, which shows a ‘figure of merit’ for the spheromak created in the flux conserver. This is roughly proportional to its total helicity. Figure 42(a) is for a 20 cm radius entrance cylinder, whilst figure 42(b) is for a 17 cm radius entrance cylinder. It is clear from these figures that there is no significant helicity input into the flux conserver unless ${\it\mu}_{1}>13~\text{m}^{-1}$ in (a), and ${\it\mu}_{1}>14~\text{m}^{-1}$ in (b). The theoretical threshold for a long cylinder would be ${\sim}15.6~\text{m}^{-1}$ for a 20 cm entrance radius.

Figure 42. Spheromak figure of merit ${\it\kappa}$ versus ${\it\lambda}_{e}\equiv {\it\mu}$ . (Reproduced with permission from Fernández et al. (Reference Fernández, Wright, Marklin, Platts and Jarboe1989). Copyright 1989, AIP Publishing LLC.)

13.1. Spheromak injection into tokamak

A remarkable application of helicity injection is the attempt (Brown & Bellan Reference Brown and Bellan1990) to influence the current in a tokamak by injecting helicity in the form of small spheromaks from a coaxial gun in a side arm into the main torus! This is illustrated in figure 43. As in the previous discussions, a small spheromak is created in the mouth of the gun and moves into the main tokamak chamber. To understand the observations, one must recall that helicity is a pseudo-scalar and therefore can be right-handed or left-handed. The handedness is related to the sign of ${\it\mu}$ , that is, to the relative directions of $\boldsymbol{J}$ and $\boldsymbol{B}$ .

Figure 43. Schematic of spheromak injection into ENCORE tokamak at Caltech. (Reproduced with permission from Brown & Bellan (Reference Brown and Bellan1990). Copyright 1990, AIP Publishing LLC.)

One can set up the main tokamak discharge with either right- or left-handed helicity by reversing the direction of the toroidal magnetic field, and one can create spheromaks with either right- or left-handed helicity by reversing the polarity of the gun voltage. Therefore, there are four permutations of the experiment:

1. (a) L-hand spheromak injected into L-hand tokamak,

2. (b) L-hand spheromak injected into R-hand tokamak,

3. (c) R-hand spheromak injected into R-hand tokamak,

4. (d) R-hand spheromak injected into L-hand tokamak.

The effect on the tokamak current of injecting the spheromak for the four cases is shown in figure 44. The dotted line in (a) shows the behaviour when there is no spheromak injection. It can be seen that when the injected spheromak and the tokamak have the same chirality (cases (a) and (c)), there is a small spike in the tokamak current. Conversely, when the spheromak and tokamak have opposite helicity (cases (b) and (d)), there is a small drop in the tokamak current. In all four cases the initial rise or drop in current is followed by a rapid decay. This is thought to be due to the cooling of the main plasma by the large amount of cold gas that accompanies the spheromak from the coaxial gun.

Figure 44. Effect on tokamak current of spheromak injection – see text for details. (Reproduced with permission from Brown & Bellan (Reference Brown and Bellan1990). Copyright 1990, AIP Publishing LLC.)

13.2. Wave injection

Another, untried, method of helicity injection is by plasma waves (Mett & Tataronis Reference Mett and Tataronis1989; Taylor Reference Taylor1989; Mett & Taylor Reference Mett and Taylor1992). One might envisage doing this by means of an antenna similar to those used for plasma heating.

Helicity injection by waves is hard to visualise (but see Chan et al. Reference Chan, Miller and Ohkawa1990). We have previously emphasised that helicity is a non-local quantity – the linking of lines of force – and it is difficult to relate this to wave motion. In fact, this process is more usually described as ‘current drive’, but it is interesting to see the connection with helicity. There are, of course, many types of plasma wave, but we will consider only the simplest – waves in a uniform plasma immersed in a uniform magnetic field $\boldsymbol{B}=(0,0,\boldsymbol{B}_{0})$ . However, this already illustrates some of the conceptual problems.

We start from what looks like a continuity equation for $\boldsymbol{A}\boldsymbol{\cdot }\boldsymbol{B}$ ,

where $\boldsymbol{Q}={\it\phi}\boldsymbol{B}+\boldsymbol{E}\times \boldsymbol{A}$ . For a plasma with an Ohm’s law $\boldsymbol{E}+\boldsymbol{v}\times \boldsymbol{B}={\it\eta}\,\boldsymbol{j}$ ,

Neither $\boldsymbol{A}\boldsymbol{\cdot }\boldsymbol{B}$ nor $\boldsymbol{Q}$ is gauge-invariant, but it is clear that $\boldsymbol{Q}$ must represent a flux of helicity. If we consider a small-amplitude wave so that $\boldsymbol{B}=\boldsymbol{B}_{0}+\boldsymbol{b}\text{e}^{\text{i}{\it\omega}t}$ , $\boldsymbol{Q}$ becomes

and we can therefore introduce a flux

which has the same divergence as $\boldsymbol{q}$ but is manifestly gauge-invariant. This form also brings out the connection between helicity and polarisation: a plane-polarised wave carries no helicity. (Actually this is obvious, as a plane-polarised wave has complete mirror symmetry whereas helicity is antisymmetric under reflection.)

However, a circularly polarised wave, propagating parallel to $\boldsymbol{B}_{0}$ , with a vector potential given by the real part of

does carry helicity in the form (13.7). (Here ${\it\omega}$ is real and $k$ complex, $k=k_{0}({\it\omega})+\text{i}k_{1}({\it\omega})$ , as for waves generated by an antenna at $z=0$ .)

From (13.8)

therefore

Define $|A|^{2}\text{e}^{-2k_{1}z}=A^{2}(z)$ , the square of the local wave amplitude, then

where $v={\it\omega}/k_{0}$ is the wave velocity. This can be interpreted as the wave carrying helicity density $k_{0}A^{2}$ . Finally, then,

Now consider a specific plasma wave, an Alfvén wave in a resistive plasma, with viscosity ${\it\nu}$ . The linearised wave equations are

This has the dispersion relation

where $v_{A}$ is the Alfvén speed (see (4.5)). When ${\it\eta}$ and ${\it\nu}$ are small, $k_{0}\approx {\it\omega}/v_{A}$ and the wave damping is $k_{1}\approx ({\it\eta}{\it\mu}_{0}^{-1}+{\it\nu})k_{0}^{2}/2v_{A}$ . Therefore

Thus it appears that we can indeed inject helicity from polarised waves, and that the input is proportional to the total wave damping, $({\it\eta}{\it\mu}_{0}^{-1}+{\it\nu})$ .

As mentioned above, this form of helicity injection is conventionally interpreted as current drive. The mean current induced by the wave fluctuations can be calculated directly from Ohm’s law ${\it\eta}\langle \,j_{z}\rangle =\langle \boldsymbol{v}\times \boldsymbol{b}\rangle$ . If the wave equation (13.14) is used to express $\boldsymbol{v}$ in terms of $\boldsymbol{b}$ , then

and using $k_{1}$ as given above, the wave-driven current is

This is in agreement with the net helicity injection given by (13.4) in steady state,

where $2{\it\eta}\langle \,\boldsymbol{j}_{{\it\omega}}\boldsymbol{\cdot }\boldsymbol{b}_{{\it\omega}}\rangle$ is the decay of helicity associated with the wave fluctuations (13.8) themselves. Note the remarkable fact that, if resistivity were the only damping, the loss of helicity through the wave fluctuations would be exactly twice the input from the wave.

16.1. Digression on adiabatic invariants

Adiabatic invariants are an interesting feature of classical mechanics and were important in the ‘old’ quantum theory – where they were the classical variables that became ‘quantised’, i.e. restricted to discrete values.

An adiabatic invariant arises when the motion of a classical particle, with Hamiltonian $H(p,q,a_{i})$ , is periodic – as it must be in one dimension if it is bounded (Goldstein Reference Goldstein1957). Then if the parameters $a_{i}$ are changed slowly compared to the period of oscillation, the resulting quasi-periodic motion possesses an invariant

which acts as a ‘constant of motion’ in a sense that will be explained shortly.

The classic example of such an invariant is associated with the motion of a pendulum whose length is changed slowly. In this case the invariant is the energy $E(t)$ divided by the frequency ${\it\omega}(t)$ of the pendulum,

i.e. as the length of the pendulum is altered, the change in its energy is proportional to the change in its frequency. (This is a hint that there might be a connection with quantum states, $\{E=I{\it\omega}\}_{classical}\rightarrow \{E=(n+1/2)\hbar{\it\omega}\}_{quantum}$ .)

It is important to note that the distinguishing feature of an invariant is not, as is sometimes stated, that it varies slowly – it is that changes in the invariant remain small for an indefinite time.

Now let us look at what we mean when we say that the invariant acts as a constant of motion.

The Hamiltonian for a harmonic oscillator with unit mass is

and, when ${\it\omega}$ is constant,

To find the invariant when ${\it\omega}(t)$ varies slowly, we first transform to new variables $P$ and $Q$ , such that

(The generating function for this transformation is $F(q,Q)=({\it\omega}/2)q^{2}\cot Q$ ; see Goldstein (Reference Goldstein1957).) Then the new Hamiltonian is

and the equations of motion become

We now introduce a slow time dependence of the oscillation frequency, via the small parameter ${\it\epsilon}$ , by writing ${\it\omega}={\it\omega}({\it\epsilon}t)$ , and look for an invariant as a power series in ${\it\epsilon}$ :

Then requiring

and using (16.11), we obtain, order by order, the recurrence relation

This allows us to generate the invariant term by term. The first term $J_{0}$ must be a function of $P$ only, and if we choose $J_{0}\equiv P$ then, using (16.8), we can confirm that this is indeed the integral $I$ introduced earlier (equation (16.3)),

The next few terms in the series are

and this pattern is followed by subsequent terms: odd-numbered terms ${\sim}\sin 2Q$ and even-numbered terms ${\sim}\cos 2Q$ plus a quantity independent of $Q$ .

So, now we see that $I=P=\oint p\,\text{d}q$ is actually only the first term of an expansion of the invariant $J$ in powers of ${\it\epsilon}$ . However, the important point is that all the higher-order terms in this series are small and vanish when ${\it\omega}(t)$ is constant. So if ${\it\omega}(t)$ starts from a steady value at $t\rightarrow -\infty$ , then changes slowly before eventually returning to a steady value at $t\rightarrow +\infty$ , the initial and final values of $I=\oint p\,\text{d}q$ should be the same. Hence, in this sense the series (16.12) is an exact constant of motion.

Perhaps surprisingly, this is not what actually happens. In fact, the initial and final values of $I$ differ by a small quantity that does not appear in the power series – because it cannot be expanded in ${\it\epsilon}$ . This difference, often called the ‘non-adiabatic change’, is of the form (Hastie et al. Reference Hastie, Hobbs and Taylor1969)

where $A$ and ${\it\alpha}$ are constants as ${\it\epsilon}\rightarrow 0$ and ${\it\phi}$ is the phase of the oscillation at some reference point (usually the point at which the phase is changing most slowly). This situation leads to the statement that an adiabatic invariant is ‘constant to all orders – but not a constant’.

(A hint that there must be a non-adiabatic contribution comes from the correspondence between a classical invariant and quantum states mentioned earlier. If the frequency of a quantum oscillator is altered slowly and then returned to its original value, the oscillator would not return to its original state. Instead, it would have an exponentially small probability of being in ‘excited’ states. The classical limit of this is the non-adiabatic change.)

A rigorous calculation of the non-adiabatic change (which, of course, cannot be obtained from an expansion or iteration in ${\it\epsilon}$ ) is very lengthy, but a strong indication of its exponential form can be given rather easily. If we return to (16.11),

and in the equation for ${\dot{P}}$ simply replace $Q$ by its first approximation,

we have

If ${\it\omega}({\it\tau})$ is analytic, this can be evaluated by contour integration, with the result

where $A$ and ${\it\alpha}$ are related to the order and position of the (complex) zero of ${\it\omega}({\it\tau})$ closest to the real axis.

16.2. Invariants and particle orbits

The motion of a charged particle in a time-dependent uniform field $B(t)$ is exactly equivalent to a harmonic oscillator. To see this, it is convenient to take the vector potential

Then the Hamiltonian is

from which we get the equations of motion in the form

where $z=x+\text{i}y$ . (The term in ${\dot{z}}$ is the Lorentz, $\boldsymbol{v}\times \boldsymbol{B}$ , term and that in ${\dot{B}}$ is due to the induced electric field.) Writing

reduces (16.27) to the oscillator equation

So we can carry over the results we obtained for the oscillator to the charged particle in a uniform field: then the lowest-order invariant is the magnetic moment $mv_{\bot }^{2}/B$ as we anticipated.

Unfortunately, we are interested in the motion of a particle in a static, non-uniform field, rather than in a time-dependent, uniform one. This is a far more complex problem that cannot be reduced to harmonic motion. Nevertheless, the qualitative features of the adiabatic invariant are similar to those for the harmonic oscillator. Thus there is a power-series expansion for the invariant, now in powers of the ratio of Larmor radius to the scale of the variation in magnetic field. Since this is proportional to $m/e$ , we can formally use ${\it\lambda}\equiv m/e$ as an expansion parameter; then the first term is $mv_{\bot }^{2}/2B$ and there is a non-adiabatic contribution of the form $A\exp (-{\it\alpha}/{\it\epsilon})\cos {\it\phi}$ .

A convenient starting point for calculating the power series is the Vlasov equation for a distribution of particles (see § 1.3). (Its connection with the invariants is, of course, that if $I_{j}$ are invariants, then $f(I_{j})$ is a stationary distribution.)

Introducing the gyro-phase angle ${\it\phi}$ about the magnetic field, the energy ${\it\epsilon}$ and the magnetic moment ${\it\mu}$ as velocity coordinates, and $({\it\alpha},{\it\beta},s)$ as spatial coordinates ( ${\it\alpha}$ and ${\it\beta}$ label a field line as in § 2.1 and $s$ is the distance along it), the Vlasov equation

can be written as

Here $D$ is a partial differential operator in $({\it\phi},{\it\mu},{\it\alpha},{\it\beta},s)$ , with coefficients depending on the gradients, curvature and torsion of the field lines.

We now seek a solution of (16.31),

that is stationary on ever increasing time scales, $t,t/{\it\lambda},t/{\it\lambda}^{2},\dots$  . It turns out (as detailed in Hastie et al. Reference Hastie, Taylor and Haas1967) that a distribution is stationary on all scales if it depends only on three variables – the energy ${\it\epsilon}$ , the magnetic moment ${\it\mu}_{0}=mv_{\bot }^{2}/2B$ , and a longitudinal invariant

The latter is associated with periodic motion along the field (between magnetic mirrors for trapped particles or around closed field lines for passing particles in a torus).

If we now write the true invariants $\hat{{\it\mu}}$ and ${\hat{J}}$ as

and identify the terms in the series calculated for a stationary distribution with the corresponding terms of the expansion

we obtain a recursion equation for the higher-order terms in the invariants $\hat{{\it\mu}}$ and ${\hat{J}}$ :

Here $D$ is the Vlasov operator (16.31), $L$ and $H$ are iterates of it, and ${\it\sigma}=\pm 1$ distinguishes the two possible directions of $v_{\Vert }$ . The first few terms, ${\it\mu}_{1}$ , ${\it\mu}_{2}$ and $J_{1}$ , of the invariants are given in Hastie et al. (Reference Hastie, Taylor and Haas1967). An unexpected feature is that they depend on the direction ${\it\sigma}$ of the particle motion along the magnetic field.

It is not clear that a rigorous calculation of the non-adiabatic contribution to ${\it\mu}$ exists for particles in a spatially varying field. Some calculations of it essentially replace the motion in a spatially varying field $B(\boldsymbol{r})$ by motion in a fictitious time-dependent field $B^{\ast }(t)$ – taken to be the field at the guiding centre of the particle, i.e.  $\text{d}B^{\ast }/\text{d}t=v_{\Vert }\,\partial B/\partial s$ . This, of course, leads to the exponential form (16.24) dependent on a complex zero of $B^{\ast }(t)$ . (The uncertainty about such a calculation lies in the fact that we are trying to calculate an exponentially small quantity when the guiding centre orbit is only accurate to order ${\it\epsilon}$ , a larger quantity. Moreover, whereas a particle in a slowly time-varying field experiences only the slow variation in $B(t)$ , a particle in a spatially varying field also experiences a ‘fast’ variation around the Larmor orbit.)

16.3. Orbit calculations

The magnetic moment is clearly the most important invariant for plasma confinement in mirrors, and it is interesting therefore to see its behaviour in particle orbits computed in a non-uniform field. This is illustrated in figure 48(a–c), reproduced from Hastie et al. (Reference Hastie, Hobbs and Taylor1969). (Actually, these are for orbits in a quadrupole field, as shown in figure 47, whose strength is given by

so that ${\it\Phi}={\rm\pi}$ is the weak-field plane of symmetry, $2a$ is the separation between the current filaments and $q$ labels field lines. In such a field there are both trapped orbits, near the weak-field plane, and untrapped circulating orbits.)

Figure 47. Quadrupole field configuration (from Hastie Reference Hastie and Dendy1993).

Figure 48(a) shows the change in the magnetic moment invariant, calculated from the approximations ${\it\mu}_{0},~({\it\mu}_{0}+{\it\lambda}{\it\mu}_{1})$ and $({\it\mu}_{0}+{\it\lambda}{\it\mu}_{1}+{\it\lambda}^{2}{\it\mu}_{2})$ to the infinite power series (16.34). This orbit starts at the weak-field plane and makes two reflections and two transits of the weak-field plane. The lowest approximation ${\it\mu}_{0}=mv_{\bot }^{2}/2B$ shows significant oscillations between transits and larger ‘jumps’ as it crosses the weak-field plane. As higher-order terms are added to the power series, the oscillations disappear, but the ‘jumps’ remain.

Figure 48. (a) Evolution of successive approximations to the particle magnetic moment as it moves in a mirror field. (b) Exponential form of the non-adiabatic jump at each transit as a function of $a/r_{L}\propto 1/v$ . (c) Long-time behaviour of a non-adiabatic particle. Here $q=1.19$ , ${\it\Phi}={\rm\pi}$ and $r/a=0.075$ . (Reproduced with permission from Hastie et al. (Reference Hastie, Hobbs and Taylor1969).)

The exponential form of these non-adiabatic jumps is confirmed by figure 48(b). (More extensive calculations of these jumps can be found in Howard (Reference Howard1971).)

Figure 48(c) shows the non-adiabatic changes over many transits of the weak-field plane. Owing to the large number of Larmor cycles that occur between successive transits, the phase at each transit is effectively random. After several transits, the particle is scattered from a trapped orbit to a passing one and back again (the boundary between trapped and passing orbits is shown by the dashed line).

It should be clearly understood that the parameters for these calculations were chosen to illustrate the points of interest and are not representative of typical mirror machines.

Finally, we should note that there is actually another invariant, in addition to $\hat{{\it\mu}}$ and ${\hat{J}}$ , for a particle in a mirror machine, associated with the periodic orbit of the guiding centre as it drifts around the axis of symmetry. The invariant in this case is the magnetic flux enclosed by the drift orbit. However, it is seldom invoked, as the drift motion can be studied directly (see § 20).

16.4. Summary

Adiabatic invariants are subtle quantities, and in the case of charged particles in a non-uniform magnetic field they are complicated functions of the field gradient, curvature and torsion. However, in practice, provided the Larmor radius is everywhere much smaller than the local field gradient $R=(\boldsymbol{{\rm\nabla}}B/B)^{-1}$ , the magnetic moment is effectively constant – apart from the effect of collisions.