2.1. Description of toroidal volumes and surfaces

Throughout our discussion of toroidal geometries, we make use of the standard, right-handed cylindrical coordinates $(r, \phi, z)$. At a point with toroidal angle $\phi$, we write the orthonormal unit vectors as $\boldsymbol {e}_r(\phi )$, $\boldsymbol {e}_\phi (\phi )$ and $\boldsymbol {e}_z$. With this notation, we emphasise the fact that the radial and azimuthal unit vectors depend on the toroidal angle.

In this article, we focus on axisymmetric geometries, which means that we only consider surfaces and volumes of revolution. We take the $z$-axis as the axis of revolution and define a simple closed curve $\gamma$ in the $(r,z)$ plane. By rotating this curve about the $z$-axis through the toroidal angle $\phi \in [0, 2{\rm \pi} ]$, we obtain a closed surface of revolution $\varGamma$. Its interior $\varOmega$ is the corresponding volume of revolution. We refer to $\gamma$ as the generating curve of $\varGamma$. It is parameterised by a single variable $t$, which we assume has period $L$. We denote the components of $\gamma$ in the $(r,z)$ plane by $(r(t), z(t))$, and we identify a point $\boldsymbol {y} \in \varGamma$ by its toroidal revolution angle $\phi$ and its generating curve parameter $t$. Correspondingly, we often write $\boldsymbol {y} = \boldsymbol {y}(\phi, t)$ to stress this parameterisation. Moreover, we assume that $\gamma$ is a $C^1$ curve which does not intersect the $z$-axis, in the sense that the derivatives $r'(t)$ and $z'(t)$ are continuous on $[0,L]$ and there exists $R_{\rm min}>0$ for which $r(t) \geqslant R_{\rm min}$ on $[0,L]$. Finally, we assume that $\gamma$ is oriented so that the vector $\boldsymbol {n}(\boldsymbol {y}(\phi, t)) = (\partial \boldsymbol {y} / \partial \phi ) \times (\partial \boldsymbol {y} / \partial t) / \mathcal {J}(t)$ is the unit outward normal to $\varOmega$ at $\boldsymbol {y}$. The quantity $\mathcal {J}(t)=\left \| (\partial \boldsymbol {y} / \partial \phi ) \times (\partial \boldsymbol {y} / \partial t) \right \|$ is the Jacobian of the parameterisation.

2.2. Single-layer and double-layer potentials for axisymmetric geometries

Layer potentials are fundamental tools in representing solutions to the partial differential equations that arise in magnetostatic and magnetohydrodynamic calculations for magnetic confinement fusion (Merkel Reference Merkel1986; Chance Reference Chance1997; Ludwig et al. Reference Ludwig, Bosco, Ferreira and Berni2006, Reference Ludwig, Rodrigues and Bizarro2013; Landreman & Boozer Reference Landreman and Boozer2016; Drevlak et al. Reference Drevlak, Beidler, Geiger, Helander and Turkin2018). Given a surface $\varGamma$ and a free-space Green's function $(\boldsymbol {x}, \boldsymbol {y}) \mapsto G(\boldsymbol {x}, \boldsymbol {y})$ for a partial differential equation, the single-layer operator $\mathcal {S}$ and the double-layer operator $\mathcal {D}$ are defined by (Guenther & Lee Reference Guenther and Lee1996)

(2.1a,b)\begin{align} [\mathcal{S} \sigma](\boldsymbol{x}) = \iint_\varGamma G(\boldsymbol{x}, \boldsymbol{y}) \sigma(\boldsymbol{y}) \,\mathrm{d} \varGamma(\boldsymbol{y}) \quad \text{and} \quad [\mathcal{D} \sigma](\boldsymbol{x}) = \iint_\varGamma \frac{\partial G(\boldsymbol{x}, \boldsymbol{y})}{\partial \boldsymbol{n}(\boldsymbol{y})} \sigma(\boldsymbol{y}) \,\mathrm{d} \varGamma(\boldsymbol{y}), \end{align}

respectively. In the double-layer representation, the quantity $\partial G(\boldsymbol {x}, \boldsymbol {y}) / \partial \boldsymbol {n}(\boldsymbol {y}) = \boldsymbol {n}(\boldsymbol {y}) \boldsymbol {\cdot } \boldsymbol \nabla _{\boldsymbol {y}} G(\boldsymbol {x}, \boldsymbol {y})$ is the derivative of the Green's function in the outward normal direction at $\boldsymbol {y}$. The function $\sigma$ is called the density function in this representation. Functions expressible as single-layer and double-layer potentials automatically satisfy the partial differential equation associated with the Green's function everywhere except the boundary $\varGamma$.

The case of Laplace's equation in three dimensions is particularly prevalent in magnetostatic and magnetohydrodynamic settings (Merkel Reference Merkel1986; Chance Reference Chance1997; Landreman & Boozer Reference Landreman and Boozer2016; Drevlak et al. Reference Drevlak, Beidler, Geiger, Helander and Turkin2018). Here, the Green's function is

(2.2)\begin{equation} G(\boldsymbol{x}, \boldsymbol{y}) = \frac 1{4 {\rm \pi}\left\| \boldsymbol{x} - \boldsymbol{y} \right\|}, \end{equation}

