2 Anisotropic heat conduction

We proceed from an anisotropic heat conduction equation, with a temperature-dependent sink that models radiation losses in and around a magnetic island. As we shall see, the full 3-D heat conduction equation can, in the limit of highly anisotropic heat conductivity, be reduced to an ordinary differential equation. This reduction is achieved by flux-surface averaging the heat conduction equation, and results in a 1-D description with an effective heat conduction coefficient that carries information about the magnetic geometry.

We consider a heat conduction equation of the form

(2.1)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol q} + R(T) = 0, \end{equation}

where ${\boldsymbol q}$ denotes the heat flux and $R(T)$ the energy loss due to radiation, which depends on the temperature $T$. For simplicity we ignore any dependence on density, and assume that the plasma is optically thin to the radiation, so that the latter simply acts as an energy sink. In the classical theory of collisional transport in magnetised plasmas (Braginskii Reference Braginskii1965; Helander & Sigmar Reference Helander and Sigmar2005), the electron heat flux is equal to

(2.2)\begin{equation} {\boldsymbol q} ={-} \kappa_\| \boldsymbol{\nabla}_\| T - \kappa_\wedge {\boldsymbol b} \times \boldsymbol{\nabla} T - \kappa_\perp \boldsymbol{\nabla}_\perp T, \end{equation}

where $T$ is the electron temperature, ${\boldsymbol b} = {\boldsymbol B} / B$ a unit vector pointing in the direction of the magnetic field, $\boldsymbol {\nabla }_{\|}= \boldsymbol b \boldsymbol b \boldsymbol {\cdot } \boldsymbol {\nabla }$, and $\boldsymbol {\nabla }_{\perp }{T} = \boldsymbol {\nabla }{T} -\boldsymbol { b} \boldsymbol { b} \boldsymbol {\cdot } \boldsymbol {\nabla }{T}$. The three different heat conductivities that appear in this relation are of very different magnitude in a strongly magnetised plasma. Heat conduction is far more effective in the direction along the field than in the two directions perpendicular to it, so we adopt the ordering

(2.3)\begin{equation} \kappa_\| \gg \kappa_\wedge{\sim} \kappa_\perp. \end{equation}

In Braginskii's theory, the diamagnetic heat flux is much larger than the perpendicular one, $\kappa _\wedge \gg \kappa _\perp$, but in practice the latter is never as small as this theory predicts because of additional turbulent transport. We therefore take the perpendicular and diamagnetic transport to be comparable and expand the temperature $T = T_0 + T_1 + \cdots$ accordingly. At zeroth order, we then find

(2.4)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \left( \kappa_\| \boldsymbol{\nabla}_\| T_0 \right) = R(T_0). \end{equation}

We assume that the magnetic field traces out topologically toroidal flux surfaces. If magnetic islands are present, these surfaces are not simply nested, but surfaces are still formed by the field lines unless the latter are chaotic due to island overlap. We assume that only a single chain of islands is present and that the field is thus integrable. (Even if it is not, similar considerations may apply thanks to the presence of so-called ghost surfaces (Hudson & Breslau Reference Hudson and Breslau2008).) It is customary to define the flux-surface average $\langle \cdots \rangle$ of any function as the volume average of this function over the region between two neighbouring flux surfaces (Helander Reference Helander2014). If $\psi$ is some flux-surface label, i.e. a coordinate that is constant on each flux surface and varies smoothly between them, then the flux-surface average of the function $f({\boldsymbol r})$ is thus

(2.5)\begin{equation} \langle f \rangle =\left. \int_{dU} f({\boldsymbol r}) \,{\rm d}{\boldsymbol r} \right/ \int_{dU} \,{\rm d}{\boldsymbol r}, \end{equation}

where ${\rm d}U$ denotes the region between the flux surfaces labelled by $\psi$ and $\psi + {\rm d}\psi$, where ${\rm d}\psi \rightarrow 0$. It follows from Gauss’ law that $\langle \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol F} \rangle = 0$ vanishes for any single-valued vector field ${\boldsymbol F}({\boldsymbol r})$ that is tangential to flux surfaces. The flux-surface average of (2.4) thus gives $\langle R(T_0) \rangle =0$, which indicates that $R(T)$ should be relegated to higher order. The zeroth-order equation then becomes

(2.6)\begin{equation} {\boldsymbol B} \boldsymbol{\cdot} \boldsymbol{\nabla} \left( \frac{\boldsymbol{\nabla}_\| T_0}{B} \right) = 0, \end{equation}

and it follows that $T_0$ cannot vary over flux surfaces, i.e. $T_0$ is some function of $\psi$ alone that is determined by the higher-order equations. We thus proceed to the next order, where we find

(2.7)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \left( \kappa_\| \boldsymbol{\nabla}_\| T_1 + \kappa_\wedge {\boldsymbol b} \times \boldsymbol{\nabla} T_0 + \kappa_\perp \boldsymbol{\nabla}_\perp T_0 \right) = R(T_0). \end{equation}

A flux-surface average of this equation conveniently annihilates the first two terms, and thus produces an equation,

