1. Necessity of surface currents in stepped-pressure equilibria

In their article ‘Stepped pressure profile equilibria in cylindrical plasmas via partial Taylor relaxation’ (Hole, Hudson & Dewar Reference Hole, Hudson and Dewar2006), Hole, Hudson & Dewar write that their example equilibrium ‘demonstrates the existence of multi-interface, tokamak-like solutions, which do not require the existence of surface currents’. This statement is incorrect, since stepped pressure equilibria necessarily require the existence of surface currents, if Maxwell's equations are to be satisfied. We will demonstrate this in a brief manner, following elementary theory of magnetostatics.

In the multiregion, relaxed magnetohydrodynamic (MRXMHD) model (Hudson et al. Reference Hudson, Dewar, Dennis, Hole, McGann, von Nessi and Lazerson2012) and in the particular case of the cylindrical stepped-pressure equilibrium presented by Hole, Hudson & Dewar and discussed in this comment, pressure jumps are allowed at each ideal interface. At an interface with a pressure jump, force balance then requires

(1.1)\begin{equation} \left\langle p+\frac{B^2}{2\mu_{0}}\right\rangle=0,\end{equation}

where $\mu _{0}$ is the permeability of free space, and $\langle x\rangle = x_{i+1}-x_{i}$ denotes the change in quantity $x$ across the interface $I_{i}$, as the authors write in their (2.3). Now, if there is a pressure jump, force balance (1.1) necessarily also requires a jump in the magnitude $B$ of the magnetic field. One concludes that at a pressure jump, at least one of the components of the magnetic field has a discontinuity. Then, Ampere's law in integral form applied to an Amperian loop straddling the interface immediately implies the existence of a surface current density $\boldsymbol {K}_{i}$, in $\textrm {A} \times \textrm {m}^{-1}$, given by (Griffiths Reference Griffiths2017)

(1.2)\begin{equation} \boldsymbol{n}\times(\boldsymbol{B}_{i+1}-\boldsymbol{B}_{i})= \mu_{0}\boldsymbol{K}_{i},\end{equation}

where $\boldsymbol {n}$ is the unit normal vector to the interface. In the example equilibrium given by Hole, Hudson & Dewar, there are five interfaces at which the pressure has a jump. We conclude that there is a surface current density $\boldsymbol {K}_{i}$ at each of the interfaces, contrary to their claim.

2. Surface currents and continuity of the safety factor

After a careful read of the article, and after analysing figure 1 in it, which shows continuous profiles for the current density, we believe that what the authors mean by the absence of surface currents is the fact that they are able to construct an equilibrium for which, at each plasma–plasma interface, the following is satisfied:

(2.1)\begin{equation} \boldsymbol{J_{i+1}}-\boldsymbol{J_{i}}=\boldsymbol{0}.\end{equation}

To state it with words, the equilibrium they construct is such that at each interface, the current density immediately on the inside of the interface is equal to the current density immediately on the outside of the interface. As we showed in the previous section, this is not equivalent to the absence of surface currents.

Figure 1. Contour plot of $F(\mu _{1},\mu _{2})$ as defined in (3.1) over the domain $(0,160]\times (0,160]$.

Furthermore, even taking into account the misleading definition of a surface current given by the authors, they write the following additional misleading statement in the article: ‘This particular example has been chosen with no change in $q$ across the interfaces, and hence no surface currents’. As we discussed previously, all the interfaces of their equilibrium have surface currents. In a cylinder, no change in $q$ across an interface is a necessary (but insufficient) condition for no surface current: through (1.2) one also has to have continuity of the magnitude $B$ of the magnetic field across the interface. The relationship between $q$, current density and surface currents can be clarified as follows.

The safety factor, $q$, is defined by

(2.2)\begin{equation} q(r)=\frac{2{\rm \pi} r}{L}\frac{B_{z}(r)}{B_{\theta}(r)},\end{equation}