and functions expressible as $\mathcal {S} \sigma$ or $\mathcal {D} \sigma$ are harmonic on $\mathbb {R}^3 \setminus \varGamma$. Assuming the density $\sigma$ is also axisymmetric in the sense that $\partial \sigma / \partial \phi = 0$, one can analytically compute the part of the surface integral over the revolution angle $\phi \in [0, 2{\rm \pi} ]$. The resulting single-layer and double-layer integrals then are one-dimensional line integrals, and the resulting Green's function in the integrand can be expressed in terms of complete elliptic integrals (Ludwig et al. Reference Ludwig, Bosco, Ferreira and Berni2006, Reference Ludwig, Rodrigues and Bizarro2013; Jardin Reference Jardin2010). We provide an explicit expression in § 4.1, as we treat in detail the application of our method to the calculation of the virtual-casing principle. At this point, we just highlight the fact that the singularity in the Green's function when $\boldsymbol {x}=\boldsymbol {y}$ requires the use of specialised quadrature when the target $\boldsymbol {x}$ is located on the surface $\varGamma$, or regularisation methods (Freidberg et al. Reference Freidberg, Grossmann and Haas1976; Merkel Reference Merkel1986; Chance Reference Chance1997; Landreman & Boozer Reference Landreman and Boozer2016; Drevlak et al. Reference Drevlak, Beidler, Geiger, Helander and Turkin2018; Malhotra et al. Reference Malhotra, Cerfon, O'Neil and Toler2019b). The purpose of this article is to show that for applications in axisymmetric geometries, the Kapur–Rokhlin quadrature scheme is simpler to implement than the known regularisation methods used in the magnetic confinement community, and leads to high-order convergence.

2.3. Kapur–Rokhlin quadrature

The Kapur–Rokhlin quadrature rules (Kapur & Rokhlin Reference Kapur and Rokhlin1997) are a collection of high-order schemes for computing

(2.3)\begin{equation} \int_0^b f(t) \,\mathrm{d} t \end{equation}

when $f$ has an integrable singularity at the origin of the form $\log t$ or $t^\lambda$ for $\lambda >-1$. We focus on the specific Kapur–Rokhlin scheme for a logarithmic singularity of the form $f(t) = p(t) \log t + q(t)$, where $p$ and $q$ need not have known formulae, but are assumed sufficiently smooth.

2.3.1. A Kapur–Rokhlin quadrature scheme for non-periodic functions

The Kapur–Rokhlin quadrature rules are corrections to the trapezoidal rule. For the standard trapezoidal rule with equal spacing $h = b/M$, the quadrature nodes for a non-singular integrand would be $t_i = i h$ for $i=0, \dots, M$. However, we omit the quadrature node $t_0=0$ because the integrand $f$ is singular there. This yields the punctured trapezoidal rule

(2.4)\begin{equation} \int_0^b f(t) \,\mathrm{d} t \approx h \left[ \sum_{i=1}^{M-1} f_i + \frac 12 f_M \right], \end{equation}

where we have written the shorthand $f_i = f(t_i)$. The punctured trapezoidal rule is low-order accurate when $f$ is singular. For example, when $f$ is logarithmically singular at the origin, the punctured trapezoidal rule error typically decays, according to Martinsson (Reference Martinsson2019), as $O(|h\log h|)$.

The Kapur–Rokhlin corrections place additional quadrature nodes outside the integration domain $[0,b]$. Specifically, the corrections depend on a convergence rate parameter $n$ and a smoothness parameter $m$. One chooses $m \geqslant 3$ an odd integer subject to the constraint that $p$ and $q$ are $m$ times continuously differentiable on a wider interval $[-nh,b+mh]$. Under these conditions, the Kapur–Rokhlin scheme sets constants $\gamma _j$ for $j=\pm 1, \dots, \pm n$ and $\beta _l$ for $l = 1, \dots, (m-1)/2$ and defines the quadrature rule

(2.5)\begin{equation} T_{m,n}^M(f) = h \left[\sum_{i=1}^{M-1} f_i + \frac 12 f_M \right] + h \sum_{l=1}^{(m-1)/2} \beta_l (f_{M-l} - f_{M+l}) + h \sum_{1 \leqslant |j| \leqslant n} \gamma_j f_j \end{equation}

for $M \geqslant n + (m-1)/2$. The error obeys the asymptotic rate

(2.6)\begin{equation} \left| T_{m,n}^M(f) - \int_0^b f(t) \,\mathrm{d} t \right| = O(h^n) \end{equation}

as $h \to 0$. The endpoint correction coefficients $\{\gamma _j\}$ and $\{\beta _l\}$ are derived by analysing the Euler–Maclaurin formula for quadrature error and solving a linear system for correction coefficients to obtain high-order accuracy. Figure 1 illustrates this quadrature scheme and compares it with the punctured trapezoidal rule.

Figure 1. Nodes and weights for the Kapur–Rokhlin and punctured trapezoidal quadrature rules for estimating $\int _0^1 f(t) \,\mathrm {d} t$ with $f(t) = \cos (4{\rm \pi} t) \log |t| + t$. We have used $M=10$, $n=2$, and $m=9$. Note that $\beta _4 \approx -3 \times 10^{-4}$, so the weight corrections corresponding to $t_{6}=0.6$ and $t_{14}=1.4$ are not visually discernible.

2.3.2. A simplified quadrature for periodic functions

This subsection follows an argument nearly verbatim from Hao et al. (Reference Hao, Barnett, Martinsson and Young2014). We explain how the Kapur–Rokhlin scheme we just presented simplifies when computing

(2.7)\begin{equation} \int_{{-}b}^b f(t)\, \mathrm{d} t \end{equation}

when $f$ is $2b$-periodic and logarithmically singular at the origin. We may express these assumed properties of $f$ through the form

(2.8)\begin{equation} f(t) = p(t) \log \left|\sin \frac{{\rm \pi} t}{2b} \right| + q(t), \end{equation}

for $2b$-periodic functions $p$ and $q$.

We sum two applications of the original Kapur–Rokhlin scheme: one for

(2.9)\begin{equation} I_1 = \int_0^b f(t) \,\mathrm{d} t \end{equation}

and another for

(2.10)\begin{equation} I_2 = \int_{{-}b}^0 f(t) \,\mathrm{d} t = \int_0^b f({-}t) \,\mathrm{d} t. \end{equation}

We assume that $p, q \in C^m[-b, b]$ and generate $2M-1$ equispaced quadrature nodes with spacing $h=b/M$, defined by $t_i = ih$ for $i=\pm 1, \dots, \pm (M-1), M$. The corrected scheme for $I_1$ is

(2.11)\begin{equation} I_1 = h \left[\sum_{i=1}^{M-1} f_i + \frac 12 f_M \right] + h \sum_{l=1}^{(m-1)/2} \beta_l (f_{M-l} - f_{M+l}) + h \sum_{1 \leqslant |j| \leqslant n} \gamma_j f_j + O(h^n), \end{equation}

and the scheme for $I_2$ is

(2.12)\begin{equation} I_2 = h \left[\sum_{i=1}^{M-1} f_{{-}i} + \frac 12 f_{{-}M} \right] + h \sum_{l=1}^{(m-1)/2} \beta_l (f_{{-}M+l} - f_{{-}M-l}) + h \sum_{1 \leqslant |j| \leqslant n} \gamma_j f_{{-}j} + O(h^n). \end{equation}

By periodicity, we may identify $f_i$ with $f_{i+2M}$ for all $i$. It follows that an $n$th-order quadrature for $I_1 + I_2$ is

(2.13)\begin{align} I_1 + I_2 & = h \left[ \sum_{1 \leqslant |i| \leqslant M-1} f_i + f_M \right] + h \sum_{1 \leqslant |j| \leqslant n} \gamma_j (f_j + f_{{-}j})+ O(h^n) \nonumber\\ & = \sum_{\substack{j={-}M+1\nonumber\\ j \ne 0}}^M w_j f_j + O(h^n) \end{align}

with

(2.14)\begin{equation} w_j = \begin{cases} h(1+\gamma_j + \gamma_{{-}j}), & 1 \leqslant |j| \leqslant n, \\ h, & \text{otherwise}. \end{cases} \end{equation}

We note that, by periodicity, the endpoint corrections corresponding to the constants $\beta _l$ exactly cancel. Moreover, we no longer need to place additional quadrature nodes beyond $[-b,b]$ because periodicity identifies the new quadrature nodes with existing nodes.

Given a table of $\gamma _j$ weights, this quadrature scheme is as easy to implement as the trapezoidal rule and yields high-order convergence, with a quadrature error that is $O(h^n)$. The necessary $\gamma _j$ weights can be found in Kapur & Rokhlin (Reference Kapur and Rokhlin1997, table 6) for $n=2,6,10$. Figure 2 can be compared with figure 1 to view the simplifications for periodic integrands.

Figure 2. Nodes and weights for the Kapur–Rokhlin and punctured trapezoidal quadrature rules for estimating $\int _{t_0-{\rm \pi} }^{t_0+{\rm \pi} } f(t) \,\mathrm {d} t$. The function $f(t)$ is the $2{\rm \pi}$-periodic integrand in (6.5) for one of our later numerical tests with a logarithmic singularity at $t_0=1$. We have used $M=10$ and $n=2$.

3.1. Formulation of the virtual-casing principle

For axisymmetric confinement devices, the virtual-casing principle is most often used to compute the poloidal flux or the poloidal magnetic field due to the toroidal current flowing in the plasma (Shafranov & Zakharov Reference Shafranov and Zakharov1972; Zakharov Reference Zakharov1973; Hirshman & Neilson Reference Hirshman and Neilson1986; Zakharov & Pletzer Reference Zakharov and Pletzer1999). A poloidal magnetic field ${\boldsymbol {B}}^{{\rm pol}}$ at any point $\boldsymbol {y}(\phi, t) \in \varGamma$ can be expressed in terms of its poloidal flux function $\psi (r,z)$ and the parameterisation $(\phi,t) \mapsto \boldsymbol {y}(\phi, t)$ by (Freidberg Reference Freidberg2014)

(3.1)\begin{align} \boldsymbol{B}^{{\rm pol}}(\boldsymbol{y}(\phi,t)) & = \boldsymbol\nabla \psi(r(t),z(t)) \times \boldsymbol\nabla \phi \nonumber\\ & = \boldsymbol\nabla \psi(r(t),z(t)) \times \left(\frac {\boldsymbol{e}_\phi(\phi)}{r(t)} \right) \nonumber\\ & ={-}\frac 1{r(t)} \dfrac{\partial \psi}{\partial z}(r(t), z(t)) \boldsymbol{e}_r(\phi) + \frac 1{r(t)} \dfrac{\partial \psi}{\partial r}(r(t), z(t)) \boldsymbol{e}_z. \end{align}

Consider an axisymmetric plasma confined by external coils in equilibrium. The poloidal field $\boldsymbol {B}^{{\rm pol}}$ at any location is the sum of the poloidal field $\boldsymbol {B}_{{\rm ext}}^{{\rm pol}}$ due to the external coils and of the poloidal field $\boldsymbol {B}_V^{{\rm pol}}$ due to the plasma current. The field $\boldsymbol {B}_V^{{\rm pol}}$ is given for all $\boldsymbol {x} \in \mathbb {R}^3$ by the Biot–Savart law:

(3.2)\begin{equation} \boldsymbol{B}_V^{\rm pol}(\boldsymbol{x})=\frac{\mu_{0}}{4{\rm \pi}} \iiint_{\varOmega}J^{{\rm tor}}_{V}(\boldsymbol{y}) \boldsymbol{e}_\phi(\boldsymbol{y}) \times \frac{\boldsymbol{x} - \boldsymbol{y}}{\left\| \boldsymbol{x} - \boldsymbol{y} \right\|^3} \,\mathrm{d} \boldsymbol{y}, \end{equation}

where $\mu _{0}$ is the permeability of free space and $J^{{\rm tor}}_{V}$ is the toroidal current density in the plasma. Equation (3.2) is a volume integral, which is expensive to evaluate numerically. The virtual-casing principle gives a formula for $\boldsymbol {B}^{{\rm pol}}_V$ that depends only on the full field $\boldsymbol {B}^{{\rm pol}}$ at the plasma boundary, and only requires the evaluation of a surface integral (i.e. line integral for axisymmetric domains) (Shafranov & Zakharov Reference Shafranov and Zakharov1972; Zakharov Reference Zakharov1973; Hanson Reference Hanson2015). In experimental settings, such a representation is useful because only the total magnetic field $\boldsymbol {B}^{{\rm pol}}$ may be directly measurable (Hutchinson Reference Hutchinson2002). In theoretical settings, the total magnetic field may be computed by solving the Grad–Shafranov equation (Grad & Rubin Reference Grad and Rubin1958; Shafranov Reference Shafranov1958) numerically with a fixed-boundary solver. Specifically, the virtual-casing principle states that if $\varGamma$ is the flux surface bounding the plasma, then $\boldsymbol {B}_V^{{\rm pol}}$ can be written in terms of a field generated by the toroidal surface current $\boldsymbol {J}^{{\rm tor}}_S$ such that $\mu _0 \boldsymbol {J}^{{\rm tor}}_S = -\boldsymbol {n} \times \boldsymbol {B}^{{\rm pol}}$, according to Hanson (Reference Hanson2015):

(3.3)\begin{equation} \boldsymbol{B}^{{\rm pol}}_V(\boldsymbol{x}) = \frac 1{4{\rm \pi}} \iint_{\varGamma} \left[ \frac{(\boldsymbol{n}(\boldsymbol{y}) \times \boldsymbol{B}^{{\rm pol}}(\boldsymbol{y})) \times (\boldsymbol{x} - \boldsymbol{y})}{\left\| \boldsymbol{x} - \boldsymbol{y} \right\|^3} \right] \mathrm{d} \varGamma(\boldsymbol{y}) + \begin{cases} \boldsymbol{B}^{{\rm pol}}(\boldsymbol{x}) & \boldsymbol{x} \in \varOmega \\ \boldsymbol{B}^{{\rm pol}}(\boldsymbol{x}) / 2 & \boldsymbol{x} \in \varGamma \\ \boldsymbol{0} & \boldsymbol{x} \notin \bar{\varOmega}. \end{cases} \end{equation}

For certain applications, one is only interested in the normal component of the poloidal magnetic field (Merkel Reference Merkel1986; Hirshman et al. Reference Hirshman, van Rij and Merkel1986; Merkel Reference Merkel1987; Landreman Reference Landreman2017; Zhu et al. Reference Zhu, Hudson, Song and Wan2018). The previous equation then leads to the more compact form

(3.4)\begin{equation} \boldsymbol{n}(\boldsymbol{x}) \boldsymbol{\cdot} \boldsymbol{B}^{{\rm pol}}_V(\boldsymbol{x}) = \frac 1{4{\rm \pi}} \boldsymbol{n}(\boldsymbol{x}) \boldsymbol{\cdot} \iint_{\varGamma} \left[ \frac{(\boldsymbol{n}(\boldsymbol{y}) \times \boldsymbol{B}^{{\rm pol}}(\boldsymbol{y})) \times (\boldsymbol{x} - \boldsymbol{y})}{\left\| \boldsymbol{x} - \boldsymbol{y} \right\|^3} \right] \mathrm{d} \varGamma(\boldsymbol{y}) \end{equation}

for $\boldsymbol {x} \in \varGamma$.

The reduction of the integral necessary to compute the field or its normal component from a volume integral to a surface integral is convenient from the point of view of the limited number of values that need to be specified as inputs, and also from the point of view of the computational cost of the integration (Lazerson et al. Reference Lazerson, Sakakibara and Suzuki2013). The surface integral in (3.3) is significantly faster to evaluate than the volume integral (3.2), although certain codes still choose to compute the latter (Hanson et al. Reference Hanson, Hirshman, Knowlton, Lao, Lazarus and Shields2009; Marx & Lütjens Reference Marx and Lütjens2017).

For axisymmetric situations, one may further take advantage of the axisymmetry of $\boldsymbol {B}^{{\rm pol}}$ to integrate with respect to $\phi$ analytically, and reduce (3.3) and (3.4) to one-dimensional integrals, which are even less computationally expensive. However, one encounters a mathematical and computational difficulty if one does so, because the surface integrals in (3.3) and (3.4) are in fact improper integrals, which must be understood in the Cauchy principal value sense. This is what we discuss in the following section.

3.2. Numerical evaluation of the normal component of the virtual-casing magnetic field in axisymmetric systems

3.2.1. Circumventing integrals in the principal value sense

Most applications in magnetic confinement fusion rely on version (3.4) of the virtual-casing principle, in which one wants to compute the normal component of the poloidal magnetic field $\boldsymbol {B}_V^{{\rm pol}}$ at a boundary point $\boldsymbol {x}\in \varGamma$. One could, in principle, calculate this normal component by first computing all the components of $\boldsymbol {B}_V^{{\rm pol}}$ by the virtual-casing principle (3.3), and then computing $\boldsymbol {n} \boldsymbol {\cdot } \boldsymbol {B}_V^{{\rm pol}}$ by a straightforward inner product. In other words, one first evaluates the double integral in (3.4), and then takes the dot product with the normal vector $\boldsymbol {n}$ to the surface $\varGamma$ at the point $\boldsymbol {x}$ of interest. This method has the advantage that it produces all components of the poloidal magnetic field $\boldsymbol {B}_V^{{\rm pol}}$ as intermediary results. However, it also has the disadvantage that one must use a careful principal value integration procedure to interpret the virtual-casing principle for $\boldsymbol {x} \in \varGamma$, as discussed in Theorem A.3 and its proof in the Appendix. The high-order singularity cancellation quadrature scheme we recently proposed for singular integrals on general non-axisymmetric surfaces (Malhotra et al. Reference Malhotra, Cerfon, O'Neil and Toler2019b) automatically yields the appropriate principal value of the integral. However, it does so thanks to the intrinsic two-dimensional nature of the integral. It does not reduce to a simple and efficient one-dimensional quadrature scheme for the principal value of the integral, as is needed for axisymmetric applications. Similarly, we are not aware of a version of the Kapur–Rokhlin quadrature scheme designed to calculate the Cauchy principal value of the virtual-casing integral.

To address this difficulty, we consider an alternative method to compute the normal poloidal field $\boldsymbol {n} \boldsymbol {\cdot } \boldsymbol {B}_V^{{\rm pol}}$, based on calculating the vector potential $\boldsymbol {A}_S$ produced by the surface current $\boldsymbol {J}^{{\rm tor}}_S$, and then obtaining $\boldsymbol {n} \boldsymbol {\cdot } \boldsymbol {B}_V^{\rm pol}$ as the tangential derivative of the poloidal flux $\psi _S$ caused by the surface currents, which is easily expressed in terms of $\boldsymbol {A}_S$. This method is also used, in a slightly different form, by the stellarator optimisation code ROSE (Drevlak et al. Reference Drevlak, Beidler, Geiger, Helander and Turkin2018), but ROSE does not combine it with a high-order quadrature scheme. The vector potential is defined by

(3.5)\begin{equation} \boldsymbol{A}_S(\boldsymbol{x}) = \frac {\mu_0}{4 {\rm \pi}} \iint_\varGamma \frac{\boldsymbol{J}^{{\rm tor}}_S(\boldsymbol{y})}{\left\| \boldsymbol{x} - \boldsymbol{y} \right\|} \,\mathrm{d} \varGamma(\boldsymbol{y}) \end{equation}

and we recognise this expression as a vector of component-wise single-layer potentials for Laplace's equation in three dimensions, as introduced in § 2.2.

In contrast to the first approach, our method does not produce all components of $\boldsymbol {B}^{{\rm pol}}_V$. On the other hand, it only requires weakly singular integrals, because it only requires evaluations of single-layer potentials in (3.5), as we prove in § 4. In the next section, we describe this alternative method for computing the normal magnetic field in detail.

3.2.2. Normal component of the magnetic field as a tangential derivative

As in the derivation of the virtual-casing principle (Shafranov & Zakharov Reference Shafranov and Zakharov1972; Zakharov Reference Zakharov1973; Hanson Reference Hanson2015), one can interpret the plasma boundary $\varGamma$ of the axisymmetric plasma as a perfectly conducting shell which contributes to confining the plasma via the poloidal magnetic field $\boldsymbol {B}^{{\rm pol}}_S$ generated by the toroidal surface current density $\boldsymbol {J}^{{\rm tor}}_S$ flowing in $\varGamma$. At equilibrium, the poloidal surface current field $\boldsymbol {B}^{{\rm pol}}_S$ matches the poloidal field $\boldsymbol {B}_{{\rm ext}}^{{\rm pol}}$ from the external confining coils. The decomposition $\boldsymbol {B}^{{\rm pol}} = \boldsymbol {B}^{{\rm pol}}_V + \boldsymbol {B}^{{\rm pol}}_S$ then immediately yields $\boldsymbol {n} \boldsymbol {\cdot } \boldsymbol {B}^{{\rm pol}}_V = -\boldsymbol {n} \boldsymbol {\cdot } \boldsymbol {B}^{{\rm pol}}_S$. Now, recall that we may represent $\boldsymbol {B}^{\rm pol}_S$ in axisymmetry via (3.1) to express its normal component at $(r, \phi, z)$ as

(3.6)\begin{equation} \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{B}^{{\rm pol}}_S = \boldsymbol{n} \boldsymbol{\cdot} (\boldsymbol\nabla \psi_S \times \boldsymbol\nabla \phi) ={-}\frac{1}{r} (\boldsymbol{n} \times \boldsymbol{e}_\phi(\phi)) \boldsymbol{\cdot} \boldsymbol\nabla \psi_S \end{equation}

by the circular shift identity of the vector triple product. Let $\boldsymbol {x}$ correspond to the parameter pair $(\phi _0, t_0)$ with $\boldsymbol {x} = (R, \phi _0, Z) = (r(t_0), \phi _0, z(t_0))$. We define the tangent vector $\boldsymbol {t}$ at a point $\boldsymbol {x} \in \varGamma$ as

(3.7)\begin{align} \boldsymbol{t}(\boldsymbol{x}) & = \boldsymbol{n}(\boldsymbol{x}) \times \boldsymbol{e}_\phi(\phi_0) \nonumber\\ & = \frac{(z'(t_0) \boldsymbol{e}_r(\phi_0) - r'(t_0) \boldsymbol{e}_z ) \times \boldsymbol{e}_\phi(\phi_0)} {\sqrt{r'(t_0)^2 + z'(t_0)^2}} \nonumber\\ & = \frac{ r'(t_0) \boldsymbol{e}_r(\phi_0) + z'(t_0) \boldsymbol{e}_z } {\sqrt{r'(t_0)^2 + z'(t_0)^2}}. \end{align}

It follows that

(3.8)\begin{align} \boldsymbol{n}(\boldsymbol{x}) \boldsymbol{\cdot} \boldsymbol{B}^{\rm pol}_V(\boldsymbol{x}) & = \frac{1}{r(t_0)} \boldsymbol{t}(\boldsymbol{x}) \boldsymbol{\cdot} \boldsymbol\nabla \psi_S(r(t_0), z(t_0)) \nonumber\\ & = \frac{1}{r(t_0)\sqrt{r'(t_0)^2 + z'(t_0)^2}} \left( \frac{\partial{\psi_S}}{\partial{r}} r'(t_0) + \frac{\partial{\psi_S}}{\partial{z}} z'(t_0) \right)\nonumber\\ & = \frac{1}{\mathcal{J}(t_0)} \frac{\mathrm{d} \psi_S(r(t_0), z(t_0))}{\mathrm{d} t_0}, \end{align}

where $\mathcal {J}(t_0) = r(t_0)\sqrt {r'(t_0)^2 + z'(t_0)^2}$ is the Jacobian of the parameterisation introduced in § 2.1.

Finally, the poloidal flux function and the vector potential are related at $\boldsymbol {x}$ by the relation (Jardin Reference Jardin2010; Freidberg Reference Freidberg2014)

(3.9)\begin{equation} \psi_S(R, Z) = R \boldsymbol{e}_\phi(\phi_0) \boldsymbol{\cdot} \boldsymbol{A}_S(\boldsymbol{x}). \end{equation}

Given only point evaluations of $\boldsymbol {A}_S$ at equispaced parameters $\{t_i\}$, one can compute $\mathrm {d} \psi _S / \mathrm {d} t$ with high-order accuracy at each $t_i$ by Fourier differentiation. This leads to our high-order accurate approach to the virtual-casing principle in axisymmetric geometries. One computes $\boldsymbol {n} \boldsymbol {\cdot } \boldsymbol {B}^{{\rm pol}}_V$ on $\varGamma$ by evaluating a weakly singular integral with a high-order accurate quadrature rule, and then Fourier differentiating. We show in § 5 that Kapur–Rokhlin quadrature provides a simple way to obtain high-order accuracy for the weakly singular integral. Before we do so, we prove in the next section that it indeed is a weakly singular integral, whose properties satisfy the requirements to obtain high accuracy with the Kapur–Rokhlin quadrature rule.

4.1. Vector potential singularity

By the virtual-casing principle, we have $\mu _0 \boldsymbol {J}^{{\rm tor}}_S = - \boldsymbol {n} \times \boldsymbol {B}^{{\rm pol}}$, so we may equivalently write the vector potential as

(4.1)\begin{equation} \boldsymbol{A}_S(\boldsymbol{x}) ={-}\frac{1}{4{\rm \pi}} \iint_\varGamma \frac{\boldsymbol{n}(\boldsymbol{y}) \times \boldsymbol{B}^{{\rm pol}}(\boldsymbol{y})} {\left\| \boldsymbol{x} - \boldsymbol{y} \right\|} \mathrm{d} \varGamma(\boldsymbol{y}), \end{equation}

which satisfies $\boldsymbol {B}^{{\rm pol}}_S = \boldsymbol { \nabla } \times \boldsymbol {A}_S$. As mentioned earlier, we may view this expression as a vector of component-wise single-layer potentials for Laplace's equation in three dimensions. Physically, it is clear that the only non-zero component of the vector potential is in the $\boldsymbol {e}_\phi$ direction. Furthermore, Chance (Reference Chance1997, § V) shows that after integrating in the toroidal angle $\phi$, the one-dimensional integrands of both the single-layer potential $\mathcal {S}[\mu _0\boldsymbol {J}^{{\rm tor}}_S]$ and the double-layer potential $\mathcal {D}[\mu _0\boldsymbol {J}^{{\rm tor}}_S]$ have integrable, logarithmic singularities. Chance shows this in a modified coordinate system using the poloidal flux function as a coordinate. Next, we verify these results for the virtual-casing principle in standard cylindrical coordinates.

4.2. Analytic reduction to a line integral

We analytically simplify the integral in (4.1) by integrating over the toroidal angle. When we do so, we introduce the complete elliptic integrals of the first and second kind, which are defined as

(4.2a,b)\begin{equation} K(k^2) = \int_0^{{\rm \pi}/2} \frac{\mathrm{d} \theta}{\sqrt{1 - k^2 \sin^2\theta}} \quad \text{and} \quad E(k^2) = \int_0^{{\rm \pi}/2} \sqrt{1 - k^2 \sin^2\theta} \,\mathrm{d} \theta, \end{equation}

respectively. Our main analytical result gives a formula for the vector potential of the virtual-casing principle as a one-dimensional integral.

Theorem 4.1 Let $\varGamma$ be a surface of revolution with a generating curve $\gamma$ that satisfies the assumptions of § 2.1. Consider the vector potential

(4.3)\begin{equation} \boldsymbol{A}_S(\boldsymbol{x}) ={-}\frac 1{4{\rm \pi}} \iint_\varGamma \frac{\boldsymbol{n}(\boldsymbol{y}) \times \boldsymbol{B}^{\rm pol}(\boldsymbol{y})} {\left\| \boldsymbol{x} - \boldsymbol{y} \right\|} \,\mathrm{d} \varGamma(\boldsymbol{y}) \end{equation}

for $\boldsymbol {x} = (R,\phi _0,Z) = (r(t_0), \phi _0, z(t_0))$ in cylindrical coordinates. Define the quantities

(4.4a,b)\begin{equation} \alpha = \alpha(t; \boldsymbol{x}) = r(t)^2 + R^2 + (Z - z(t))^2 \quad \text{and} \quad \beta = \beta(t; \boldsymbol{x}) = 2 R r(t) \end{equation}

and set the modulus

(4.5)\begin{equation} k^2 = k(t; \boldsymbol{x})^2 = \frac{2\beta}{\alpha + \beta} = \frac{4 R r(t)}{(R + r(t))^2 + (Z - z(t))^2}. \end{equation}

Then the vector potential can be expressed as

(4.6)\begin{equation} \boldsymbol{A}_S(\boldsymbol{x}) ={-}\boldsymbol{e}_\phi(\phi_0) \int_0^L \mathcal{A}(t; \boldsymbol{x}) \,\mathrm{d} t \end{equation}

with the scalar-valued integrand

(4.7)\begin{equation} \mathcal{A}(t; \boldsymbol{x}) = \frac{1}{4{\rm \pi}} \left( \frac{4}{\sqrt{\alpha + \beta}} \right) \left[\dfrac{\partial \psi}{\partial z} r'(t) - \dfrac{\partial \psi}{\partial r} z'(t) \right] \left( \frac{2}{k^2} \left(K(k^2) - E(k^2) \right) - K(k^2) \right). \end{equation}

Proof. The unit outward normal vector of $\varGamma$ at $\boldsymbol {y}$ is given (by our assumption on the orientation of $\gamma$) by

(4.8)\begin{align} \boldsymbol{n}(\boldsymbol{y}(\phi,t)) & =\left. \frac{\partial{\boldsymbol{y}}}{\partial{\phi}}\times\frac{\partial{\boldsymbol{y}}}{\partial{t}} \right/ \left\| \frac{\partial{\boldsymbol{y}}}{\partial{\phi}}\times\frac{\partial{\boldsymbol{y}}}{\partial{t}} \right\| \nonumber\\ & = \frac{z'(t) \boldsymbol{e}_r(\phi) - r'(t) \boldsymbol{e}_z} {\sqrt{r'(t)^2 + z'(t)^2}}. \end{align}

From the representation (3.1) of $\boldsymbol {B}^{{\rm pol}}$ in axisymmetry, it follows that

(4.9)\begin{equation} (\boldsymbol{n}(\boldsymbol{y}(\phi,t)) \times \boldsymbol{B}^{{\rm pol}}(\boldsymbol{y}(\phi,t))) \mathcal{J}(t) = \left[\dfrac{\partial \psi}{\partial z} r'(t) - \dfrac{\partial \psi}{\partial r} z'(t) \right] \boldsymbol{e}_\phi(\phi). \end{equation}

Next, we consider the difference

(4.10)\begin{align} \boldsymbol{x} - \boldsymbol{y}(\phi, t) & = [R\boldsymbol{e}_r(\phi_0) + Z\boldsymbol{e}_z] - [r(t)\boldsymbol{e}_r(\phi)+ z(t)\boldsymbol{e}_z] \nonumber\\ & = ( R\cos\phi_0 - r(t)\cos\phi ) \boldsymbol{e}_x + ( R\sin\phi_0 - r(t)\sin\phi ) \boldsymbol{e}_y + (Z- z(t))\boldsymbol{e}_z, \end{align}

which we have expressed in both cylindrical and rectangular coordinates. Recall that the unit vector $\boldsymbol {e}_z$ is identical in both coordinate systems, and the remaining rectangular unit vectors $\boldsymbol {e}_x$ and $\boldsymbol {e}_y$ are related to standard cylindrical unit vectors by

(4.11a,b)\begin{equation} \boldsymbol{e}_r(\phi) = \cos\phi \, \boldsymbol{e}_x + \sin\phi \, \boldsymbol{e}_y \quad \text{and} \quad \boldsymbol{e}_\phi(\phi) ={-}\sin\phi \, \boldsymbol{e}_x + \cos\phi \, \boldsymbol{e}_y. \end{equation}

From the representation in rectangular coordinates, we use the trigonometric identity $\cos (\phi _0-\phi )=\cos \phi _0\cos \phi + \sin \phi _0\sin \phi$ and immediately obtain

(4.12)\begin{align} \left\| \boldsymbol{x} - \boldsymbol{y}(\theta, t) \right\|^2 = R^2 + r(t)^2 + (Z - z(t))^2 - 2 R r(t) \cos(\phi_0-\phi) = \alpha - \beta \cos(\phi_0 - \phi). \end{align}

We have now shown that the surface integral (4.1) for the vector potential is equivalent to the following double integral over a rectangle in the parameter domain:

(4.13)\begin{equation} \boldsymbol{A}_S(\boldsymbol{x}) ={-}\frac{1}{4{\rm \pi}} \int_0^L \left[\dfrac{\partial \psi}{\partial z} r'(t) - \dfrac{\partial \psi}{\partial r} z'(t) \right] \int_0^{2{\rm \pi}} \frac{\boldsymbol{e}_\phi(\phi)}{\sqrt{\alpha(t) - \beta(t) \cos(\phi_0 - \phi)}} \,\mathrm{d} \phi \,\mathrm{d} t. \end{equation}

We may analytically evaluate the inner integral using trigonometric identities and known integral formulae. We use the identities

(4.14)\begin{equation} \begin{cases} \sin(\phi+\phi_0) = \sin\phi\cos\phi_0 + \cos\phi\sin\phi_0\\ \cos(\phi+\phi_0) = \cos\phi\cos\phi_0 - \sin\phi\sin\phi_0 \end{cases} \end{equation}

to compute that

(4.15)\begin{align} \int_0^{2{\rm \pi}} \frac{-\sin\phi \, \boldsymbol{e}_x + \cos\phi \, \boldsymbol{e}_y} {\sqrt{\alpha - \beta\cos(\phi_0 - \phi)}} \,\mathrm{d} \phi & = \int_{0}^{2{\rm \pi}} \frac{-\sin(\phi+\phi_0) \, \boldsymbol{e}_x + \cos(\phi+\phi_0) \, \boldsymbol{e}_y} {\sqrt{\alpha - \beta\cos\phi}} \,\mathrm{d} \phi \nonumber\\ & = \boldsymbol{e}_\phi(\phi_0) \int_0^{2{\rm \pi}} \frac{\cos\phi \, \mathrm{d} \phi}{\sqrt{\alpha - \beta\cos\phi}} - \boldsymbol{e}_r(\phi_0) \int_0^{2{\rm \pi}} \frac{\sin\phi \, \mathrm{d} \phi}{\sqrt{\alpha - \beta\cos\phi}} \nonumber\\ & = \boldsymbol{e}_\phi(\phi_0) \int_0^{2{\rm \pi}} \frac{\cos\phi \, \mathrm{d} \phi}{\sqrt{\alpha - \beta\cos\phi}}. \end{align}

Here, we have used the fact that

(4.16)\begin{equation} \int_0^{2{\rm \pi}} \frac{\sin\phi \, \mathrm{d} \phi}{\sqrt{\alpha - \beta \cos\phi}} = 0, \end{equation}

which holds because it is the integral of an odd function over a single period.

Now, the remaining integral can be expressed in terms of the complete elliptic integrals, through the identity

(4.17)\begin{align} \int_0^{2{\rm \pi}} \frac{\cos\phi \, \mathrm{d} \phi}{\sqrt{\alpha - \beta\cos\phi}} & = 2\int_{-{\rm \pi}/2}^{{\rm \pi}/2} \frac{\cos (2\phi+{\rm \pi}) \, \mathrm{d} \phi}{\sqrt{\alpha - \beta\cos (2\phi+{\rm \pi})}} \nonumber\\ & = 2\int_{-{\rm \pi}/2}^{{\rm \pi}/2} \frac{-(1 - 2\sin^2\phi) \, \mathrm{d} \phi}{\sqrt{\alpha + \beta(1 - 2\sin^2\phi)}} \nonumber\\ & = \frac{-4}{\sqrt{\alpha+\beta}} \int_0^{{\rm \pi}/2} \frac{1-2\sin^2\phi} {\sqrt{1 - k^2 \sin^2\phi}} \mathrm{d} \phi \nonumber\\ & = \frac{-4}{\sqrt{\alpha+\beta}} \left( K(k^2) - \frac 2{k^2} (K(k^2) - E(k^2))\right). \end{align}

We obtain the last equality (4.17) from the definition of $K$ and from a formula of Gradshteyn & Ryzhik (Reference Gradshteyn and Ryzhik2014, 2.584-4). The desired result now immediately follows.

The Kapur–Rokhlin quadrature rule that we analysed in § 2.3.2 applies to integrands with a logarithmic singularity. With the result we just obtained, we can verify, by analysing the behaviour of the complete elliptic integrals, that the integrand $\mathcal {A}$ is indeed logarithmically singular as $t \to t_0$, in agreement with the results from Chance (Reference Chance1997, § V). Specifically, as the modulus $k$ tends to $1$, the second-kind integral $E$ is continuous and bounded, and the first-kind elliptic integral $K$ is logarithmically singular. The mapping $t \mapsto k(t; \boldsymbol {x})^2$ is continuous, so $E(k(t;\boldsymbol {x})^2)$ is continuous and $K(k(t;\boldsymbol {x})^2)$ is logarithmically singular as $t \to t_0$. Readers interested in more detail regarding these results are referred to Lemma A.2 in the Appendix.

6.1. Double-layer identity

For our first numerical verification, we consider an identity associated with harmonic functions. Consider the Green's function $G(\boldsymbol {x}, \boldsymbol {y}) = (4 {\rm \pi}\left \| \boldsymbol {x} - \boldsymbol {y} \right \|)^{-1}$ for Laplace's equation in three dimensions. It satisfies the double-layer jump condition (Malhotra et al. Reference Malhotra, Cerfon, O'Neil and Toler2019b)

(6.3)\begin{equation} \iint_\varGamma \frac{\partial G(\boldsymbol{x}, \boldsymbol{y})}{\partial \boldsymbol{n}(\boldsymbol{y})} \,\mathrm{d} \varGamma(\boldsymbol{y}) = \frac 1{4{\rm \pi}} \iint_\varGamma \frac{\boldsymbol{n}(\boldsymbol{y}) \boldsymbol{ \cdot} (\boldsymbol{x} - \boldsymbol{y})} {\left\| \boldsymbol{x} - \boldsymbol{y} \right\|^3} \,\mathrm{d} \varGamma(\boldsymbol{y}) ={-} \frac 12 \end{equation}

for $\boldsymbol {x} \in \varGamma$. It follows that

(6.4)\begin{equation} 1 + \frac 1{2{\rm \pi}} \iint_\varGamma \frac{\boldsymbol{n}(\boldsymbol{y}) \boldsymbol{ \cdot} (\boldsymbol{x} - \boldsymbol{y})} {\left\| \boldsymbol{x} - \boldsymbol{y} \right\|^3} \,\mathrm{d} \varGamma(\boldsymbol{y}) = 0, \end{equation}

again for $\boldsymbol {x} \in \varGamma$. Following identical methodology to the proof of Theorem 4.1, we integrate out the toroidal angle analytically to obtain the one-dimensional integral identity

(6.5)\begin{align} 1 + \frac 1{2{\rm \pi}} \int_0^L \frac{4r}{(\alpha+\beta)^{3/2}} \left\{ -\frac{2 z' R}{k^2} K(k^2) + \left( \frac{2 z' R}{k^2} + \frac{z'(R-r) - r'(Z-z)}{1-k^2} \right) E(k^2) \right\} \mathrm{d} t = 0. \end{align}

In this expression, we have suppressed the dependence of $\{r,z,r',z',\alpha,\beta,k^2\}$ on $t$ for clarity. As usual, we have also used the identification $\boldsymbol {x} = (R, \phi _0, Z) = (r(t_0), \phi _0, z(t_0))$. The integrand is logarithmically singular because of the presence of the singular elliptic integral $K(k^2)$, and because the other seemingly singular coefficient is actually bounded when $r''(t_0)$ and $z''(t_0)$ exist (which is the case in our example), with

(6.6)\begin{equation} \lim_{t \to t_0} \frac{z'(t)(R-r(t)) - r'(t)(Z-z(t))}{1-k(t)^2} = 2 R^2 \left( \frac{r''(t_0) z'(t_0) - r'(t_0) z''(t_0)}{r'(t_0)^2 + z'(t_0)^2} \right). \end{equation}

We conclude that the integral in (6.5) is of the Kapur–Rokhlin form from § 2.3.2.

Figure 3 illustrates the performance of the 10th-order periodic Kapur–Rokhlin quadrature scheme for verifying the identity (6.5). We compare this method with the alternating trapezoidal rule, which uses common quadrature weight $h={\rm \pi} /M$ and quadrature nodes $\tilde {t}_i = t_0 + (i - \frac 12) h$ for $i=0, \pm 1, \dots, \pm (M-1), M$ that straddle the singularity at $t_0$. We find that the Kapur–Rokhlin scheme achieves the theoretical 10th-order accuracy, and that we obtain the full accuracy of the quadrature scheme using about 175 quadrature nodes. Because the behaviour of the integrand is not generally symmetric about $t=t_0$, the alternating trapezoidal rule performs poorly. In order to avoid unfairly representing the performance of singularity subtraction schemes currently used in plasma physics applications, we have not attempted to code our own versions of them to compare their accuracy and their convergence rate in this figure. However, we explain in Malhotra et al. (Reference Malhotra, Cerfon, O'Neil and Toler2019b) why all these singularity subtraction schemes are expected have low-order convergence, with second-order convergence expected for most of them.

Figure 3. Error convergence for the double-layer identity (6.5) at $t_0 = 1$.

6.2. Virtual-casing principle

For our second test, we evaluate the accuracy of our method for calculating the normal component of the magnetic field due to the plasma current on the plasma boundary, i.e. $\boldsymbol {n} \boldsymbol {\cdot } \boldsymbol {B}^{{\rm pol}}_V$ on $\varGamma$. The plasma equilibrium we consider is the Solov'ev equilibrium described previously (Lütjens et al. Reference Lütjens, Bondeson and Sauter1996; Lee & Cerfon Reference Lee and Cerfon2015). As we do not know the analytic solution to this problem, we compare our implementation with an existing high-order accurate implementation from Malhotra et al. (Reference Malhotra, Cerfon, O'Neil and Toler2019b). The method used in that implementation is different from the approach presented here in several ways, making it appropriate for our verification. Specifically, the code presented in Malhotra et al. (Reference Malhotra, Cerfon, O'Neil and Toler2019b) views the plasma equilibrium as a fully three-dimensional equilibrium, and does not assume axisymmetry. Furthermore, Malhotra et al. (Reference Malhotra, Cerfon, O'Neil and Toler2019b) obtain $\boldsymbol {n} \boldsymbol {\cdot } \boldsymbol {B}^{{\rm pol}}_V$ on $\varGamma$ directly via a direct evaluation of (3.3), as opposed to first computing the vector potential. Finally, the Cauchy principal value of (3.3) is numerically evaluated via a partition of unity scheme to handle the singulary of the integrand. We take the result from a high-resolution calculation of Malhotra et al. (Reference Malhotra, Cerfon, O'Neil and Toler2019b) for this problem as the ground truth, against which we test the accuracy of our approach as a function of the number of quadrature points.

Figure 4 illustrates our results when using the 10th-order Kapur–Rokhlin scheme for this problem. Again, we find that the Kapur–Rokhlin scheme achieves the theoretical convergence rate, and that the accuracy has converged by about 400 quadrature nodes. In both this example and the double-layer identity example, we observe that the Kapur–Rokhlin scheme errors do not converge all the way to machine precision. This is one known drawback of the Kapur–Rokhlin methods, driven partially by the instabilities caused by the correction weights. Nonetheless, in contexts where full machine precision is not necessary, this scheme provides high accuracy and is easy to implement by reading the correction weights $\gamma _j$ from a table of pre-computed values.

Figure 4. Comparison with existing code to compute the normal component of the magnetic field.