(2.8)\begin{equation} \langle \boldsymbol{\nabla} \boldsymbol{\cdot} \left( \kappa_\perp \boldsymbol{\nabla}_\perp T_0 \right) \rangle = R(T_0), \end{equation}

which, together with the boundary conditions, determines the lowest-order temperature $T_0(\psi )$. If the magnetic surfaces do not exist, as is the case for chaotic magnetic fields (Hudson & Breslau Reference Hudson and Breslau2008), we note that the flux-surface average can be replaced by an appropriate average over ghost surfaces. As a result, the term representing parallel heat conduction is then not entirely annihilated. Since the loss function $R$ is assumed only to depend on the temperature $T$, its flux-surface average is independent of the flux surface on which the radiation occurs, i.e. $\langle R(T_0) \rangle =R(T_0)$. Thanks to the large anisotropy of the heat conduction, we have thus reduced the original 3-D partial differential equation (2.1) to an ordinary 1-D differential equation. Indeed, if we supplement our flux-label coordinate $\psi$ with two $2{\rm \pi}$-periodic coordinates $(\theta,\varphi )$ that specify the location within each magnetic surface, then (Hazeltine & Meiss Reference Hazeltine and Meiss2003; Helander Reference Helander2014; Helander, Hudson & Paul Reference Helander, Hudson and Paul2022; Paul et al. Reference Paul, Hudson and Helander2022)

(2.9)\begin{equation} \langle \boldsymbol{\nabla} \boldsymbol{\cdot} \left( \kappa_\perp \boldsymbol{\nabla}_\perp T_0 \right) \rangle = \frac{1}{V'(\psi)} \frac{{\rm d}}{{\rm d}\psi} \left( V'(\psi) \langle \kappa_\perp |\boldsymbol{\nabla} \psi|^2 \rangle \frac{{\rm d} T_0}{{\rm d}\psi} \right), \end{equation}

where

(2.10)\begin{equation} V'(\psi) = \int_0^{2 {\rm \pi}} \int_0^{2 {\rm \pi}} \frac{{\rm d}\theta\, {\rm d}\varphi}{(\boldsymbol{\nabla} \psi \times \boldsymbol{\nabla} \theta) \boldsymbol{\cdot} \boldsymbol{\nabla} \varphi} \end{equation}

is the derivative of the enclosed flux surface volume $V(\psi )$. Specifically, if we choose $V(\psi )$ as our flux-surface label instead of $\psi$, our equation (2.8) becomes simply

(2.11)\begin{equation} \frac{{\rm d}}{{\rm d}V} \left( \langle \kappa_\perp |\boldsymbol{\nabla} V|^2 \rangle \frac{{\rm d} T}{{\rm d}V} \right) = R(T), \end{equation}

where we have dropped the subscript on $T_0$. If $\kappa _\perp (T)$ is temperature-dependent, this dependence be eliminated from $\langle \kappa _\perp |\boldsymbol {\nabla } V|^2 \rangle$ though the definition of a new temperature

(2.12)\begin{equation} \tilde{T}=\int_{0}^{T} \kappa_\perp (T^{\prime}) \,{\rm d} T^{\prime}, \end{equation}

and a redefinition of $R(T)$.

In the 1-D conduction equation (2.11), the function $\kappa (V)=\langle |\boldsymbol {\nabla } V|^2 \rangle$ plays the role of an effective heat conduction coefficient, which encapsulates the influence of magnetic flux expansion on the perpendicular heat conduction. Our next step will be to calculate $\kappa (V)$ in a field with a single magnetic island chain. As we shall see, the fact that the topology of the magnetic field lines changes across the island separatrix implies a non-trivial topology of the set in which (2.11) is solved. Although this is an ordinary differential equation in terms of the independent variable $V$ (the enclosed flux-surface volume), the latter does not have the usual topology of the real line.

3 Heat conduction equation near a magnetic island

3.1 Magnetic island geometry

Magnetic islands can arise in toroidal magnetic fields with nested flux surfaces in places where the rotational transform is a rational number. A resonant perturbation to the field-line Hamiltonian $\chi$ describing the field without the island in general produces such an island. If the perturbation only consists of a single Fourier harmonic, the perturbed field is described through a flux function similar to that of a pendulum Hamiltonian. From this analytical description of the island, it is possible to derive an expression for $\kappa (V)$ in terms of complete elliptic integrals.

A toroidal magnetic field can be described by

(3.1)\begin{equation} \boldsymbol{B}=\boldsymbol{\nabla} \psi \times \boldsymbol{\nabla} \theta+\boldsymbol{\nabla} \varphi \times \boldsymbol{\nabla} \chi_{0}. \end{equation}