where $r$ is the minor radius of the cylinder, $L$ the length of the (periodic) cylinder, $B_{z}$ is the $z$-component of the magnetic field, $B_{\theta }$ is the $\theta$-component of the magnetic field and $(r,\theta ,z)$ is the cylindrical coordinate system naturally associated with the cylindrical geometry. Let us consider two neighbouring regions $i$ and $i+1$ which are each in a Taylor state with finite current. We therefore have $\boldsymbol {\nabla }\times \boldsymbol {B}_{i}=\mu _{i}\boldsymbol {B}_{i}$ and $\boldsymbol {\nabla }\times \boldsymbol {B}_{i+1}=\mu _{i+1}\boldsymbol {B}_{i+1}$, with $\mu _{i}\neq 0$ and $\mu _{i+1}\neq 0$, so that $\boldsymbol {B}_{i}=({\mu _{0}}/{\mu _{i}})\boldsymbol {J}_{i}$ and $\boldsymbol {B}_{i+1}=({\mu _{0}}/{\mu _{i+1}})\boldsymbol {J}_{i+1}$. Let $r_{i}$ be the radius of the interface between the $i$th region and the $(i+1)$th region. The condition that there be no change in $q$ across the interface is

(2.3)\begin{equation} \frac{B_{z,i}(r_{i})}{B_{\theta,i}(r_{i})}= \frac{B_{z,i+1}(r_{i})}{B_{\theta,i+1}(r_{i})} \Leftrightarrow \frac{J_{z,i}(r_{i})}{J_{\theta,i}(r_{i})}= \frac{J_{z,i+1}(r_{i})}{J_{\theta,i+1}(r_{i})}. \end{equation}

Equation (2.3) shows that the absence of a jump in $q$ across the interface implies the ratio of current densities on either side of the interface is the same, but does not imply the absence of jumps in the current densities at the interfaces. The conclusion of this section is that if one desires to numerically construct stepped-pressure equilibria which do not have jumps in the current densities at the ideal interfaces, solving for equilibria which have a continuous safety factor profile $q$ may not be a satisfactory solution.

3. Impossibility of constructing the desired equilibrium with the data given

In the article, the authors characterize the equilibrium they constructed by providing the radii $r_{i}$ of the interfaces, the coefficients $k_{i}$ and $d_{i}$ in front of the Bessel functions for the magnetic field, and the magnitude of the components $B_{\theta ,V}$ and $B_{z,V}$ of the vacuum field. Since the authors do not provide the values of the Beltrami parameters $\mu _{i}$, the set of parameters prescribing the equilibrium in their (3.4) is not completely specified. We therefore attempted to compute the values of the $\mu _{i}$ in order to fully specify the equilibrium, by enforcing the constraint that there be no jump in the current density at the interfaces, in agreement with figure 1 and with the statements of the authors in the article.

With the information given by the authors, the problem can be solved interface by interface. At the first interface, with radius $r_{1}$, only the Beltrami parameters $\mu _{1}$ and $\mu _{2}$ are unknown. We can therefore look for all the zeros of the function

(3.1)\begin{equation} F(\mu_{1},\mu_{2})=(J_{z,11}(\mu_{1})-J_{z,21}(\mu_{2}))^2 +(J_{\theta,11}(\mu_{1})-J_{\theta,21}(\mu_{2}))^2,\end{equation}

