On heat conduction in an irregular magnetic field. Part 1

Abstract

Anisotropic heat conduction in a plasma embedded in a magnetic field with irregular, possibly chaotic, field lines is discussed. If the collisional mean free path exceeds the electron gyroradius, the heat conductivity is much larger along the field lines than across them, and this enhances the transport across a domain where good flux surfaces do not exist. Recognising that anisotropic heat conduction may be cast in a variational form, and by constructing increasingly sophisticated trial functions that are based on invariant and almost-invariant structures under the magnetic field-line flow, bounds are derived on this enhancement and on the temperature variation along the magnetic field. In this way, remarkably accurate approximations for the temperature can be rapidly constructed without solving the diffusion equation, even in the small perpendicular-diffusion limit when the solution for the temperature is dominated by the fractal structure the magnetic field lines.

1. Introduction

In plasmas where the collisional mean free path, $\lambda$, is much longer than the electron gyroradius, $\rho$, heat is more easily conducted along the magnetic field lines than across them. The ratio of the classical heat conductivities perpendicular and parallel to the magnetic field is of the order of

(1.1)\begin{equation} \epsilon = \frac{\kappa_\perp}{\kappa_\|} \sim \left( \frac{\rho}{\lambda} \right)^{2} \ll 1, \end{equation}

and is extremely small in plasmas of interest in fusion research and astrophysics. For this reason, the confining quality of a magnetic field is very sensitive to the question of whether the field lines trace out magnetic surfaces. If they do not, e.g. if the magnetic field lines are chaotic, the dominant heat conduction along the field can significantly enhance the energy losses from the plasma.

There is a large body of literature discussing the extent of this enhancement and the underlying physical mechanisms in chaotic magnetic fields, see, e.g., Rechester & Rosenbluth (1978), Kadomtsev & Pogutse (1979), Krommes, Oberman & Kleva (1983) and Myra et al. (1993). In these works, scalings for the effective cross-field diffusion coefficient are derived in terms of local diffusion coefficients and the statistical properties of the magnetic field lines, such as Kolmogorov or correlation lengths.

Such information is meaningful insofar as these properties are approximately constant over the region under consideration. However, chaotic magnetic fields often contain an astonishing amount of structure, such as remnant flux surfaces, magnetic islands and partial transport barriers, which impart a similar structure, or variability, on the transport. For instance, the temperature is almost constant across a magnetic island if $\kappa _\| \gg \kappa _\perp$ and the island is sufficiently wide (Fitzpatrick 1995), and the temperature gradient is relatively large in the vicinity of an intact flux surface or partial barrier (Hudson & Breslau 2008). In other words, if the magnetic field is not completely chaotic, then the transport exhibits strong spatial variation that is not accounted for in the traditional theory of transport in chaotic fields.

Figure 1 illustrates this point, showing a Poincaré plot of field lines in an incompletely chaotic magnetic field, upon which isotherms from a heat conduction equation [see (1.2) with $\epsilon = 10^{-8}$] are superimposed in the left half of the figure. The distance between the isotherms is highly variable, and it is clear that the net transport cannot be described by anything like a constant ‘effective’ heat conductivity, such as that of Rechester and Rosenbluth. The latter may be appropriate in the limit of a strongly chaotic field, which is however rarely realised in practice.

Figure 1. A Poincaré puncture plot of a partially chaotic magnetic field and, in the left panel, isotherms of a solution of the anisotropic heat conduction equation (1.2).

Our aim is to understand the relation between magnetic-field geometry and transport. Following Hudson & Breslau (2008), we consider the simplest problem of anistotropic heat conduction in a magnetic field,

(1.2)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{q} = 0, \end{equation}

where

(1.3)\begin{equation} \boldsymbol{q} ={-} (\kappa_\| \boldsymbol{\nabla}_\| T +\kappa_\perp \boldsymbol{\nabla}_\perp T) \end{equation}

denotes the heat flux. Here $\boldsymbol {\nabla }_\| T = \boldsymbol {bb} \boldsymbol {\cdot }\boldsymbol {\nabla } T$ and $\boldsymbol {\nabla }_\perp T = \boldsymbol {\nabla } T - \boldsymbol {\nabla }_\| T$, with $\boldsymbol {b} = \boldsymbol {B}/B$ denoting the unit vector along the magnetic field, ${\bf B}$. This elliptic partial differential equation requires boundary conditions, which we shall take to be

(1.4)\begin{gather} T(\boldsymbol{r}) = T_0,\quad \boldsymbol{r} \in S_0, \end{gather}

(1.5)\begin{gather}T(\boldsymbol{r}) = T_\infty,\quad \boldsymbol{r} \in S_\infty, \end{gather}

where $S_0$ and $S_\infty$ are surfaces bounding the domain $\varOmega$ in which the heat conduction equation is solved. These surfaces are, thus, isothermals, and our primary problem consists of determining the other isothermals and deriving bounds on the net heat flux across the domain, defined by

(1.6)\begin{equation} Q = \int_{S_\infty} \boldsymbol{q} \boldsymbol{\cdot} \, {{\rm d}}\boldsymbol{S}={-}\int_{S_0} \boldsymbol{q} \boldsymbol{\cdot} \, {{\rm d}}\boldsymbol{S}, \end{equation}

where ${{\rm d}}\boldsymbol {S} = \boldsymbol {n} \,{{\rm d}}S$ is the surface element multiplied by the unit normal vector $\boldsymbol {n}$, pointing outward.

We shall assume that one of the boundary surfaces, $S_0$ say, is a magnetic surface, so that $\boldsymbol {b} \boldsymbol {\cdot } \boldsymbol {n} = 0$ everywhere on $S_0$. This is, for instance, the case of actual interest in a tokamak or stellarator with a so-called ergodic divertor, in which case $S_0$ could represent one of the intact magnetic surfaces and $S_\infty$ a surface further out in the plasma edge.

In classical transport theory (Braginskii 1965; Helander & Sigmar 2002), the conductivities $\kappa _\|$ and $\kappa _\perp$ depend on the temperature and the magnetic-field strength. However, these dependencies do not qualitatively affect the nature of the problem, particularly not in the limit in which the ratio $\epsilon = \kappa _\perp / \kappa _\|$ is very small, which is our main concern. In the next few sections, we therefore neglect this complication and instead take $\epsilon$ to be constant. We also neglect components of the classical heat flux other than those included in (1.3). In principle, one should also include the diamagnetic heat flux, which can be important. Its divergence is, however, relatively small if the magnetic field strength is almost constant. In the electron channel, there is also frictional heating as well as a component of the heat flux that is proportional to the relative velocity between the electron and ion fluids, but these terms in the energy equation are unimportant if the plasma current is sufficiently small. In low-collisionality plasmas, energy transport caused by wave–particle resonances can be more important than that due to chaotic field lines, but such transport is also ignored here.

2. Variational principle

The problem (1.2) is equivalent to finding the global minimum of the functional

(2.1)\begin{equation} D[T_{\rm trial}] = \frac{1}{2}\int_\varOmega (\kappa_\| | \boldsymbol{\nabla}_\| T_{\rm trial} |^{2} + \kappa_\perp |\boldsymbol{\nabla}_\perp T_{\rm trial} |^{2} ) \, {{\rm d}}V, \end{equation}