If all flux surfaces are simply nested, then $\chi _0$ is a function of $\psi$ alone. We can thus write $\boldsymbol {\nabla } \chi _{0}=\iota (\psi ) \boldsymbol {\nabla }\psi$, where $\iota (\psi ) = d\chi _0/d\psi$, and Clebsch coordinates can be defined by $\alpha = \theta -\iota \varphi$, so that ${\boldsymbol B} = \boldsymbol {\nabla } \psi \times \boldsymbol {\nabla } \alpha$. The rotational transform $\iota (\psi )$ describes the number of poloidal turns per toroidal turn along a magnetic field line. Its derivative $\iota ^{\prime }(\psi )$ is known as the magnetic shear and describes how the rotational transform varies between flux surfaces. As already mentioned, magnetic islands on rational surfaces arise in response to a resonant perturbation $\delta \chi$ to the Hamiltonian system, and it is convenient for the description of a magnetic island to define a new effective Hamiltonian $\chi =\chi _{0} + \delta \chi -\iota \psi$. The effective Hamiltonian facilitates the analysis by simplifying the system in a way similar to the description in Cary & Hanson (Reference Cary and Hanson1991) and Lichtenberg & Lieberman (Reference Lichtenberg and Lieberman2013). Expanding the perturbation $\delta \chi$ in a Fourier series in the angular variables $\theta$ and $\varphi$, and keeping only one term in the expansion gives

(3.2)\begin{equation} \chi=\frac{1}{2}\iota^{\prime} x^{2} - F \cos\alpha, \end{equation}

where $w$ measures the island width, $x$ the deviation in $\psi$ from the rational surface and $F={w^{2}\iota ^{\prime }}/{4}$. For a detailed derivation see (Helander Reference Helander2014) and (Hazeltine & Meiss Reference Hazeltine and Meiss2003). The contours of $\chi$ are shown in figure 1. Next, we introduce the normalised island flux $m=\tfrac {1}{2}({\chi }/{F}+1 )$:

(3.3)\begin{equation} m=\frac{1}{4 \varOmega^{2}} x^{2} + \sin^{2}\vartheta. \end{equation}

Here we have introduced the island half-width $\varOmega =({F}/{\iota ^{\prime }})^{{1}/{2}}= {w}/{2}$ and the new angle $\vartheta ={\alpha }/{2}$. The normalised flux has the important property of lying between $0< m<1$ inside the island and $0<{1}/{m}<1$ outside the island. For simplicity, we assume that unperturbed magnetic field is homogeneous and neglect any spatial variation of $|\boldsymbol {\nabla } x|$, thus effectively considering the geometry of a plasma slab with a magnetic island.

Figure 1. Contours of the magnetic flux $\chi$ for a magnetic field with an island and $\varOmega ={1}/{4}$.

3.2 Effective heat conductivity

We proceed to determine the effective heat conduction coefficient in the magnetic island topology. From the 1-D heat conduction equation (2.11), it is clear that an expression for

(3.4)\begin{equation} {\left|\boldsymbol{\nabla}{V}\right|}^{2}=\left( \frac{{\rm d}V}{{\rm d}m} \right)^{2} {\left|\boldsymbol{\nabla}{m}\right|}^{2}, \end{equation}

needs to be determined from our normalised pendulum Hamiltonian (3.3). We use $m$ as a flux-surface label and find, with the definition of the $V^\prime$ given earlier,

(3.5)\begin{equation} \frac{{\rm d}V}{{\rm d}m} = \begin{cases} \dfrac{4 \varOmega}{\sqrt{m}}K(m^{{-}1}), & m>1 \\ 8 \varOmega K(m) , & m<1 , \end{cases} \end{equation}

where $K$ and $E$ are complete elliptic integrals of the first and second kind, respectively. The argument may take any value between zero and unity. At the island O-point and far from the island $m \rightarrow 0$, whereas $m \rightarrow 1$ at the separatrix. To determine $|\boldsymbol {\nabla } m|^{2}$ we express $x^{2}$ in (3.3) as a function of $\vartheta$ and $m$ itself, and thus find

(3.6)\begin{equation} {\left|\boldsymbol{\nabla}{m}\right|}^{2} = \frac{1}{\varOmega^{2}} \left( m -\sin^{2}\vartheta\right) +\frac{1}{4} \sin^{2} 2\vartheta. \end{equation}

Combining the expressions for ${\left |\boldsymbol {\nabla }{ m}\right |}^{2}$ and ${{\rm d} V}/{{\rm d} m}$ using $\sin ^{2} 2\vartheta = 4 ( \sin ^{2} \vartheta - \sin ^{4}\vartheta )$, an expression for ${\left |\boldsymbol {\nabla }{V}\right |}^{2}$ inside the island and outside the island can be found. The integrals in the flux surface averaging can again be reduced to combinations of $K$ and $E$. For details on elliptic integrals see Byrd & Friedman (Reference Byrd and Friedman2012) and Gradshteyn & Ryzhik (Reference Gradshteyn and Ryzhik2014). Proceeding in this fashion for $m<1$ and $m>1$, respectively, analytical expressions for the heat conduction coefficient $\kappa$ can be found as a function of $m$. The result,