where $J_{z,11}(\mu _{1})=\mu _{1}k_{1}J_{0}(|\mu _{1}|r_{1})$, $J_{\theta ,11}(\mu _{1})=|\mu _{1}|k_{1}J_{1}(|\mu _{1}|r_{1})$, $J_{z,21}(\mu _{2})=\mu _{2}(k_{2}J_{0}(|\mu _{2}| r_{1}) +d_{2}Y_{0}(|\mu _{2}| r_{1}))$, and $J_{\theta ,21}(\mu _{2})=|\mu _{2}|(k_{2}J_{1}(|\mu _{2}|r_{1}) +d_{2}Y_{1}(|\mu _{2}|r_{1}))$, with $k_{1}, k_{2}, d_{2}$ scalar coefficients given in the article, and $J_{0}, J_{1}$ and $Y_{0}, Y_{1}$ Bessel functions of the first kind of order 0, 1 and second kind of order 0, 1, respectively. Note that $F$ is such that $F(-\mu _{1},-\mu _{2})=F(\mu _{1},\mu _{2})$. Thanks to this symmetry with respect to the origin in $(\mu _{1},\mu _{2})$ space, it suffices to look for zeros of $F$ for $\mu _{1}\in \mathbb {R}^{*}, \mu _{2}>0$ (neither $\mu _{1}=0$ nor $\mu _{2}=0$ are allowed since we know from the article that the current densities are finite in regions 1 and 2), and all zeros of $F$ can then be obtained without further computation. Furthermore, the authors write that the obtained equilibrium is tokamak-like. We therefore assume that there is no reversal of the magnetic field. Since the current densities also do not change sign in figure 1 of the article, we can restrict our search to the region in which $\mu _{1}$ and $\mu _{2}$ have the same sign, namely $\mu _{1}>0, \mu _{2}>0$.

In figure 1, we show the contours of $F$ in the domain $(0,160]\times (0,160]$. The sinusoidal nature of the Bessel functions leads to a clear landscape of alternating valleys and ridges. We can search for the global minima of $F$ by looking for the minima in each valley. Doing so, we numerically find two global minima in this domain. To confirm the existence of these two global minima, we use Newton's method to find the zeros of the vector function whose two components are the jump in the $J_{z}$ current density and the jump in the $J_{\theta }$ current density at the interface, and taking the global minima we found as initial guesses for these Newton solves. We indeed find two solutions: $(\mu _{1},\mu _{2})\approx (3.066135, 2.574780)$ and $(\mu _{1},\mu _{2})\approx (141.8329,110.2088)$. We note that by continuing the search for minima inside the valleys beyond the limits of the domain $(0,160]\times (0,160]$, one finds additional global minima. However, they correspond to higher values of both $\mu _{1}$ and $\mu _{2}$. Both the solution $(\mu _{1},\mu _{2})\approx (141.8329,110.2088)$ and these additional minima correspond to highly oscillatory magnetic fields and safety factor, which change sign within the regions in which they are defined. They can be discarded since they do not correspond to the profiles shown in the article. We conclude that there are only two acceptable solutions without a jump in the current densities between region 1 and region 2: $(\mu _{1},\mu _{2})\approx (3.066135, 2.574780)$ and $(\mu _{1},\mu _{2})\approx (-3.066135, -2.574780)$.

For these two values of $\mu _{2}$, we can then look for the values of $\mu _{3}$ such that the current densities do not have a jump at the interface between region 2 and region 3. We define the function

(3.2)\begin{equation} G(\mu_{3})=(J_{z,22}(\mu_{2})-J_{z,32}(\mu_{3}))^2+(J_{\theta,22}(\mu_{2}) -J_{\theta,32}(\mu_{3}))^2,\end{equation}

