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.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 (Reference Hudson and Breslau2008). 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 Reference Meiss1992). Magnetic field lines pass through partial barriers, unlike flux surfaces, but they do so slowly in some sense.Footnote ¹ 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 (Reference Dewar, Hudson and Price1994) and Hudson & Dewar (Reference Hudson and Dewar2009).Footnote ² 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).

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. Reference Amestoy, Duff, Koster and L'Excellent2001, Reference Amestoy, Buttari, L'Excellent and Mary2019) 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. Reference Dewar, Hudson and Price1994; Hudson & Suzuki Reference Hudson and Suzuki2014). 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 Reference Hudson and Suzuki2014). 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.