(3.7)\begin{equation} \kappa(m ) = \begin{cases} \begin{aligned} & \dfrac{16 \varOmega^{2} K(m^{{-}1})}{3} \\ & \left[ 2\left( 1-m\right) K(m^{{-}1}) +\left( 2m -1+\frac{3}{\varOmega^{2}}\right) E(m^{{-}1}) \right]\end{aligned} & m>1,\\ \begin{aligned} & \dfrac{8^{2} \varOmega^{2} K(m)}{3} \\ & \left[\left( m-1\right)\left( \frac{3}{\varOmega^{2}}-1 \right)K(m)+\left(2 m +\frac{3}{\varOmega^{2}}-1 \right)E(m) \right] \end{aligned} & m<1, \\ \end{cases} \end{equation}

can be combined into a single expression that is valid for all values of $m$:

(3.8)\begin{equation} \kappa(m)=\frac{1}{6} \left[ \left(\frac{3}{\varOmega^{2}}+2 m-1\right) \frac{{\rm d}V^{2}}{{\rm d}m}+ 4\left(m-m^{2} \right)\left(\frac{{\rm d}V}{{\rm d}m} \right)^{2}\right]. \end{equation}

As we shall see below, inside the island $\kappa (m)$ is well approximated by the first term in the square brackets, whereas both terms are required for an accurate description far outside the island. As is evident from the numerical plot given in figure 2, when the separatrix is approached from either side, $\kappa$ goes to infinity because of the dependence on $K(m)$.

Figure 2. The equivalent 1-D heat conduction coefficient over the range of $0< m<3$. Outside the separatrix, $\kappa$ quickly approaches the limiting value for $m\gg 1$. Here $\kappa$ is normalised to $\kappa _{\infty }$.

An analytical expression can also be given for $V(m)$,

(3.9)\begin{equation} V = \begin{cases} 8 \varOmega \sqrt{m} E(m^{{-}1}), & m>1, \\ 16 \varOmega \left[ E(m)+(m-1)K(m) \right], & m<1 , \end{cases} \end{equation}

which is plotted in figure 3. The enclosed flux surface volume for the pendulum Hamiltonian is related to the action integral $J = V/2{\rm \pi}$, which is well known in Hamiltonian mechanics (Rechester & Stix Reference Rechester and Stix1979; Lichtenberg & Lieberman Reference Lichtenberg and Lieberman2013).

Figure 3. The enclosed flux-surface volume in and around the island and $V^{\prime }={{\rm d}V}/{{\rm d}m}$ with $\varOmega ={1}/{4}$. The discontinuity in $V$ at the separatrix ($m=1$) is due to the fact that there are two regions ‘above’ and ‘below’ the island for $m>1$, each enclosing half of the total enclosed flux surface volume around the island. To understand the impact of competing $V$ and $V^{\prime }$ terms, $V V^{\prime }$ is also shown.

At this point, it may be helpful to recapitulate the various coordinate transformations that have been performed. The unperturbed magnetic field (3.1) was originally written in general toroidal coordinates $(\psi,\theta,\varphi )$, and a magnetic island was added in a standard way through the introduction of the helical flux $\chi$ in (3.2). On the perturbed flux surfaces, which contain a periodic chain of magnetic islands, it proved useful to introduce the surface label $m$ through (3.3) and the enclosed volume $V(m)$ = (3.9), in terms of which the effective heat conductivity (3.8) can be expressed.

For most choices of the radiation function $R(T)$, the differential equation (2.11) can only be solved numerically. To simplify the problem, we shall take $R(T)$ to be piecewise constant, in which case analytical progress is possible.

Regardless of the choice of $R(T)$, however, this equation has the interesting feature of operating in a domain of non-trivial topology. Although it is an ordinary differential equation, $T(V)$ is not to be understood as a function along a line but on the set with the topology shown in figure 5. The topology is similar to a rod with a cooling fin, where the latter corresponds to the interior of the island and the heat conductivity depends on the position along the rod and the cooling fin.

3.3 Limiting cases

There are certain regions in which the heat conduction equation becomes particularly simple, namely, in the vicinity of the island O-point, at the separatrix and far away from the island. As one would intuitively expect, far from the island the equation reduces to an equation describing perpendicular heat conduction across flux surfaces without an island. Close to the O-point, the heat conduction equation is equivalent to that in an ellipse, and, as we shall see, the effective heat conduction coefficient $\kappa$ reaches a maximum at the separatrix.

For the case far away from the island, we note that $m\gg 1$ and thus $\kappa \to \kappa _{\infty }=4{\rm \pi} ^2$ becomes independent of $V$, as expected in a straight magnetic field. Equation (2.11) then reduces to

(3.10)\begin{equation} \frac{{\rm d}^{2}T}{{\rm d}V^{2}} = R(T), \end{equation}

where we have absorbed $4{\rm \pi} ^{2}$ into $R(T)$.

In order to treat the region close to the O-point, we expand $\kappa$ and $V$ for small $m$, $\kappa (V) =2 {\rm \pi}\varOmega V ( 1+ {1}/{\varOmega ^{2}}) +O{m^{2}}$, and note that (2.11) reduces to a heat conduction equation in elliptical geometry,

(3.11)\begin{equation} 2 {\rm \pi}\varOmega \left( 1+ \frac{1}{\varOmega^{2}} \right) \frac{{\rm d}}{{\rm d}V} \left( V\frac{{\rm d} T}{{\rm d}V} \right) = R(T). \end{equation}