over all differentiable trial functions $T_{\rm trial}(\boldsymbol {r})$ satisfying the boundary conditions. These trial functions need not be particularly smooth and $\boldsymbol {\nabla } T_{\rm trial}$ is, for instance, allowed to be discontinuous, but should be square integrable. This quadratic form measures the total dissipation, which is minimised by the temperature distribution, $T(\boldsymbol {r})$, satisfying the differential equation (1.2). To see this, it is helpful to write the heat flux in matrix form,

(2.2)\begin{equation} q_i ={-} \kappa_{ij} {\partial}_j T, \end{equation}

and note that the conductivity matrix is symmetric, $\kappa _{ij} = \kappa _{ji}$, so that the first-order variation of

(2.3)\begin{equation} D[T] = \frac{1}{2} \int_\varOmega \kappa_{ij}({\partial}_i T) ({\partial}_j T) \,{{\rm d}}V, \end{equation}

becomes

(2.4)\begin{equation} \delta D[T] = \int_\varOmega \kappa_{ij} ({\partial}_i T) ({\partial}_j \delta T) \,{{\rm d}}V ={-} \int_\varOmega {\partial}_j (\kappa_{ij} {\partial}_i T) \delta T \,{{\rm d}}V. \end{equation}

if $\delta T(\boldsymbol {r})$ vanishes on the boundary. It is now clear that the differential equation (1.2) is equivalent to the variational principle $\delta D[T] = 0$, where the temperature is to be held constant on the boundaries. The functional $D$ has no maximum, and it follows that any minimum is attained by a solution to (1.2). That the stationary point is indeed a minimum follows from the fact that the second variation $\delta ^{2} D[\delta T,\delta T]$ is positive for all $\delta T \ne 0$ independently of $T$.

This minimum is related to the net heat flux by

(2.5)\begin{equation} Q = \frac{2 D[T]}{T_0 - T_\infty}, \end{equation}

as follows from

(2.6)\begin{equation} 0 = \int_\varOmega T \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{q} \,{{\rm d}}V ={-} \int_\varOmega \boldsymbol{q} \boldsymbol{\cdot}\boldsymbol{\nabla} T \,{{\rm d}}V + \int_{S_0 \cup S_\infty} T \boldsymbol{q} \boldsymbol{\cdot} {{\rm d}}\boldsymbol{S} = 2 D[T] - Q (T_0 - T_\infty). \end{equation}

3. Net heat flux

The variational property of $D[T]$ and the relation (2.5) are useful for deriving upper and lower bounds on the net transport. A convenient dimensionless measure of the latter is

(3.1)\begin{equation} \hat Q = \frac{Q}{\kappa_\perp (T_0 - T_\infty) A_0/L}, \end{equation}

Here $L$ denotes an average distance between $S_0$ and $S_\infty$, which we define as

(3.2)\begin{equation} L = \frac{2V}{A_0 + A_\infty}, \end{equation}

where $A_0$ and $A_\infty$ denote the areas of the boundaries and $V$ the volume of the region between them.

3.1. Lower bound

A lower bound is easily obtained by noting that the solution $T$ to (1.2) satisfies

(3.3)\begin{equation} 2 D[T] \ge \kappa_\perp \int_\varOmega |\boldsymbol{\nabla} T|^{2} \,{{\rm d}}V, \end{equation}

where the integral on the right-hand side necessarily exceeds its minimum taken over all functions satisfying the boundary conditions,

(3.4)\begin{equation} \int_\varOmega |\boldsymbol{\nabla} T|^{2} \,{{\rm d}}V \ge \inf_{T(\boldsymbol{r})} \int_\varOmega |\boldsymbol{\nabla} T|^{2} \,{{\rm d}}V. \end{equation}

Equality is realised by the (unique) harmonic function, $\nabla ^{2} T = 0$, that equals $T_0$ on $S_0$ and $T_\infty$ on $S_\infty$. If we define a dimensionless constant $K$ characterising the geometry of the domain by

(3.5)\begin{equation} \inf_{T(\boldsymbol{r})} \int_\varOmega |\boldsymbol{\nabla} T|^{2} \,{{\rm d}}V = \frac{(T_0 - T_\infty)^{2} A_0}{L}\cdot K, \end{equation}

then it can be concluded from (2.5), (3.1) and (3.3) that the net transport is bounded from below by

(3.6)\begin{equation} \hat Q \ge K. \end{equation}

This result is not surprising: it merely says that the net transport will at least be that expected from isotropic heat conduction with the conductivity $\kappa _\perp$.

3.2. Upper bound

A more interesting question is to what extent the transport is enhanced by parallel heat conduction. An upper bound on this enhancement can be obtained from the fact that the total dissipation, (2.1), associated with the solution of (1.2) is smaller than, or equal to, that from any trial function $T_{\rm trial}$,

(3.7)\begin{equation} D[T] \le D[T_{\rm trial}]. \end{equation}

Of course, in order to make this upper bound as ‘tight’ as possible, the trial function should be chosen optimally. In the case of interest here, where the boundary surface $S_0$ is tangential to the magnetic field, the normal component of the field is small in the vicinity of this boundary. Indeed, if the magnetic field is continuously differentiable, then for some length $\ell _0 > 0$ the normal component of the magnetic field satisfies

(3.8)\begin{equation} |\boldsymbol{b} \boldsymbol{\cdot} \boldsymbol{n}| < \frac{z}{\ell_0}, \end{equation}

where $z$ denotes the distance from the boundary $S_0$. For the argument that follows, it is sufficient that this inequality be satisfied in a small neighbourhood of $S_0$.

A useful choice of trial function is (see figure 2)