where $J_{z,22}(\mu _{2})=\mu _{2}(k_{2}J_{0}(|\mu _{2}|r_{2}) +d_{2}Y_{0}(|\mu _{2}|r_{2}))$, $J_{\theta ,22}(\mu _{2})=|\mu _{2}|(k_{2}J_{1}(|\mu _{2}|r_{2}) +d_{2}Y_{1} (|\mu _{2}|r_{2}))$, $J_{z,32}(\mu _{3})=\mu _{3}(k_{3}J_{0}(|\mu _{3}|r_{2}) +d_{3}Y_{0}(|\mu _{3}|r_{2}))$, $J_{\theta ,32}(\mu _{3})=|\mu _{3}|(k_{3}J_{1}(|\mu _{3}|r_{2}) +d_{3}Y_{1}(|\mu _{3}|r_{2}))$, with $r_{2}$ the radius of the interface between region 2 and region 3, and $k_{3}, d_{3}$ scalar coefficients given in the article. One can show that for $\mu _{2}\approx 2.574780, G$ has a unique global minimum for $\mu _{3}\in (0,160]$, with approximate value $1.843450 \times 10^{-3}$, reached for $\mu _{3}\approx 2.176953$. Likewise, for $\mu _{2}\approx -2.574780, G$ has a unique global minimum for $\mu _{3}\in [-160,0)$, with approximate value $1.843450\times 10^{-3}$, reached for $\mu _{3}\approx -2.176953$. For these values of $\mu _{2}$ and $\mu _{3}$, the magnitude of the jump in $J_{\theta }$ at the interface is approximately $2.19350\times 10^{-2}$ and the magnitude of the jump in $J_{z}$ at the interface is approximately $3.69101\times 10^{-2}$; it is finite.

We conclude that with the parameter values given by the authors in the article, it is not possible to construct a stepped-pressure equilibrium such that the jump in the current densities is zero at each interface within the plasma. We invite the authors to provide in a corrigendum the correct values for the coefficients $k_{i}, d_{i}$ which make the equilibrium shown in figure 1 of the article realizable, and also provide the values of the Beltrami parameters $\mu _{i}$, in order to completely specify the equilibrium.

In closing, we would like to emphasize the fact that the authors of the article we comment upon are experts of MRXMHD and stepped-pressure equilibria, have published a large number of excellent articles on the topic and, as far as we know, have not repeated their incorrect statements in these other articles. We therefore do not think we are addressing a controversial question in this comment on their article. Still, given that the article has gathered a fair number of citations, indicating a fair number of reads, we thought our comment could save time for future readers, who otherwise may ponder these questions just like we have for a little while.

For that same reason, we also highlight the following typographical errors in the manuscript. In (3.1) of the article, the magnetic field should be

(3.3)\begin{equation} \boldsymbol{B}=\{0,\mbox{sign}(\mu_{1})k_{1}J_{1}(|\mu_{1}|r),k_{1}J_{0}(|\mu_{1}|r)\}. \end{equation}

In (3.2) of the article, the magnetic field should be

(3.4)\begin{equation} \boldsymbol{B}=\{0,\mbox{sign}(\mu_{i})(k_{i}J_{1}(|\mu_{i}|r) +d_{i}Y_{1}(|\mu_{i}|r),k_{i}J_{0}(|\mu_{i}|r)+d_{i}Y_{0}(|\mu_{i}|r)\}. \end{equation}

The absolute value of the Beltrami parameter inside the Bessel functions is required because when $u$ is a negative real number, $Y_{0}(u)$ and $Y_{1}(u)$ are in general complex numbers, with a non-zero imaginary part. This would correspond to components of the magnetic field which have a non-zero imaginary part, which is not physical. Mathematically, the absolute value can be introduced inside the Bessel functions regardless of the sign of the Beltrami parameters without restricting the solution space because only the squares of the Beltrami parameters, $\mu _{i}^2$, appear in the differential equations for $B_{\theta }$ and $B_{z}$. Strictly speaking, the absolute value is only required for (3.2). It is not required in (3.1) (provided one removes $\mbox {sign}(\mu _{1})$ at the same time as the absolute value) since $J_{0}(u)$ and $J_{1}(u)$ are real numbers for $u\in \mathbb {R}$, and $J_{0}$ is an even function and $J_{1}$ an odd function, as desired. However, for the sake of a consistent notation, we believe it is good to use the expressions with absolute values for (3.1) as well. Finally, the equilibrium must be specified with $4N+1$ parameters, instead of $4N+2$ parameters as stated in the article.