The equation of the corresponding ellipse can be obtained by expanding equation (3.3) in $\vartheta$, resulting in $4 m =({x}/{\varOmega })^{2}+ (2 \vartheta )^{2}+O{\vartheta ^{3}}$. Thus, the ellipticity is given by the half-width of the island $\varOmega$, the radius of the ellipse is given through $r=2\sqrt {m}$ and the enclosed volume near the O-point the enclosed flux surface volume dominates the heat conductivity.

In the vicinity of the separatrix, $\kappa$ can be expanded in $x\ll 1$ with $m=1-x$ inside the island and $m={1}/{(1-x)}$ outside the island. In both cases, we find that the effective heat conduction coefficient is given by

(3.12)\begin{equation} \kappa= \frac{1}{6}\left(\frac{3}{\varOmega^{2}}+1\right) V V^{\prime} +O{x} , \end{equation}

where $V^{\prime }={{\rm d}V}/{{\rm d}m}$.

These approximations of the effective heat conductivity are asymptotically exact in the respective limits and less accurate elsewhere. The relative error ${{\left |\kappa _{{\rm app}}-\kappa \right |}}/{\kappa }$ is shown in figure 4. If the full expressions are used for the quantities $V$ and $V^{\prime }$ in (3.12),

(3.13)\begin{equation} \kappa_{{\rm app}}=\frac{1}{6}\left(\frac{3}{\varOmega^{2}}+1\right) V V^{\prime} , \end{equation}

the resulting expression yields a remarkably good approximation everywhere in the domain, with an error less than a few per cent in figure 4. From this approximation, we conclude that the product $V V^{\prime }$ is the main contributing factor to the variation of the effective heat conductivity over the domain.

Figure 4. The deviation of the 1-D heat conduction coefficient from its asymptotic form in various limits. The expansion around the O-point and the separatrix (given as $\kappa _{O}$ and $\kappa _{{\rm sep}}$), have very large errors farther away from their expansion point. While $\kappa _{{\rm app}}$ does not deviate more than $0.04$ from $\kappa (m)$ in the relevant region in and around the island.

4 Variational principle

The effective heat conductivity across the magnetic field $\kappa (m)$ approaches infinity at the separatrix, which contradicts the ordering used to derive the 1-D heat conduction equation (2.11) in § 2. Accordingly, there is a boundary layer around the separatrix in which the temperature varies along the magnetic field and is not properly described by the ordinary differential equation (2.11). Instead, a 2-D treatment is necessary in this region, and one may ask whether this boundary layer may somehow affect the global properties of the solution. This seems unlikely since the boundary layer is narrow in the limit

(4.1)\begin{equation} \epsilon = \frac{\kappa_\perp}{\kappa_\|} \rightarrow 0, \end{equation}

but it is perhaps nevertheless of interest to demonstrate that the solution to (2.11) indeed converges to that of the original 2-D equation (2.1) in this limit.

An argument based on a variational principle can be constructed if the diamagnetic heat flux is neglected, $\kappa _\wedge = 0$. Then the trial function $T$ that minimises the functional

(4.2)\begin{equation} D[T] = \int \left[ (\boldsymbol{\nabla}_\| T)^2 + \epsilon (\boldsymbol{\nabla}_\perp T)^2 + \epsilon r(T) \right] \,{\rm d}{\boldsymbol r}, \end{equation}

with $r'(T) = 2 R(T) / (\epsilon \kappa _\|)$, solves our original partial differential equation (2.1) if the temperature or its normal derivative is prescribed on the boundary. With these boundary conditions, (2.1) is satisfied if, and only if, $\delta D[T] = 0$. In particular, the function $T({\boldsymbol r})$ that minimises the functional $D[T]$ satisfies this equation.

Equation (2.11) also corresponds to a variational principle, namely $\delta Q[T_0] = 0$, with

(4.3)\begin{equation} Q[T_0] = \int \left[ \langle |\boldsymbol{\nabla} V |^2 \rangle \left( \frac{{\rm d}T_0}{{\rm d}V} \right)^2 + k(T_0) \right] \,{\rm d}V, \end{equation}

where $k'(T_0) = 2 R(T_0) / \kappa _\perp$. This functional is closely related to $D$: for functions $T({\boldsymbol r})$ such that $\boldsymbol {\nabla }_\| T_0 = 0$ they are proportional to each other, $D[T] = \epsilon Q[T]$. In other words, whereas (2.1) is satisfied by the function $T$ that minimises $D[T]$ under the constraint that the boundary conditions are satisfied, the reduced equation (2.11) is satisfied by the function $T_0$ that minimises $D[T_0]$ under the additional constraint that $\boldsymbol {\nabla }_\| T_0 = 0$. It follows that $D[T] \le D[T_0]$.

Having established this variational principle, we now proceed to argue that the function $T({\boldsymbol r})$ that solves (2.1) in the limit of large anisotropy in the heat conductivity is constant along magnetic-field lines, i.e.