(3.9)\begin{equation} T_{\rm trial}(z) = \left\{\begin{array}{@{}ll} T_0 - \dfrac{(T_0 - T_\infty) z}{\delta}, \quad & 0 \le z < \delta, \\ T_\infty, & z > \delta \end{array}\right. \end{equation}

so that

(3.10)\begin{gather} |\boldsymbol{\nabla}_\| T_{\rm trial}| < \frac{(T_0 - T_\infty) z}{\ell_0 \delta}, \end{gather}

(3.11)\begin{gather}|\boldsymbol{\nabla}_\perp T_{\rm trial}| \le \frac{T_0 - T_\infty}{\delta}, \end{gather}

and thus

(3.12)\begin{equation} D[T_{\rm trial}] \le \frac{1}{2} \int_{z<\delta} \left( \kappa_\| | \boldsymbol{\nabla}_\| T |^{2} + \kappa_\perp |\boldsymbol{\nabla}_\perp T |^{2} \right) \, {{\rm d}}V \end{equation}

is bounded by

(3.13)\begin{equation} D[T_{\rm trial}] < \frac{\kappa_\| A_0 (T_0 - T_\infty)^{2}}{2} \left( \frac{\delta}{3 \ell_0^{2}} + \frac{\epsilon}{\delta} \right), \end{equation}

where the volume of the boundary layer is $A_0 \delta$ plus a small term of order $\delta ^{2}$, which we ignore.

Figure 2. The trial function (3.9).

The optimal choice of $\delta = \ell _0 (3 \epsilon )^{1/2}$ now gives

(3.14)\begin{equation} D[T_{\rm trial}] < \frac{\kappa_\| A_0 (T_0 - T_\infty)^{2}}{\ell_0} \sqrt{\frac{\epsilon}{3}}.\end{equation}

Recalling (2.5), (3.1) and (3.6), we thus conclude

(3.15)\begin{equation} K \le \hat Q < \frac{L}{\ell_0 \sqrt{3\epsilon}}. \end{equation}

Within a factor of order unity that depends the geometry of the domain, the lower bound is realised if the magnetic field is entirely regular, and we hypothesise that the scaling with $\epsilon$ of the upper bound applies in a strongly chaotic field, i.e. we propose that the upper bound cannot be improved except by lowering the numerical coefficient.

4. Transport barriers

In the previous section, nothing was assumed about the magnetic field except that it should be tangential to one of the boundaries. The bounds (3.15) hold irrespective of whether the magnetic field in the interior of the domain $\varOmega$ traces out nested flux surfaces, forms magnetic islands or has chaotic field lines.

4.1. One flux surface

If some good flux surfaces within $\varOmega$ do exist, the upper bound can, however, be improved. To see how this can be done, suppose that one such surface, $S_1$, lies somewhere between $S_0$ and $S_\infty$. Let $A_1$ be the area of this surface, and introduce a coordinate $z$ such that $z=z_0=0$ on $S_0$, $z=z_1$ on $S_1$, and $|\boldsymbol {\nabla } z| = 1$ in the vicinity of both $S_0$ and $S_1$. The magnetic field then satisfies $\boldsymbol {b} \boldsymbol {\cdot }\boldsymbol {\nabla } z = 0$ on $S_0$ and $S_1$, and there are numbers $\ell _0$ and $\ell _1$ such that

(4.1)\begin{equation} |\boldsymbol{b} \boldsymbol{\cdot}\boldsymbol{\nabla} z| < \left| \frac{z-z_i}{\ell_i} \right| \quad \mbox{in the vicinity of } S_i \end{equation}

for $i \in \{0,1\}$. For the following trial function (figure 3)

(4.2)\begin{equation} T_{\rm trial}(z) = \left\{\begin{array}{@{}ll} T_0 - (T_0 - T_1) \dfrac{z}{\delta_0}, & 0 < z < \delta_0,\\ T_1, & \delta_0 < z < z_1 - \delta_1,\\ T_1 - (T_1 - T_\infty) \dfrac{z - z_1 + \delta_1}{2 \delta_1}, \quad & - \delta_1 < z - z_1 < \delta_1,\\ T_\infty, & \delta_1 < z - z_1, \end{array}\right. \end{equation}

the gradients satisfy the following inequalities

(4.3)\begin{gather} | \boldsymbol{\nabla}_\| T_{\rm trial} | < \left| \frac{(T_0 - T_1) z}{\delta_0 \ell_0}\right|, \end{gather}

(4.4)\begin{gather}|\boldsymbol{\nabla}_\perp T_{\rm trial} | < \left| \frac{T_0 - T_1}{\delta_0} \right|, \end{gather}

for $0 \le z \le \delta _0$, and analogous inequalities for $-\delta _1 \le z-z_1 \le \delta _1$,

(4.5)\begin{gather} | \boldsymbol{\nabla}_\| T_{\rm trial} | < \left| \frac{(T_1 - T_\infty) z}{2 \delta_1 \ell_1} \right|, \end{gather}

(4.6)\begin{gather}|\boldsymbol{\nabla}_\perp T_{\rm trial} | < \left| \frac{T_1 - T_\infty}{2\delta_1} \right|. \end{gather}

Figure 3. The trial function (4.2).

Our variational form (2.1) is thus bounded from above by

(4.7)\begin{equation} D[T_{\rm trial}] \le \frac{\kappa_\| A_0 (T_0 - T_1)^{2}}{2} \left( \frac{\delta_0}{3 \ell_0^{2}} + \frac{\epsilon}{\delta_0} \right) + \frac{\kappa_\| A_1 (T_1 - T_\infty)^{2}}{4} \left( \frac{\delta_1}{3 \ell_1^{2}} + \frac{\epsilon}{\delta_1} \right), \end{equation}

where $A_1$ denotes the area of $S_1$. This expression is minimised by choosing $\delta _j = \ell _j (3 \epsilon )^{1/2}$ and

(4.8)\begin{equation} \frac{T_1 - T_\infty}{T_0 - T_\infty} = \frac{1}{\dfrac{\ell_0 A_1 }{2 \ell_1 A_0} + 1}. \end{equation}

The result then becomes

(4.9)\begin{equation} D[T_{\rm trial}] < \frac{\kappa_\| A_0 (T_0 - T_\infty)^{2}}{\ell_0 \left(1 + \dfrac{2 \ell_1 A_0}{\ell_0 A_1}\right)} \sqrt{\frac{\epsilon}{3}}, \end{equation}

which is smaller than (3.14), in fact three times smaller if $\ell _0 = \ell _1$ and $A_0 = A_1$.

4.2. Several flux surfaces

It is straightforward to generalise this result to the case of $n>1$ intact flux surfaces, situated at $z = z_i$, having areas $A_i$, and in whose vicinity $|\boldsymbol {\nabla } z | = 1$ and

(4.10)\begin{equation} |\boldsymbol{b}\boldsymbol{\cdot}\boldsymbol{\nabla} z| < \left| \frac{z-z_i}{\ell_i} \right|. \end{equation}

for $1 \le i \le n$. An appropriate trial function is

(4.11)\begin{equation} T_{\rm trial}(z) = \left\{\begin{array}{@{}ll} T_0 - (T_0 - T_1) \dfrac{z}{\delta_0}, & 0 < z < \delta_0\\ T_i, & z_{i-1} + \delta_{i-1} < z < z_i - \delta_i\\ T_i - (T_i - T_{i+1}) \dfrac{z - z_i + \delta_i}{2 \delta_i}, & - \delta_i < z - z_i < \delta_i\\ T_\infty, & \delta_n < z - z_n, \end{array}\right.\end{equation}

where we set $T_{n+1} = T_\infty$ and $\delta _i = \ell _i \sqrt {3\epsilon }$. Note the qualitative similarity between this trial function and the numerical solution of the heat conduction equation obtained by Hudson & Breslau (2008). With this choice of $T_{\rm trial}$, we obtain

(4.12)\begin{equation} D[T_{\rm trial}] < \kappa_\| \sqrt{\frac{\epsilon}{3}} \left[\frac{A_0 (T_0 - T_1)^{2}}{\ell_0} + \sum_{i=1}^{n} \frac{A_i (T_i - T_{i+1})^{2}}{2 \ell_i} \right], \end{equation}

which is minimised by requiring

(4.13)\begin{equation} \frac{{\partial} D}{{\partial} T_i} = 0,\quad i=1,\ldots, n. \end{equation}

This equation gives

(4.14)\begin{equation} \frac{2 A_0}{\ell_0} (T_0-T_1) = \frac{A_i}{\ell_i} (T_i - T_{i+1}) = \frac{A_n}{\ell_n} (T_n - T_\infty), \end{equation}

for $1 \le i \le n-1$, which implies

(4.15)\begin{equation} T_0 - T_\infty = (T_n -T_\infty) \frac{A_n}{\ell_n} \left( \frac{\ell_0}{2 A_0} + \sum_{i=1}^{n} \frac{\ell_i}{A_i} \right), \end{equation}

and, finally, leads to

(4.16)\begin{equation} D[T_{\rm trial}] < \left. \kappa_\|\sqrt{\frac{\epsilon}{3}} \frac{(T_0-T_\infty)^{2}}{2} \middle/ \left( \frac{\ell_0}{2 A_0} + \sum_{i=1}^{n} \frac{\ell_i}{A_i} \right). \right. \end{equation}

In particular, if all $\ell _i$ and all $A_i$ are the same, then this result is equal to our earlier result (3.14) divided by $2n+1$. Thus, if there are many intact surfaces within the domain $\varOmega$, then the upper bound on the heat flux is much smaller than the basic estimate (3.15). If one of the surfaces has a large ratio $\ell _i / A_i$, then this surface effectively acts as a bottleneck and the overall transport becomes particularly small.

4.3. Imperfect flux surfaces

Fully intact magnetic flux surfaces thus impede the transport considerably, but this is true also for a wider class of surfaces known as partial barriers in the literature on dynamical systems (Meiss 1992). Magnetic field lines pass through partial barriers, unlike flux surfaces, but they do so slowly in some sense.¹ If $S_i$ is a surface with normal vector $\boldsymbol {n}$ and area $A_i$, we shall declare it to be an imperfect flux surface if the number $R_i$ defined by

(4.17)\begin{equation} R_i = \frac{1}{A_i} \int_{S_i} (\boldsymbol{b} \boldsymbol{\cdot} \boldsymbol{n})^{2} dS \end{equation}

is much smaller than unity. Note the similarity to the quadratic-flux-minimising (QFM) surfaces of Dewar, Hudson & Price (1994) and Hudson & Dewar (2009).² Moreover, we assume that this surface is optimally positioned in the sense that the quadratic flux through neighbouring surfaces is larger than $R_i A_i$ and satisfies

(4.18)\begin{equation} \int_{z = {\rm constant}} (\boldsymbol{b} \boldsymbol{\cdot} \boldsymbol{n})^{2} \,{{\rm d}}S \le A_i \left[ R_i + \left(\frac{z-z_i}{\lambda_i} \right)^{2} \right] \end{equation}

for some number $\lambda _i$ if $|z-z_i| < \delta _i$. If the trial function is chosen in the same way (4.11) as before, then

(4.19)\begin{equation} \frac{1}{2} \int_{|z - z_i| < \delta_i} \left( \boldsymbol{\nabla}_\| T |^{2} + \epsilon |\boldsymbol{\nabla}_\perp T |^{2} \right) \, {{\rm d}}V \le \frac{A_i (T_i - T_{i+1})^{2}}{4 \delta_i} \left( R_i + \epsilon + \frac{\delta_i^{2}}{3 \lambda_i^{2}} \right), \end{equation}

where the right-hand side is minimised by choosing

(4.20)\begin{equation} \left( \frac{\delta_i}{\lambda_i}\right)^{2} = 3(R_i + \epsilon) \end{equation}

and then becomes equal to

(4.21)\begin{equation} \frac{A_i (T_i - T_{i+1})^{2}}{2 \lambda_i} \sqrt{\frac{R_i + \epsilon}{3}}. \end{equation}

In other words, if we write

(4.22)\begin{equation} \ell_i = \frac{\lambda_i}{\sqrt{1 + \dfrac{R_i}{\epsilon}}}, \end{equation}

then we have

(4.23)\begin{equation} \frac{1}{2} \int_{|z - z_i| < \delta_i} \left( |\boldsymbol{\nabla}_\| T |^{2} + \epsilon |\boldsymbol{\nabla}_\perp T |^{2} \right) {{\rm d}}V = \sqrt{\frac{\epsilon}{3}} \frac{A_i (T_i - T_{i+1})^{2}}{2 \ell_i}, \end{equation}

a result strongly reminiscent of (4.12). Hence, we conclude that if the domain $\varOmega$ is permeated by a number of surfaces $S_i$, some of which are ordinary flux surfaces and the others imperfect flux surfaces, then the dissipation is bounded from above by (4.12) where the numbers $l_i$ are defined by (4.10) for the ordinary flux surfaces and by (4.22) for the imperfect ones. As is clear from this equation, an imperfect flux surface impedes the transport roughly much as a full flux surface does if the quadratic flux through it is small enough, $R_i \le \epsilon$.

4.4. Bounds on the net transport

All the results in this section can be summarised as follows. When ordinary and imperfect flux surfaces are present, then the normalised heat flux (3.1) is subject to the following bounds:

(4.24)\begin{equation} K \le \hat Q < \left. \frac{L}{\sqrt{3\epsilon}}\middle/ \left(\ell_0 + 2 A_0 \sum_{i=1}^{n} \frac{\ell_i}{A_i} \right), \right. \end{equation}

where the numbers $\ell _i$ are defined by (4.10) and (4.22).

5. A continuous set of surfaces

We may go one step further and construct the optimal upper bound on the net heat flux given a continuous set of isotherms. Suppose we are able to guess the approximate geometry of all the isotherms by inspecting the geometry of the magnetic field, i.e. we know that the temperature is approximately constant on surfaces of constant $z(\boldsymbol {r}) \in [0,1]$, some known function of space, but we do not know the function $T(z)$. It is then possible to construct the ‘optimal’ such function in the sense of the variational principle, by minimising $D[T]$ subject to the constraint that $T$ depend only on $z$. If we introduce two other coordinates $(x,y)$ and the Jacobian $\sqrt {g} = 1/(\boldsymbol {\nabla } x \times \boldsymbol {\nabla } y) \boldsymbol {\cdot } \boldsymbol {\nabla } z > 0$, then

(5.1)\begin{equation} D[T] =\frac{\kappa_\|}{2} \int \left( |\boldsymbol{\nabla}_\| T|^{2} + \epsilon |\boldsymbol{\nabla}_\perp T|^{2} \right) \, {{\rm d}}V = \frac{1}{2} \int_0^{1} k(z) \left(\frac{{{\rm d}}T}{{{\rm d}}z} \right)^{2} \, {{\rm d}}z, \end{equation}

where

(5.2)\begin{equation} k(z) = \kappa_\| \int_{z = const.} \left[(\boldsymbol{b} \boldsymbol{\cdot} \boldsymbol{n})^{2} + \epsilon |\boldsymbol{\nabla} z|^{2} \right] \sqrt{g} \,{{\rm d}x}\,{{\rm d}y}. \end{equation}

The Euler–Lagrange equation for the function $T(z)$ that minimises (5.1) is

(5.3)\begin{equation} \frac{{{\rm d}}}{{{\rm d}}z} \left( k \frac{{{\rm d}}T}{{{\rm d}}z} \right) = 0, \end{equation}

and implies (Hudson 2009)

(5.4)\begin{equation} \frac{{{\rm d}}T}{{{\rm d}}z} = \frac{C}{k(z)}, \end{equation}

where $C$ is an integration constant, which is determined by the boundary conditions,

(5.5)\begin{equation} C \int_0^{1} \frac{{{\rm d}}z}{k(z)} = T_\infty- T_0. \end{equation}

The optimal trial function is, thus, given by

(5.6)\begin{equation} T(z) = \left. T_0 + (T_\infty - T_0) \int_0^{z} \frac{{{\rm d}}z'}{k(z')}\middle/ \int_0^{1} \frac{{{\rm d}}z}{k(z)}. \right. \end{equation}

Substituting this result in (5.1) gives

(5.7)\begin{equation} \left. D[T] = \frac{(T_0 - T_\infty)^{2}}{2} \middle/ \int_0^{1} \frac{{{\rm d}}z}{k(z)}, \right. \end{equation}

and it follows from (2.5) that

(5.8)\begin{equation} \left. Q = (T_0 - T_\infty) \middle/ \int_0^{1} \frac{{{\rm d}}z}{k(z)}. \right. \end{equation}

From our variational principle, we thus conclude that the total dissipation is given by this formula if the temperature is constant on surfaces of constant $z$ and is otherwise bounded by it from above. Our dimensionless measure of the total heat flux (3.1) thus satisfies

(5.9)\begin{equation} \hat Q \le \frac{\kappa_\| L}{ \left. \kappa_\perp A \int_0^{1} {{\rm d}}z \middle/ \int_{z = const.} \right.\left[(\boldsymbol{b} \boldsymbol{\cdot} \boldsymbol{n})^{2} + \epsilon |\boldsymbol{\nabla} z|^{2} \right] \sqrt{g} \,{{\rm d}x}\,{{\rm d}y}}. \end{equation}

6. Parallel temperature variation

The result (3.14) can be used to derive an upper limit on the temperature variation along the magnetic field. The bound is most effective in subdomains $\varPhi$ of $\varOmega$ aligned with the magnetic field, usually known as ‘flux tubes’ in the plasma physics literature. A flux tube is a region whose boundary can be written as ${\partial } \varPhi = F_a \cup S \cup F_b$, where the magnetic field is tangential to the ‘side surface’ $S$ and perpendicular to the ‘end faces’ $F_a$ and $F_b$. (It is assumed that the magnetic field does not vanish anywhere.) Because of the divergence-free nature of magnetic fields, the latter are connected by magnetic field lines, i.e. from each point on $F_a$ there is a magnetic field line to a point on $F_b$, and vice versa.

Suppose that the temperature difference between the coldest point on $F_a$ and the hottest point on $F_b$ exceeds $\Delta T > 0$. The total dissipation in the flux tube $\varPhi$ is certainly less than that in the entire domain $\varOmega$,

(6.1)\begin{equation} D_\varPhi[T(\boldsymbol{r})]\le D[T(\boldsymbol{r})], \end{equation}

and can be bounded from below by neglecting perpendicular heat conduction

(6.2)\begin{equation} D_\varPhi[T(\boldsymbol{r})] \ge \frac{1}{2} \int_\varPhi \kappa_\| | \boldsymbol{\nabla}_\| T |^{2} \,{{\rm d}}V. \end{equation}

If $l$ denotes the arc length along the field and $dS$ the area element perpendicular to it,

(6.3)\begin{equation} D_\varPhi[T(\boldsymbol{r})] \ge \frac{\kappa_\|}{2} \int {{\rm d}}S \int \left( \frac{{\partial} T}{{\partial} l} \right)^{2} \, {{\rm d}}l, \end{equation}

where the last integral is bounded from below by

(6.4)\begin{equation} \int \left( \frac{{\partial} T}{{\partial} l} \right)^{2} {{\rm d}}l \ge \frac{(\Delta T)^{2}}{L_\|}, \end{equation}

where $L_\|$ is the length of the of the longest field line within the flux tube. If we denote the area of its smallest cross section by $a$, we thus have

(6.5)\begin{equation} D_\varPhi[T(\boldsymbol{r})] \ge \frac{\kappa_\| a (\Delta T)^{2}}{2 L_\|}, \end{equation}

which we combine with (3.14) and (6.1) to establish

(6.6)\begin{equation} \left| \frac{\Delta T}{T_0 - T_\infty} \right| < \left( \frac{L_\| A}{\ell_0 a} \right)^{1/2} \left( \frac{\epsilon}{3} \right)^{1/4}. \end{equation}

It follows from this inequality that the temperature variation along a flux tube approaches zero as $\epsilon \rightarrow 0$.

7. Numerical results

We consider a magnetic field given by the vector potential

(7.1)\begin{equation} \boldsymbol{A} = \rho \boldsymbol{\nabla} \theta - \chi(\rho,\theta,\zeta) \boldsymbol{\nabla} \zeta, \end{equation}

where $(\rho,\theta,\zeta )$ are ‘toroidal’ coordinates, with $\rho$ being the radial or action coordinate, $\theta$ being the poloidal angle or position coordinate and $\zeta$ being the toroidal angle or time coordinate; and where the field-line Hamiltonian is given by

(7.2)\begin{equation} \chi(\rho,\theta,\zeta) = \chi_0(\rho) + \sum \chi_{m,n}(\rho) \cos(m\theta-n\zeta), \end{equation}

where $\chi _0(\rho )$ and $\chi _{m,n}(\rho )$ are commonly referred to as the integrable Hamiltonian and the perturbation.

Any toroidal magnetic field that is nearly integrable can be written in such a form. A transformation to straight field-line coordinates, $\boldsymbol {x}(\rho,\theta,\zeta )$ for which $\chi = \chi (\rho )$, for the integrable part of the nearly integrable magnetic field might be first required. A construction of a nearby integrable field for a given non-integrable field is given by Hudson & Dewar (1998).

The reconnected flux inside magnetic islands, the existence of flux surfaces, the flux through cantori, and the fractal structure of phase space are completely determined by the field-line Hamiltonian. The ‘average’ or ‘background’ rotational transform profile is given by

(7.3)\begin{equation} \text{$\iota\!\!$-}_0(\rho) \equiv \frac{\partial \chi_0(\rho)}{\partial \rho}. \end{equation}

The actual rotational transform profile will be influenced by the perturbations. For irregular ‘chaotic’ field lines, the rotational transform may not be defined because $\Delta \theta /\Delta \zeta$, where $\Delta \theta$ is the increment in the poloidal angle along a field line, may not converge even as $\Delta \zeta \rightarrow \infty$. The geometry of flux surfaces etc. is determined by the coordinate transformation, $\boldsymbol {x}(\rho,\theta,\zeta )$.

In this article, we are primarily concerned with understanding how the non-integrabilty of the magnetic field influences transport. For the geometry, for convenience, we consider a doubly periodic Cartesian slab, $\boldsymbol {x}(\rho,\theta,\zeta ) = \theta \boldsymbol {i} + \zeta \boldsymbol {j} + \rho \boldsymbol {k}$, with $\rho \in [0,1]$. We set the rotational transform profile to vary between $\text{$\iota\!\!$-}{(1)}=0$ and $\text{$\iota\!\!$-}{(1)}=1$ by choosing $\chi _0(\rho ) = \rho ^{2} / 2$.

For each $\chi _{m,n}(\rho ) \ne 0$, there will be an $\text{$\iota\!\!$-}=n/m$ island. For each pair of $n/m$ and $n'/m'$ islands, there will be a smaller $(n+n')/(m+m')$ island, and so on ad infinitum. If the islands are sufficiently large, then large regions of extended irregular trajectories (chaos) will emerge.

Throughout this section, we choose our domain to be bounded by $\rho =0$ and $\rho =1$ and impose boundary conditions $T(0,\theta,\zeta )=0$ and $T(1,\theta,\zeta )=1$.

7.1. Transport with one intact flux surface

In this section, we verify the lower (§ 3.1) and upper (§ 4.1) bounds on the transport in the presence of one intact flux surface. We will assume both of the bounding surfaces, $S_0$ and $S_{\infty }$, to be magnetic surfaces and the inequality (4.10) to be satisfied in the neighbourhood of each surface. In this case, the upper bound on the transport can be further reduced by assuming a trial function of the form,

(7.4)\begin{equation} T_{\rm trial}(z) = \left\{\begin{array}{@{}ll} T_0 - \dfrac{(T_0 - T_1) z}{\delta_0}, & 0 < z < \delta_0, \\ T_1, & \delta_0 < z < z_{1}-\delta_{1}, \\ T_1 - \dfrac{(T_1 - T_2)\left(z - \left(z_1 - \delta_1 \right) \right)}{2\delta_1}, & z_{1} - \delta_{1} < z < z_{1}+\delta_{1}, \\ T_{2}, & z_1 + \delta_1 < z < z_{\infty} - \delta_{\infty}, \\ T_{2} - \dfrac{(T_{2} - T_{\infty}) \left(z - (z_{\infty} - \delta_{\infty})\right)}{\delta_{\infty}}, \quad & z_{\infty}-\delta_{\infty} < z, \end{array}\right. \end{equation}

where $z_1$ is the location of the interior flux surface.

The diffusion is minimised by $\delta _j = \ell _j (3 \epsilon )^{1/2}$ and

(7.5a)\begin{gather} T_1 = \frac{2A_0 A_{\infty} \ell_1 T_0 + A_0 A_1 \ell_{\infty} T_0 + A_1 A_{\infty} \ell_0 T_{\infty}}{A_1 A_{\infty} \ell_0 + 2 A_0 A_{\infty}\ell_1 + A_0 A_1 l_{\infty}} \end{gather}

(7.5b)\begin{gather}T_2 = \frac{A_0 A_1 \ell_{\infty} T_0 + A_1 A_{\infty} \ell_0 T_{\infty} + 2 A_0 A_{\infty} \ell_1 T_{\infty}}{A_1 A_{\infty} \ell_0 + 2 A_0 A_{\infty} \ell_1 + A_0 A_1 \ell_{\infty}}, \end{gather}

resulting in the upper bound,

(7.6)\begin{equation} D[T_{\rm trial}] < \frac{(T_{\infty}-T_0)^{2} A_0 A_1 A_{\infty}}{A_1 A_{\infty}l_0 + 2 A_0 A_{\infty} \ell_1 + A_0 A_1 \ell_{\infty}} \sqrt{\frac{\epsilon}{3}} \kappa_{||} . \end{equation}

To numerically test these bounds, we consider a model Hamiltonian which possesses magnetic surfaces at $\rho = 0$, $\rho = 0.5$ and $\rho = 1$,

(7.7)\begin{equation} \chi = \chi_0(\rho) + \mathcal{E} \rho (\rho - 1) (\rho - 0.5) \sum_{n,m} \cos(m \theta - n \zeta) \end{equation}

with $n/m = [1/10,1/5,3/10,2/5,3/5,7/10,4/5,9/10]$. We take our boundaries to be at $\rho = 0$ and $\rho = 1$ and have an internal flux surface at $\rho = 0.5$. Thus, we can associate $\rho$ with the coordinate $z$ appearing in the trial function. We adjust the level of the perturbation amplitude, $\mathcal {E}$, to investigate the effect on the transport as it approaches the upper bound.

We can compute the length scales, $\ell _j$, for the radial component of the field, which satisfy the inequality,

(7.8)\begin{equation} |B_{\rho}| = \left|\frac{\partial \chi}{\partial \theta}\right| \le \mathcal{E} \rho |\rho - 1| |\rho - 0.5| \left|\sum_{m,n} m \sin (m \theta - n \zeta)\right|. \end{equation}

The maximum of the function

(7.9)\begin{equation} f(\theta,\zeta) = \left|\sum_{m,n} m \sin (m \theta - n \zeta)\right| \end{equation}

is numerically computed to be $f_{{\rm max}} \approx 54.7163$ such that

(7.10)\begin{equation} \frac{|B_{\rho}|}{B} \lesssim \frac{\mathcal{E} {\rho}|\rho - 1||\rho-0.5|f_{{\rm max}}}{\sqrt{1+\rho^{2}}}. \end{equation}

In the neighbourhood of $\rho = 0$, $|B_{\rho }|/B \lesssim 0.5\mathcal {E} \rho f_{{\rm max}}$ giving a length scale of $\ell _0 = 2/(\mathcal {E} f_{{\rm max}})$. Similarly, the length scales at $\rho = 0.5$ and $\rho = 1$ are given by $\ell _1 = 2 \sqrt {5}/(\mathcal {E} f_{{\rm max}})$ and $\ell _{\infty } = 2 \sqrt {2}/(\mathcal {E} f_{{\rm max}})$.

We compute the numerical solution to the anisotropic diffusion equation (1.2)–(1.3) using a Fourier Galerkin discretisation in $\theta$ and $\zeta$ and a fourth-order finite-difference radial discretisation. The resulting sparse linear system is solved with the PETSc library using the MUMPS package (Amestoy et al. 2001, 2019) to compute the $LU$-factorisation. In figure 4 we present Poincaré sections of our model field at $\zeta = 0$ in addition to the computed isotherms with $\epsilon = 10^{-6}$. At smaller values of the perturbation amplitude, we see the remnants of island chains in addition to barriers and chaotic regions. The isotherms conform to the structure of the separatrices of the island chains. We note that, as the perturbation amplitude is increased, the magnetic field becomes visually chaotic and the temperature gradient becomes localised near the magnetic surfaces, but away from these the isotherms display a remarkable degree of structure. For a moderately large perturbation amplitude, $\mathcal {E} = 0.1$, the field appears to be strongly stochasticised. Nonetheless, there remain several isotherms which conform to the remnants of secondary island structures near the barriers. Only in the case of the largest perturbation amplitude, $\mathcal {E} = 0.5$, do the pressure gradients become extremely localised near the boundaries. As the trial functions involve complete flattening of the temperature gradient over large regions, we expect to approach the upper transport bound only for sufficiently small $\epsilon$ and sufficiently large $\mathcal {E}$, so that the field is strongly stochastic with dominantly parallel diffusion.

Figure 4. Numerical solutions of the anisotropic diffusion equation with $\epsilon = 10^{-6}$ are overplotted on the Poincaré section of the model magnetic field (7.1) with the displayed values of $\mathcal {E}$: (a) $\mathcal {E} = 0.01$, (b) $\mathcal {E} = 0.05$, (c) $\mathcal {E} = 0.1$ and (d) $\mathcal {E} = 0.5$. For each plot, 50 equally spaced isotherms are displayed (red). The following resolution parameters were used: $m \leq 20$, $|n|\leq 20$, $N_{\rho }=250$ ($\mathcal {E} = 0.01$); $m \leq 20$, $|n|\leq 20$, $N_{\rho }=250$ ($\mathcal {E} = 0.05$); $m \leq 30$, $|n|\leq 30$, $N_{\rho }=600$ ($\mathcal {E} = 0.1$); $m \leq 5$, $|n| \leq 5$, $N_{\rho } = 10^{5}$ ($\mathcal {E} = 0.5$). Here $m$ is the poloidal mode number, $n$ is the toroidal mode number and $N_{\rho }$ is the number of radial grid points.

In figure 5 we compare the numerical value of the diffusion with the upper and lower bounds for $\epsilon = 10^{-6}$, $10^{-7}$ and $10^{-8}$ and several values of the perturbation amplitude $\mathcal {E}$. The lower bound is computed using the minimiser of the isotropic diffusion integral, $D_{{\rm lower}}[\tilde {T}] = ({\kappa _{\perp }}/{2}) \int _{\varOmega } |\boldsymbol {\nabla } \tilde {T}|^{2} \,{{\rm d}}V$, which yields $T_{{\rm lower}} = \rho (T_{\infty }-T_{0})$. For the smallest perturbation amplitude, the total diffusion is within an order of magnitude of the lower bound. As $\epsilon$ decreases, parallel diffusion becomes more dominant, and the difference between $D$ and $D_{{\rm lower}}$ increases marginally. Only for the largest perturbation amplitude ($\mathcal {E} =0.5$) does the total diffusion come within an order of magnitude of the upper bound. There are several factors controlling the difference between $D$ and $D_{{\rm upper}}$. Given the finite value of $\kappa _{\perp }$, there remains a small temperature gradient between the interfaces. The diffusion is further reduced from the upper bound due to the sinusoidal variation of the radial field. Finally, there is a further slight reduction in the diffusion due to the small angular dependence of the isotherms.

Figure 5. The numerical value of the total diffusion (2.1) is compared with the lower (3.6) and upper (7.6) bounds: (a) $\mathcal {E} = 0.01$, (b) $\mathcal {E} = 0.05$, (c) $\mathcal {E} = 0.1$ and (d) $\mathcal {E} = 0.5$.

In figure 6 we choose a specific value of the perturbation amplitude, $\mathcal {E} =0.5$, and compare the numerical solution ($\epsilon = 10^{-6}$) for the temperature profile averaged over $\theta$ and $\zeta$ to the trial function (7.4). We note that the numerical temperature profile is somewhat smoother than the trial function. Nonetheless, several of the features are captured, including the value of the temperature between the interfaces and the location of the strong gradient. We remark that the trial function coincides with the physical temperature profile only for very large perturbation amplitudes ($\mathcal {E} = 0.5$), as significant poloidal and toroidal variation of the isotherms remain at $\mathcal {E} = 0.1$ in addition to a widening of the temperature gradient length scale in comparison with the predicted $\delta _j$.

Figure 6. The numerical solution of the temperature profile (blue) averaged over $\theta$ and $\zeta$ is compared with the trial function (red). The heat diffusion equation is solved with $\epsilon = 10^{-6}$ for the magnetic field (7.1) with perturbation amplitude $\mathcal {E} = 0.5$. (a) The overall temperature profile, which exhibits strong variation in thin boundary layers at $\rho = 0$, $\rho = 0.5$ and $\rho = 1$. These are magnified in (b)–(d).

7.2. Transport with continuously nested isotherms

In this section, we build upon the ideas presented above in § 5 and illustrate the construction of accurate trial functions for the temperature when a continuous family of isotherms can be approximated. The representation of the magnetic field and the boundary conditions are the same as that described in § 7. For the magnetic perturbation, we choose $\chi _{m,n}(\rho ) = \epsilon _{m,n} \rho ^{2} (\rho - 1 )^{2}$, with $\epsilon _{2,1} = 3 \epsilon$, $\epsilon _{3,1} = \epsilon _{3,2} = \epsilon$ and $\epsilon _{4,1}= \epsilon _{4,3} = \epsilon /2$ for $\epsilon =-0.005$. These perturbations vanish on the boundary, and $\rho =0$ and $\rho =1$ are flux surfaces. A Poincaré plot of this field is shown in figure 7. In this section, we take $\epsilon =\kappa _\perp /\kappa _\parallel =10^{-7}$.

Figure 7. (a) Poincaré plot of the magnetic field shown in $(\rho,\theta )$ coordinates, the isotherms of the exact solution shown as black lines (left half only), and the selected QFM surfaces shown as blue lines (right half only). (b) Poincaré plot of the same magnetic field shown in $(\psi,\theta _s)$ coordinates, and the isotherms of the numerically computed temperature shown as black lines (left half only). To a good approximation, the isotherms coincide with surfaces $\psi =const.$

We consider three trial functions that are constructed by first constructing a continuous set of surfaces, which are assumed to approximate isotherms and which we call trial isotherms, and which may be labelled by an arbitrary surface label, $\psi$. For the first of the trial functions, we take the temperature profile to be given by the solution to the isothermal transport case. For the latter two of the trial functions, given the geometry of the trial isotherms, we solve (5.3) for the temperature profile. The accuracy of these trial functions is evaluated by constructing the ‘exact solution’ numerically.

Given that the boundaries of the domain are themselves assumed to be isotherms, the simplest construction of a family of trial isotherms is obtained by suitably interpolating between the boundaries. A particularly relevant interpolation is provided by the contours of the harmonic function $\nabla ^{2}T_{\rm iso}=0$ that equals $T_0$ on $S_0$ and $T_\infty$ on $S_\infty$, as this corresponds to the isotropic diffusion case as discussed in § 3.1. Thus, for our first family of trial isotherms, we take the isotherms of the solution that is obtained when the enhanced parallel transport is ignored, which is obtained when $\epsilon =\kappa _\perp /\kappa _\parallel =1$. For the simple geometry considered here, this solution is $T_{\rm iso}(\rho,\theta,\zeta )=\rho$, which has contours $\rho =const$. Substituting this function directly into the diffusion integral (that is, without changing the geometry of the isotherms and without changing the temperature profile), we obtain $D[T_{\rm iso}]=213.178$.

For the second trial function we assume that the isotherms remain coincident with those of the isothermal solution (in this case, surfaces $\rho =const.$), but we now solve (5.3) with $\epsilon =\kappa _\perp /\kappa _\parallel =10^{-7}$ to obtain the temperature profile. Inserting the trial function for the temperature thus obtained, which we call $\tilde T_{\rm iso}$, into the diffusion integral we obtain $D[\tilde T_{\rm iso}]=89.5803$.

More accurate trial functions can be obtained if, by some means, more accurate approximations to the isotherms can be constructed. Given that, as $\kappa _\parallel$ becomes larger (or as $\kappa _\perp$ becomes smaller), it is the spatial structure of the magnetic field that primarily determines transport. It becomes increasingly plausible that we can set more accurate (i.e. lower) upper bounds on transport across non-integrable magnetic fields by identifying structures that are invariant or almost-invariant under the magnetic field-line flow, such as KAM surfaces and cantori.

A practical numerical approach for constructing invariant and almost invariant surfaces is provided by the construction of (weighted) QFM surfaces (Dewar et al. 1994; Hudson & Suzuki 2014). An algorithm for constructing these surfaces is described in Appendix A. A set of QFM surfaces, defined by their periodicities $\boldsymbol {p}/\boldsymbol {q}=\{ 0/ 1$, $1/ 10$, $1/ 9$, $1/ 8$, $2/ 15$, $1/ 7$, $3/ 20$, $2/ 13$, $3/ 19$, $1/ 6$, $3/ 17$, $2/ 11$, $3/ 16$, $1/ 5$, $3/ 14$, $2/ 9$, $3/ 13$, $1/ 4$, $3/ 11$, $5/ 18$, $2/ 7$, $5/ 17$, $3/ 10$, $4/ 13$, $1/ 3$, $4/ 11$, $3/ 8$, $5/ 13$, $2/ 5$, $5/ 12$, $3/ 7$, $7/ 16$, $4/ 9$, $5/ 11$, $1/ 2$, $6/ 11$, $5/ 9$, $9/ 16$, $4/ 7$, $7/ 12$, $3/ 5$, $8/ 13$, $5/ 8$, $7/ 11$, $2/ 3$, $7/ 10$, $12/ 17$, $5/ 7$, $8/ 11$, $3/ 4$, $7/ 9$, $11/ 14$, $4/ 5$, $9/ 11$, $5/ 6$, $11/ 13$, $6/ 7$, $13/ 15$, $7/ 8$, $8/ 9$, $9/ 10$, $1/ 1 ) \}$, are shown as the blue curves in the left figure in figure 7. (How these periodicities were chosen is described further in the following.)) These surfaces are taken as coordinate surfaces to construct the coordinate mapping $\boldsymbol {x}(\psi,\theta _s,\zeta )$, where $\psi$ labels the enclosed toroidal flux and $\theta _s$ is a straight pseudo-field-line angle. An interpolation between the surfaces provides a continuous coordinate framework for non-integrable magnetic fields, which reduces to straight-field-line coordinates where flux surfaces exist. Much of the field-line dynamics that remains invariant after perturbation can be incorporated into the coordinates, and we call these coordinates chaotic coordinates (Hudson & Suzuki 2014). The Poincaré plot of the same magnetic field as shown in the left plot of figure 7 is shown in chaotic coordinates in the right plot of figure 7.

We take as our third family of trial isotherms to be the QFM surfaces and their radial interpolates, where the temperature profile is given by the solving (5.3). For the trial function, $T=T_{\boldsymbol {p}/\boldsymbol {q}}$, thus obtained, the diffusion is $D[T_{\boldsymbol {p}/\boldsymbol {q}}] = 24.8809$.

We may compute the exact solution using a mixed finite-difference Fourier representation for the temperature, $T(\rho _i,\theta,\zeta )=\sum _{m=0}^{M}\sum _{n=-N}^{N}T_{m,n,i} \cos (m\theta -n\zeta )$ including Fourier modes up and including $M=10$ and $N=10$ with a second-order finite-difference approximation for the radial derivatives and $256$ radial grid points, and minimise $D[T]$ with respect to the $T_{m,n,i}$. The numerically constructed temperature, $T_e$, which for convenience we consider to be the exact solution, is shown as the black curves in the right-hand figure in figure 7. Note that, visually, the isotherms coincide with coordinate surfaces to a very good approximation. With the accurate approximation of the exact temperature, $T_e$, obtained numerically, the diffusion is $D[T_e]=22.8721$.

We can compare the values of $D[T_{\rm iso}]$, $D[\tilde T_{\rm iso}]$ and $D[T_{\boldsymbol {p}/\boldsymbol {q}}]$ to $D[T_{e}]$. We see that $(D[T_{\rm iso}]-D[\tilde T_{\rm iso}])/(D[T_{\rm iso}]-D[T_e])\approx 0.65$, which means that $65\,\%$ of the total possible reduction of the diffusion integral (starting from the isothermal solution) is achieved merely by solving for the isothermal solution to obtain a family of trial isotherms and then solving (5.3) for the temperature profile. (Note that solving for the isothermal solution is a much simpler numerical task than solving for the exact solution, because the isothermal solution only depends on the geometry of the boundary and not on the magnetic field itself.) For the trial function based on the QFM surfaces and their interpolates, and solving (5.3), we obtain $(D[T_{\rm iso}]-D[T_{\boldsymbol {p}/\boldsymbol {q}}])/(D[T_{\rm iso}]-D[T_e])\approx 0.99$. Based on this, the trial function based on chaotic coordinates and solving (5.3) yields an improvement of $99\,\%$ in the sense that over $99\,\%$ of the total possible reduction of the diffusion integral resulting from including the strong parallel transport has been achieved. It is much faster to construct QFM surfaces than to solve the anisotropic diffusion equation numerically.

The selection of the QFM surfaces requires some explanation. QFM surfaces are characterised by their periodicity, and we may distinguish the following cases.

1. (i) For $(p,q)$ that are low-order rationals, e.g. $p/q=1/2, 1/3, \ldots,$ because the low-order islands are typically larger than the high-order islands and because the $(p,q)$ QFM surface pass through the $(p,q)$ island, we may expect that low-order $(p,q)$ QFM surfaces are isotherms because the temperature tends to flatten inside the larger islands.

2. (ii) For $(p,q)$ that are high-order rationals that approximate a strongly irrational number, then the QFM surface is an excellent approximation to the KAM surface if it exists, or if it does not the QFM surface effectively ‘fills in the gaps’ in the cantorus. In both cases, because complete barriers and partial barriers to magnetic field-line transport are likely to be isotherms, then these QFM surfaces are likely to be isotherms too.

3. (iii) A particularly interesting case is when $p/q$ is a high-order rational that approximates a low-order rational, e.g. $p/q=101/200$. In this case, the QFM surface lies just outside of the low-order island separatrix, and whether this QFM surface is an isotherm depends on how large the low-order island is and on $\kappa _\perp$.

The question thus arises: which QFM surfaces are most likely to coincide with isotherms? The answer partly depends on the structure of the non-integrable magnetic field, namely whether the magnetic field is completely integrable, nearly integrable or strongly chaotic, and also on how much the parallel transport dominates the perpendicular. A practical algorithm, used in this study, is provided by constructing successively sophisticated trial functions and evaluating the diffusion integral. If the diffusion integral decreases, then the trial function is a more accurate approximation of the true solution.

We start from a trial function, $T_{\boldsymbol {p}/\boldsymbol {q}}$, built on a family of trial isotherms constructed by an interpolation of a selection of a few, low-order QFM surfaces (in this example we may begin with the $\boldsymbol {p}/\boldsymbol {q} = \{0/1, 1/2, 1/1\}$ set of QFM surfaces) and solve for the temperature profile given in (5.3), and then compute $D[T_{\boldsymbol {p}/\boldsymbol {q}}]$. Next, we construct a refined trial function by augmenting the previous selection of QFM surfaces by, typically, adding mediants, to obtain $(p_{i}+p_{i+1})/(q_{i}+q_{i+1})$, to obtain, for example, the set of $\boldsymbol {p}/\boldsymbol {q} = \{0/1, 1/2, 1/2, 2/3, 1/1\}$ QFM surfaces, recompute the interpolation to construct the family of trial isotherms, re-solve (5.3) and then recompute the diffusion. Adding more and higher-order QFM surfaces effectively adds more of the fractal structure of the magnetic field to the trial function. If the diffusion integral for the refined trial function decreases, then a better (lower) upper bound on transport has been obtained and the refined trial function is a better approximation to the exact solution. Additional mediants can be included, one at a time if need be, until the diffusion integral no longer decreases, and this may occur when the scale length of the QFM surfaces is smaller than the scale length of the true solution, which is in part set by $\epsilon =\kappa _\perp /\kappa _\parallel$. In this fashion, we may systematically improve our estimate for the global temperature solution, and in so doing identify the most important barriers to anisotropic transport for non-zero $\epsilon$ in non-integrable magnetic fields. There are other potential algorithms for identifying the most suitable QFM surfaces to be used as trial isotherms, and we shall explore these ideas in future work.

8. Discussion

Magnetic fields with chaotic field lines can possess a surprising degree of structure, which is not readily visible in a Poincaré puncture plot but becomes apparent when one for instance considers anisotropic transport processes. In such situations, traditional estimates (Rechester & Rosenbluth 1978; Kadomtsev & Pogutse 1979; Krommes et al. 1983) of the ‘effective’ diffusion coefficient are unreliable and one must instead take the detailed geometry of the magnetic field into consideration. KAM surfaces, and remnants thereof, act as barriers to the diffusion, and their effect can be quantified through rigorous upper and lower bounds on the net transport. A more restrictive upper bound can be found whenever the isothermal surfaces can be guessed to good accuracy. Such surfaces can be constructed in natural by using ‘chaotic coordinates’, which also find wide application in other parts of plasma physics and in Hamiltonian dynamical systems.