(4.4)\begin{equation} \lim_{\epsilon \rightarrow 0+} \boldsymbol{\nabla}_\| T = 0. \end{equation}

If this were not the case, it would be possible to find two points, ${\boldsymbol r}_P$ and ${\boldsymbol r}_Q$, on the same field line and a number $\tau > 0$ such that $T({\boldsymbol r}_P) - T({\boldsymbol r}_Q) > \tau$ for all sufficiently small values of $\epsilon >0$. By continuity, the same would be true for points in a neighbourhood of ${\boldsymbol r}_P$ connected by field lines to a neighbourhood of ${\boldsymbol r}_Q$. Consider a ‘flux tube’ $F$ connecting these neighbourhoods; then

(4.5)\begin{equation} D[T] \ge \int_F (\boldsymbol{\nabla}_\| T)^2 \,{\rm d}{\boldsymbol r}, \end{equation}

where the right-hand side can be bounded from below (Helander et al. Reference Helander, Hudson and Paul2022). The bound is proportional to $\tau$ and depends on the geometry of the flux tube but is independent of $\epsilon$, and it follows that $D[T]$ is similarly bounded from below in the limit $\epsilon \rightarrow 0+$.Footnote ¹ But this leads to a contradiction since then $D[T] > D[T_0]$ for the function $T_0$ satisfying (2.11). We thus conclude that $T({\boldsymbol r})$ cannot vary along flux tubes of non-vanishing volume in the limit $\epsilon \rightarrow 0+$, and that $T({\boldsymbol r}) \rightarrow T_0({\boldsymbol r})$.

5 Heat transport in the equivalent 1-D geometry

We now proceed to solve the heat conduction equation (2.11) in a region containing the magnetic island and part of the surrounding plasma. We begin by appropriately normalising our variables and introducing a specific radiation function $R(T)$. This function is defined as follows:

(5.1)\begin{equation} R(T)=\begin{cases} R_{0}, & T_{0}< T< T_{1},\\ 0 , & \text{otherwise}, \end{cases} \end{equation}

and ensures that radiation occurs only within a well-defined zone, bounded by the temperatures $T_{0}$ (cold side) and $T_{1}$ (hot side). The choice of this radiation function is driven by simplicity and aligns with the behaviour of line radiation from impurities in a typical plasma edge. The constant $R_{0}$ controls the radiation amplitude and is proportional to the impurity fraction and the square of the plasma density $n_{e}^{2}$.

There are several possibilities for the radiation front's position, as depicted in figure 5. The radiation zone can be situated on one side of the island, around the separatrix, below the island and between the separatrix and the O-point, as well as directly at the O-point. The region outside the island is taken to be bounded in one direction by a ‘wall’ with zero temperature and to extend some distance in the other direction, which in a tokamak or a stellarator would correspond to the far interior of the core plasma.Footnote ² As already described, the topology of the 1-D set on which differential equation is to be solved is that of a rod with a cooling fin, see figure 5. Boundary conditions need to be prescribed at the wall, at infinity and at the O-point of the island. At the latter position, the heat flux must vanish, while the temperature is unconstrained. At the wall we impose the boundary condition $T=0$, as the wall temperature is significantly lower than that of the plasma. At the other (hot) end of the domain, the incoming heat flux is prescribed.

Figure 5. Equivalent 1-D topology of the heat conduction in a magnetic field with an island and the different positions of the radiation front. The topology is similar to a rod with a cooling fin. Here the interior of the island corresponds to the cooling fin and the heat conductivity depends on the position along the rod and inside the cooling fin. The radiation zone can be situated on one side of the island, around the separatrix, below the island and between the separatrix and the O-point, as well as directly at the O-point.

It is convenient to normalise the variables in our 1-D heat conduction equation by writing

(5.2a–e)\begin{equation} \hat{T}=\frac{T}{T_{1}} {,}\quad \hat{x}=\frac{V}{w} {,}\quad \hat{\kappa}(x)=\frac{\kappa(V)}{\kappa({\rm d} w)} {,}\quad \hat{q}_{{\rm in}}=\frac{q_{{\rm in}}} w{\kappa({\rm d} w) T_{1}} {,}\quad \widehat{R_{0}}=\frac{R_{0} w^2}{\kappa({\rm d} w) T_{1}} {,} \end{equation}

where $q_{{\rm in}}$ denotes the incoming heat flux, i.e. the heat flux at $x_d$, and we have chosen $w=4$ in the example shown above. Through the definition of $\hat x$, the enclosed flux-surface volume $V$ is normalised to the island width $w$ defined earlier, while in the definition of $\hat \kappa$ the heat conduction coefficient $\kappa (V)$ is normalised to its value at the hot end of the computational domain, which is situated at $\hat x = d$. Additionally, we introduce $x_{w}$ as the position of the wall at the cold end of the domain.

When written in terms of these normalised variables, our heat conduction equation (2.11) and boundary conditions become

(5.3)\begin{gather} \frac{{\rm d}}{{\rm d} x} \left( \kappa(x) \frac{{\rm d}T}{{\rm d} x} \right) = R_0 \left[ \varTheta(T - \tau) - \varTheta(T-1) \right], \end{gather}

(5.4)\begin{gather} T(x_w) = 0, \end{gather}

(5.5)\begin{gather} T'(d) = q_{{\rm in}}, \end{gather}

where $\varTheta$ denotes the Heaviside step function, we have introduced the parameter $\tau ={T_{0}}/{T_{1}}$, and we have dropped carets over all quantities. By integrating this equation once, we find the relations

(5.6a–c)\begin{equation} \frac{\tau}{q_{{\rm out}}} =\int_{x_{w}}^{x_{0}}\frac{1}{\kappa(x)}\,{\rm d}{x} {,}\quad 1= \tau + \int_{x_{0}}^{x_{1}}\frac{R_{0}(x-x_{0})+q_{{\rm out}}}{\kappa(x)}\,{\rm d}{x} {,}\quad T_{{\rm up}} =1 +q_{{\rm in}}\int_{x_{1}}^{d}\frac{1}{\kappa(x)}\,{\rm d}{x}. \end{equation}

It is evident from the dimensionless form of the equation that there are five free parameters, for instance $q_{{\rm in}}$, $\tau$, $d$, $R_0$ and $x_w$. We proceed by computing a number of illustrative solutions where four of these are held fixed while the remaining one is varied. We are particularly interested in seeing how the position and topology of the radiation zone then varies, which can be monitored by plotting $x_0$ and $x_1$, which are defined as the values of $x$ at the cold and hot ends of the radiation belt, respectively, i.e. $T(x_0) = \tau$ and $T(x_1) = 1$. The total radiated power is $q_{{\rm in}} - q_{{\rm out}} = R_0(x_1 - x_0)$ where $q_{{\rm out}}$ represents the heat flux to the wall, and $x_1 - x_0$ is proportional to the radiating volume.

We do not endeavour to survey all possible solutions, but limit our attention to two different parameter scans in which $\tau$ and $d$ are kept constant. In the first case, $q_{{\rm in}}$ is varied while $R_{0}=1$ is kept fixed, whereas in the second case we keep the incoming heat flux fixed and vary $R_{0}$. In addition, in the second case, solutions for different choices of $x_{w}$ are shown as well. To highlight the influence of the island, we also present a solution using the normalised coordinates for a case without an island in the same normalised coordinates outside the island, with a constant heat conduction $\kappa (x)=1$.

Results from the first parameter scan are shown in figure 6, where $x_{0}$ and $x_{1}$ are given as functions of the incoming heat flux $q_{{\rm in}}$. Figure 6(a) shows the position of the radiation zone outside the island, i.e. in the equivalent 1-D geometry between the wall and the hot end, while in figure 6(b) the radiation zone inside the island is shown. It is useful to note that $q_{{\rm in}}$ is not only the heat flux across the upper end of the domain but also represents the temperature gradient between the upper end and the radiation zone at $x_{1}$, where the incoming heat flux starts falling due to the radiated power and the temperature gradient therefore decreases. Thus, if $q_{{\rm in}}$ is large and the radiation zone is situated at the wall, the temperature gradient must be large. The temperature thus changes rapidly with $x$ and reaches the value $\tau$ close to the wall, which is located at $x_w = -4$ in this example. In the limit $q_w \rightarrow \infty$ we thus expect $x_0 \rightarrow x_w$ and $x_1 \rightarrow x_w$, making the radiation band narrow and the radiated power fraction small. In other words, only a small amount of the incoming heat flux is radiated away, $q_{{\rm rad}} / q_{{\rm in}} \rightarrow 0$, and most of the incoming heat flux is deposited on the wall, $q_{{\rm out}} \simeq q_{{\rm in}}$. If $q_{{\rm in}}$ is reduced, the radiation zone recedes from the wall and approaches the island. Below a certain threshold, a multitude of possibilities arise: there is not one, but several solutions to the heat conduction equation, and radiation is possible both inside and outside the island.Footnote ³ Possibilities of radiation outside the island are illustrated in figure 6(a), which shows three pairs of curves. One observes a sudden jump in the sustainable heat influx $q_{{\rm in}}$ (see figure 6a) when radiation inside the island becomes possible. As the width of the radiation zone outside the island remains nearly constant in the vicinity of the island, all the changes in the heat flux are a result from changes inside the island. Only in this close vicinity of the island does the dependence of $x_{1}$ and $x_{0}$ differ substantially from the case without an island. When the radiation zone is close to the O-point, the radiation inside the island remains almost constant. Once the radiation in the island leaves the vicinity of the O-point, the radiation area increases (figure 6b) until the upper end of the radiation zone in the island $x_{1i}$ reaches the separatrix. At this point, the greatest radiation intensity is achieved. Incidentally, in this position, the radiation zone also encloses the area with the greatest change in the enclosed volume, the separatrix. From this point onward, the radiation zone inside the island decreases with decreasing heat flux $q_{{\rm in}}$. When the radiation zone is situated above the island but in close proximity to it, an increase in the radiated heat flux compared with the case without the island is observed. Here, bifurcated solutions are observed, as multiple solutions are possible for a single $q_{{\rm in}}$.

Figure 6. As the incoming heat flux $q_{{\rm in}}$ at $x=d$ is varied, the radiation zone moves around. Multiple positions are possible, both inside and outside the island, and the radiation zone may also straddle the separatrix. The curves show the location of the possible radiation zones (a) outside and (b) inside the island, and figures depict the geometry of the radiation zones in red. The wall is positioned at $x_{w}=-4$ the separatrix at $x_{s}=0.5$, $R_{0}=1$ and $\tau =0.8$. (a) Locations $x_0$ and $x_1$ defining the boundaries of the radiation band outside the island as functions of the incoming heat flux $q_{{\rm in}}$. Also shown as dotted curves are corresponding results for the case without an island. (b) Locations $x_{0i}$ and $x_{1i}$ defining the boundaries of the radiation band inside the island as functions of the incoming heat flux $q_{{\rm in}}$.

It is evident from (3.8) that the rapid change $V'(m)$ in the enclosed volume makes $\kappa$ particularly large close to the separatrix, implying a relatively small temperature gradient in this region and creating a greater capacity to sustain radiation. As mentioned earlier, the island opens up the possibility for multiple positions of the radiation front for the same $q_{{\rm in}}$. Figure 6(a) shows five different possibilities. Some of these may be thermally unstable, but we do not explore this question further.

We now proceed to the second parameter scan, where we keep $q_{{\rm in}}$ fixed but vary the radiation amplitude $R_0$, which mimics the effect of varying the impurity concentration in a fusion device at constant heating power (Stroth et al. Reference Stroth, Bernert, Brida, Cavedon, Dux, Huett, Lunt, Pan and Wischmeier2022). The radiated power fraction $f_{{\rm rad}}={q_{{\rm rad}}}/{q_{{\rm in}}}$ of course increases with increasing $R_0$ but does so in a non-trivial and interesting way. In figure 7, the radiation fraction $f_{{\rm rad}}={q_{{\rm rad}}}/{q_{{\rm in}}}$ is shown as a function of $R_{0}$. To compare $f_{{\rm rad}}$ for different positions of the wall $x_{w}$, $R_{0}$ is normalised to the $R_{0}$ value at $f_{{\rm rad}}=0.4$.

Figure 7. The radiation fraction $f_{{\rm rad}}$ as a function of $R_{0}$. Analytical solutions for a case without an island are also shown as dotted lines for comparison. The influence of position of the wall $x_{w}$ and of the island on the radiation fraction is seen in the different position and gradients of the curves. The plot is normalised to $R_{0}(f_{{\rm rad}}=0.4)$.

The cases without the island show a steady increase in the radiation fraction $f_{{\rm rad}}$ with $R_{0}$, connected to movement of the radiation zone. The cases with an island are very different. Most strikingly, $f_{{\rm rad}}$ is almost constant over a wide range of $R_{0}$. This range corresponds to situations in which radiation occurs from within the island and close to the separatrix. It is thus evident that the presence of an island greatly affects the radiated power. The latter is much larger than in the case without an island, and the radiation has a tendency to linger in the vicinity of the separatrix. For small values of $R_0$, radiation mostly occurs close to the O-point and the radiation zone depends sensitively on $R_{0}$. If $x_{w}$ is decreased, so that the island is brought closer to the wall, the heat flux domain deposited on the wall increases and the range in $R_{0}$ over which $f_{{\rm rad}}$ is constant decreases, as is evident from a comparison of the cases $x_{w}=1.9$ and $x_{w}=2.1$ in figure 7.

Finally, we explore the difference in effects resulting from the non-trivial topology induced by the island and the fact that the 1-D heat conductivity $\kappa (x)$ varies with $x$. To this end, we compare the temperature at the hot boundary of the region in two scenarios: firstly, a case featuring an island, as depicted in figure 6; secondly, a hypothetical scenario with identical topology but with constant heat conductivity $\kappa =1$ throughout the domain (figure 8). This comparison shows that the ability to sustain high radiated heat fluxes in the vicinity of the island can mainly be attributed to the underlying topology: the additional volume made available through the island enhances the dissipation of thermal power. However, it is also evident that the broader range of potential radiation zone positions within the island does not emerge when $\kappa$ is constant, indicating that part of the multiplicity is solely attributable to the spatial variation of $\kappa (x)$. The broader range of solutions stems from this variation within the island. For radiation between the separatrix and the O-point, the higher $\kappa (x)$ near the separatrix enables regions closer to it to radiate more power, as illustrated in figure 8. Similarly, the significantly lower heat conductivity near the O-point, combined with a radiation zone around it, is responsible for the third branch of solutions. This explains the additional branches in $T_{{\rm up}}(\kappa (x))$ shown in figure 8, as compared with the hypothetical scenario with identical topology but constant heat conductivity.

Figure 8. The temperature at the hot boundary of the region as a function of $q_{{\rm in}}$. Shown are the curves for a case featuring an island in black and a hypothetical scenario with identical topology but with constant heat conductivity $\kappa =1$ as a dotted lines for comparison.