Direct, simple and efficient computation of all components of the virtual-casing magnetic field in axisymmetric geometries with Kapur–Rokhlin quadrature

Abstract

In a recent publication (Toler et al., J. Plasma Phys., vol. 89, issue 2, 2023, p. 905890210), we demonstrated that for axisymmetric geometries, the Kapur–Rokhlin quadrature rule provided an efficient and high-order accurate method for computing the normal component, on the plasma surface, of the magnetic field due to the toroidal current flowing in the plasma, via the virtual-casing principle. The calculation was indirect, as it required the prior computation of the magnetic vector potential from the virtual-casing principle, followed by the computation of its tangential derivative by Fourier differentiation, to obtain the normal component of the magnetic field. Our approach did not provide the other components of the virtual-casing magnetic field. In this letter, we show that a more direct and more general approach is available for the computation of the virtual-casing magnetic field. The Kapur–Rokhlin quadrature rule accurately calculates the principal value integrals in the expression for all the components of the magnetic field on the plasma boundary, and the numerical error converges at a rate nearly as high as the indirect method we presented previously.

1. Introduction

The virtual-casing principle (Shafranov & Zakharov 1972; Zakharov 1973; Hanson 2015) is an effective tool for computing the magnetic field due to the currents flowing in the plasma, since it reduces the dimensionality of the integrals involved in this computation, from volume to surface integrals for non-axisymmetric geometries (Lazerson, Sakakibara & Suzuki 2013; Malhotra et al. 2019; Kappel, Landreman & Malhotra 2023) and from surface to line integrals for axisymmetric settings (Zakharov 1973; Zakharov & Pletzer 1999; Toler, Cerfon & Malhotra 2023). Computing the magnetic field with the virtual-casing principle is nevertheless challenging when the evaluation point is on the plasma boundary, because it involves the evaluation of integrals with singular integrands.

Singularity subtraction schemes have often been used by the plasma physics community to address this difficulty (Drevlak et al. 2018). They are robust, but are typically low-order accurate schemes (Malhotra et al. 2019), which are also tedious to implement. In contrast, there exist quadrature weights and abscissae in the applied mathematics literature which are specifically designed to give high accuracy for the singularities encountered in plasma physics applications, and are simple to implement. In recent work (Toler et al. 2023), we demonstrated that for axisymmetric geometries, the Kapur–Rokhlin quadrature scheme (Kapur & Rokhlin 1997) provided high-order accuracy for the computation of the normal component on the plasma surface of the virtual-casing magnetic field. Since the Kapur–Rokhlin scheme is designed for integrable singularities with known asymptotic behaviour, we constructed this calculation in two steps. We first computed the virtual-casing vector potential, which can be expressed in terms of an integral with a known logarithmic singularity. Then, we calculated the normal component of the magnetic field on the plasma surface by computing the tangential derivative of the vector potential, by Fourier differentiation.

In the present letter, we show that the Kapur–Rokhlin quadrature scheme can be used in a more direct and more general way for the virtual-casing principle in axisymmetric geometries. Specifically, we demonstrate that the Kapur–Rokhlin scheme accurately computes the Cauchy principal value integrals arising in the expressions for the magnetic field in the virtual-casing principle, allowing us to compute all components of the magnetic field directly anywhere on the plasma boundary, via line integrals that are straightforward to implement and do not require any manipulation of the integrand.

Our previous article and the present letter address different needs associated with different applications. In certain situations, only the virtual-casing vector potential, or equivalently the virtual-casing poloidal magnetic flux function, is required (Lao et al. 2005; Marx & Lütjens 2017; Pustovitov & Chukashev 2021), in which case the method we presented in our first article is the most appropriate method. In other situations, knowledge of the component of the magnetic field normal to the plasma surface is desired (Landreman & Boozer 2016; Landreman 2017). In these cases, both our previous article and the present letter offer satisfactory solutions, although the approach presented here has the advantage of yielding the normal component of the magnetic field directly, bypassing the Fourier differentiation of the vector potential. Finally, there exist situations in which the tangential component of the magnetic field on the plasma surface is required, for example, as a Neumann boundary condition for the poloidal magnetic flux function (Zaitsev et al. 2011; Blum, Boulbe & Faugeras 2012; Ricketson et al. 2016). The method we present here is uniquely well suited for these situations.

The letter is structured as follows. In § 2, we review the virtual-casing principle for axisymmetric settings and highlight the relevant line integrals for computing the magnetic field due to the plasma, evaluated on the plasma boundary. Section 3 describes numerical quadrature schemes for evaluating the principal value integrals that appear in these line integrals. In § 4, we illustrate the performance of the quadrature schemes on an example determined by an analytic Grad–Shafranov equilibrium, and demonstrate the high order convergence and excellent accuracy of the Kapur–Rokhlin scheme. We conclude in § 5.

2. Axisymmetric virtual-casing principle

2.1. Cylindrical coordinates for toroidally axisymmetric surfaces

Throughout the letter, we will be concerned with plasmas with toroidally axisymmetric boundaries. We shall therefore make use of the standard, right-handed cylindrical coordinates $(r, \phi, z)$ naturally associated with the toroidal geometry, where the $z$-axis is the axis of revolution. 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 emphasize the fact that the radial and azimuthal unit vectors depend on the toroidal angle.

Let $\gamma$ be the simple closed curve in the $(r,z)$ plane such that by rotating $\gamma$ about the $z$-axis for $\phi \in [0, 2{\rm \pi} ]$, we obtain the closed surface of revolution $\varGamma$ corresponding to the plasma boundary. The generating curve $\gamma$ is parametrized 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 parametrization. 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_{\min }>0$ for which $r(t) \ge R_{\min }$ on $[0,L]$.

2.2. Virtual-casing principle for axisymmetric plasmas

Consider an axisymmetric plasma confined by external coils in equilibrium. Let $\boldsymbol {B}$ denote the total magnetic field. The poloidal field $\boldsymbol {B}^{{\rm pol}} = \boldsymbol {B} - B_{\phi } \boldsymbol {e}_{\phi }$ 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:

(2.1)\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(\phi(\boldsymbol{y})) \times \frac{\boldsymbol{x} - \boldsymbol{y}}{\| \boldsymbol{x} - \boldsymbol{y}\|^3} \mathrm{d} \boldsymbol{y}, \end{equation}

where $\mu _{0}$ is the permeability of free space, the integration volume $\varOmega$ is the plasma domain and $J^{{\rm tor}}_{V}$ is the toroidal current density in the plasma. The vector $\boldsymbol {e}_{\phi }(\phi (\boldsymbol {y}))$ is the unit vector in the toroidal direction at $\boldsymbol {y}$, where the argument $\phi (\boldsymbol {y})$ is the toroidal angle of the point $\boldsymbol {y}$.

Equation (2.1) 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 1972; Zakharov 1973; Hanson 2015). 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}}$, where $\boldsymbol {n}$ is the outward unit normal vector to $\varGamma$. This relation is given according to Hanson (2015):

(2.2)\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})}{\| \boldsymbol{x} - \boldsymbol{y}\|^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}

In this letter, we are interested in the most common application of the virtual-casing principle, in which we want to compute the magnetic field due to the plasma currents on the plasma boundary, i.e. for $\boldsymbol {x}\in \varGamma$. That motivates why we focus on the poloidal magnetic field. Indeed, for axisymmetric plasma boundaries and magnetic fields, the toroidal magnetic field on the plasma boundary is entirely generated by external coils. To see this, let $C_{r,z}$ be any circle of radius $r$ in a plane perpendicular to the $z$-axis, at altitude $z$, that is tangent to the axisymmetric plasma boundary. Applying Stokes’ theorem along this circle, we can write

(2.3)\begin{equation} \oint_{C_{r,z}}B_{\phi}(r,z)r\,\mathrm{d}\phi = 2{\rm \pi} r B_{\phi}(r,z) = \mu_{0}\iint_{D_{r,z}} \boldsymbol{J} \boldsymbol{\cdot} \mathrm{d}\boldsymbol{S} = \mu_{0} I_{\mathrm{coils}}, \end{equation}

where $D_{r,z}$ is the flat disc perpendicular to the $z$-axis bounded by the circle $C_{r,z}$, and $I_{\mathrm {coils}}$ is the total coil-induced current. In the last equality, we have used the fact that the plasma currents are divergence free. We obtain the well-known expression for the vacuum toroidal field, $B_{\phi }(r)=\mu _{0}I_{\mathrm {coils}}/(2{\rm \pi} r)$, which holds all the way to the plasma boundary and confirms that for axisymmetric plasmas, plasma currents do not contribute to the toroidal magnetic on the plasma boundary.

2.3. Reduction of the virtual-casing principle to line integrals

Since we are considering an axisymmetric surface $\varGamma$, the surface integral in (2.2) can be expressed in terms of a line integral by integrating over the toroidal angle analytically.

The 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 parametrization $(\phi,t) \mapsto \boldsymbol {y}(\phi, t)$ by (Freidberg 2014)

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

Inserting this expression in (2.2) and integrating it analytically with respect to the toroidal angle, we may evaluate $\boldsymbol {B}_V^{{\rm pol}}$ directly through the remaining univariate integral. Presently, we consider evaluating on the plasma boundary $\varGamma$. Without loss of generality, we assume that the evaluation point lies on the cross-section of the plasma boundary with azimuthal angle $\phi =0$, which has cylindrical coordinate unit vectors that align with the rectangular unit vectors $\{\boldsymbol {e}_x, \boldsymbol {e}_y, \boldsymbol {e}_z\}$. We denote such an evaluation point as $(R,0,Z)$ in cylindrical coordinates.

Upon integrating the toroidal angle, the univariate expression for $\boldsymbol {B}_V^{{\rm pol}}$ on the plasma boundary is

(2.5)\begin{equation} \boldsymbol{B}_V^{{\rm pol}}(R,Z)= \frac{1}{4{\rm \pi}} \int_{0}^{L} \boldsymbol{f}(t)\,\mathrm{d} t + \frac{1}{2} \boldsymbol{B}^{{\rm pol}}(R,Z),\quad (R,0,Z) \in \varGamma, \end{equation}

where (Toler et al. 2023)

(2.6)\begin{align} \boldsymbol{f}(t) & = \frac {2} {r(t) \sqrt{\alpha + \beta}} \left[ \frac{\partial \psi}{\partial z} r'(t) - \frac{\partial \psi}{\partial r} z'(t) \right] \left\{\frac {Z - z(t)}{R}\left(-K(k^2)+ \frac{\alpha}{\alpha - \beta} E(k^2)\right) \right.\boldsymbol{e}_x \nonumber\\ & \quad \left.+ \left(K(k^2)+ \frac{r(t)^2 - R^2 - (Z-z(t))^2}{\alpha - \beta} E(k^2)\right)\boldsymbol{e}_z\right\}\!. \end{align}

Here, the functions $K$ and $E$ are the complete elliptic integrals of the first and second kind, primes denote derivatives with respect to the parametrization variable $t$, and we have introduced the quantities

(2.7)\begin{equation} \begin{cases} \alpha = \alpha(t; R,Z) = R^2 + r(t)^2 + (Z-z(t))^2, \\ \beta = \beta(t; R,Z) = 2 R r(t), \\ k^2 = k(t; R,Z)^2 = \dfrac{2\beta}{\alpha+\beta}. \end{cases} \end{equation}

Since we are evaluating on the plasma boundary, the evaluation point can be identified by its generating curve parameter $t$, say $t_0$, so $r(t_0)=R$ and $z(t_0)=Z$. There are two sources of singularity in $\boldsymbol {f}(t)$ near $t_0$. The mapping $t \mapsto K(k(t)^2)$ is $O(\log |t-t_0|)$ near $t_0$ (Gradshteyn & Ryzhik 2014; Toler et al. 2023), and the errors from Taylor series for $r(t)$ and $z(t)$ yield $1/(\alpha -\beta ) = O(1/(t-t_0)^2)$.

3. High-order principal value quadrature for the axisymmetric virtual-casing principle

Previously (Toler et al. 2023), we showed that the integral in (2.5) is not convergent. That means that this integral must be interpreted in the Cauchy principal value sense for all components of $\boldsymbol {B}_V^{{\rm pol}}(R,Z)$ in (2.5) to be well defined. Since the observation point $(R,Z)$ lies on the plasma boundary, we may identify this point with a value of the boundary parameter $t$, say $t_0$.

To evaluate all the components of $\boldsymbol {B}_V^{{\rm pol}}(R,Z)$ via (2.5), we must evaluate the integral in this expression with the following substitution:

(3.1)\begin{equation} \int_{0}^{L} \boldsymbol{f}(t)\,\mathrm{d} t\quad \xrightarrow{\qquad} \quad \lim_{\epsilon \to 0^+} \left(\int_{0}^{t_0-\epsilon} \boldsymbol{f}(t)\,\mathrm{d} t + \int_{t_0+\epsilon}^{L} \boldsymbol{f}(t)\,\mathrm{d} t \right)\!, \end{equation}

where $\boldsymbol {f}$ is given in (2.6). (If $t_0=0$, the integration interval within the limit should be $[\epsilon, L-\epsilon ]$.) Let us discuss two quadrature rules to do so: the alternating trapezoidal rule and the periodic Kapur–Rokhlin rule (Kapur & Rokhlin 1997; Toler et al. 2023).

We begin by considering the alternating trapezoidal rule. We recall that the standard $N$-point periodic trapezoidal rule for the integral $\int _{0}^{L} g(t)\,\mathrm{d} t$ for an $L$-periodic integrand $g$ is defined by equispaced quadrature nodes with uniform weights $h=L/N$. The alternating trapezoidal rule for the case when $g$ has a singularity at $t_0$ is a specific case of the trapezoidal rule, where the quadrature nodes $\{s_i\}_{i=1}^{N}$ are placed to straddle $t_0$. That is, the nodes are $s_i = t_0 + (i - \frac {1}{2}) h$ for $i=1, \ldots, N$. The alternating trapezoidal rule for some types of Cauchy principal value integrals is supported by mathematical theory. In particular, it is known to converge spectrally accurately for principal value integrals of the form $\int _{-1}^1 (q(t)/t)\,\mathrm{d} t$ when $q$ is an analytic function (Kress & Martensen 1970; Trefethen & Weideman 2014). Additionally, if $q \in C^{P}[0,L]$ but is not analytic, then the quadrature converges with order $P$ (Sidi & Israeli 1988). By staggering quadrature nodes with equal weights on each side of the singularity, the alternating trapezoidal rule captures the cancellation that is necessary for a well-defined Cauchy principal value. For our application, we note that the asymptotic behaviour of the integrand is not $q(t)/t$, but $q_1(t)/t + q_2(t) \log |t| + q_3(t)$ for bounded functions $q_1, q_2, q_3$. As a result of the logarithmically singular term, we do not expect high-order convergence of the quadrature error. Rather, we expect low-order convergence, as we will verify numerically in the next section.

In contrast, the Kapur–Rokhlin rule is precisely designed to integrate logarithmic singularities with high-order accuracy (Toler et al. 2023). We recall that the $(N-1)$-point rule for the periodic setting of interest in our context can be written as

(3.2)\begin{equation} \int_{t_0-L/2}^{t_0+L/2} g(t)\,\mathrm{d} t = \sum_{\substack{j=-(N-1)/2 \\ j \ne 0}}^{(N-1)/2} w_j g_j + O(h^n), \end{equation}

where $g$ is logarithmically singular at $t_0$, with

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

and $g_j \equiv g(t_0 + j h)$. The convergence order $n$ may be chosen by the user and the weights $\{ \gamma _j \}$ may be read from a table of values (Kapur & Rokhlin 1997). The spacing $h=L/N$ between quadrature points is uniform (except for the doubled space between nodes $-1$ and $1$), and we are using notation so that the meaning of $h$ is the same for both the alternating trapezoidal rule and the Kapur–Rokhlin rule. In (3.2), we assume that $N$ is odd to simplify notation. If $N$ is even, the points may be indexed by $j=\pm 1, \ldots, \pm (N/2-1), N/2$.

Like the alternating trapezoidal rule, the periodic Kapur–Rokhlin rule staggers nodes around the singularity with equal weights, so it also captures the principal value cancellation when integrating $\boldsymbol {f}$. Indeed, the weights in (3.3) satisfy $w_{j} = w_{-j}$. However, the quadrature error converges much faster than for the alternating trapezoidal rule, since the quadrature rule was also designed for the logarithmic singularity of $\boldsymbol {f}$. We next discuss practical performance of these numerical methods in the context of a virtual-casing principle calculation.

4. Numerical results

To compare and verify the numerical methods for calculating the principal value of the integral in (2.5), we use the same reference data as Toler et al. (2023) on the same geometry. Specifically, the axisymmetric plasma boundary is given by the level set ${\{ \psi =0 \}}$ of the poloidal flux function given by

(4.1)\begin{equation} \psi(r,z) = \frac{\kappa F_B}{2 R_0^3 q_0} \left[\frac 14 (r^2 - R_0^2)^2 + \frac 1{\kappa^2} r^2 z^2 - a^2 R_0^2 \right]\!, \end{equation}

which solves the Grad–Shafranov equation with the Solov'ev profiles $\mu _{0}p(\psi )=-[F_{B}(\kappa + 1/\kappa )/( R_0^{3}q_0)]\psi$ and $F(\psi )=F_B$, where $p(\psi )$ is the plasma pressure profile, and $F(\psi )=rB_{\phi }$, with $B_\phi$ the toroidal magnetic field (Lütjens, Bondeson & Sauter 1996; Lee & Cerfon 2015). The parameters $R_0$ and $q_0$ may be interpreted as the major radius and safety factor at the magnetic axis, and $\kappa$ and $a$ as the elongation and minor radius of the plasma boundary, respectively. Throughout, we use the fusion relevant values $F_B = R_0 = q_0 = 1$, and $\kappa = 1.7$ and $a = 1/3$. The level set $\{ \psi = 0 \}$ may be parametrized by the functions (Lütjens et al. 1996; Lee & Cerfon 2015)

(4.2a,b)\begin{equation} (r(t))^2 = R_0^2 + 2 a R_0 \cos t\quad \text{and}\quad z(t) = \kappa a \frac{R_0}{r(t)} \sin t \end{equation}

for $t \in [0, L]$ with $L=2{\rm \pi}$.

We take as ground truth the plasma-induced magnetic field from a high-resolution computation by the approach of Malhotra et al. (2019). The method used in that implementation is different from the one-dimensional integral approach presented here in several ways, making it appropriate for our verification. Specifically, the code presented by Malhotra et al. (2019) views the plasma equilibrium as a fully three-dimensional equilibrium and does not assume axisymmetry. Furthermore, Malhotra et al. (2019) obtain $\boldsymbol {B}^{{\rm pol}}_V$ on the plasma boundary directly by evaluating the virtual-casing principle as a surface integral, as opposed to first reducing to a univariate integral by axisymmetry. Finally, the Cauchy principal value of the surface integral is numerically evaluated via a partition of unity scheme to handle the singularity of the integrand. The reference solution was computed on a fine grid with 1,580 toroidal angle values and 1,200 poloidal points on each cross-section. In a self-convergence comparison with a solution from finer resolution, these data were found to be accurate to 11 digits of precision. Figure 1 illustrates the geometry and reference solution for this problem. This approach by Malhotra et al. (2019) is more robust by treating fully three-dimensional volumes with two-dimensional surfaces, but axisymmetric specializations like those we present here require discretization of one less dimension and result in faster runtimes.

Figure 1. Reference solution for the Solov'ev profiles of § 4. Colours represent $\| \boldsymbol {B}\|$ on the plasma surface.

With these reference values, we shall assess the numerical accuracy of the alternating trapezoidal rule and the Kapur–Rokhlin rule, which we discussed in the previous section. For each quadrature scheme, we compute both the radial component $\boldsymbol {B}_{V,R}^{{\rm pol}}$ of $\boldsymbol {B}_V^{{\rm pol}}$ and the vertical component $\boldsymbol {B}_{V,Z}^{{\rm pol}}$ at 1,200 points on the plasma boundary which are equispaced in the generating curve parameter $t$. We measure the relative quadrature error for each component as the maximum of the pointwise absolute difference between the computed field component and the reference values. We normalize both errors by the maximum absolute value of the 2,400 reference values for $\boldsymbol {B}_{V,R}^{{\rm pol}}$ and $\boldsymbol {B}_{V,Z}^{{\rm pol}}$.

Figure 2 illustrates that the alternating trapezoidal rule indeed converges to the same values found with the scheme presented by Malhotra et al. (2019), but does so with a low order of convergence. The convergence rates are $O(h^3)$ for the radial component $\boldsymbol {B}_{V,R}^{{\rm pol}}$ and $O(h)$ for the vertical component $\boldsymbol {B}_{V,Z}^{{\rm pol}}$. The radial component converges quicker because its integrand has more regularity due to the following asymptotic behaviour. Recall from § 2.3 that $t \mapsto K(k(t)^2)$ is $O(\log |t-t_0|)$ near $t_0$. Hence, in the integrand $\boldsymbol {f}$ of (2.6), the product $(Z-z(t))K(k^2)$ has only a removable discontinuity at $t_0$ and not a logarithmic singularity. The only singularity in the radial component is the principal value singularity of $(Z-z(t))/(\alpha -\beta )$. No such product dampens the singularity of $K$ in the vertical component. Nevertheless, the convergence rates for both components of $\boldsymbol {B}_{V}^{{\rm pol}}$ are better than those guaranteed by the quadrature theory we discussed in § 3.

Figure 2. Maximum pointwise error using the alternating trapezoidal rule to compute $\boldsymbol {B}_{V,R}^{{\rm pol}}$ and $\boldsymbol {B}_{V,Z}^{{\rm pol}}$ via principal value.

We next consider the periodic Kapur–Rokhlin rule for the principal value quadrature. We view each component of the integrand $\boldsymbol {f}$ as the sum of a logarithmically singular term and a remainder term:

(4.3)\begin{equation} \begin{cases} \boldsymbol{e}_x \boldsymbol{\cdot} \boldsymbol{f}(t) = p_R(t) \log|t-t_0| + q_R(t), \\ \boldsymbol{e}_z \boldsymbol{\cdot} \boldsymbol{f}(t) = p_Z(t) \log|t-t_0| + q_Z(t). \end{cases} \end{equation}

In particular, we have

(4.4)\begin{equation} p_R(t) \log|t-t_0| = \frac {2} {r(t) \sqrt{\alpha + \beta}} \left[ \frac{\partial \psi}{\partial z} r'(t) - \frac{\partial \psi}{\partial r}z'(t) \right] \frac {Z - z(t)}{R} (-K(k^2)) \end{equation}

and

(4.5)\begin{equation} p_Z(t) \log|t-t_0| = \frac {2} {r(t) \sqrt{\alpha + \beta}} \left[\frac{\partial \psi}{\partial z} r'(t) - \frac{\partial \psi}{\partial r} z'(t) \right]K(k^2), \end{equation}

and we set $q_R$ and $q_Z$ to be the remaining corresponding terms in (2.6). In figures 3, 4 and 5, we show the results of applying the second-, sixth- and tenth-order Kapur–Rokhlin rules, respectively. These are the theoretical convergence rates guaranteed for smooth coefficient functions $p_R, p_Z, q_R, q_Z$. In this example, however, we see convergence rates that deviate slightly from these values. The sixth-order rule converges at the expected rate of $O(h^6)$. However, the second-order rule converges faster than predicted at rates of approximately $O(h^{2.71})$ and $O(h^{2.52})$ for the radial and vertical components, respectively. Additionally, the tenth-order rule converges at approximate rates of $O(h^{8.74})$ and $O(h^{8.73})$. We compute these approximate rates by performing an ordinary least squares line of best fit to the data measurements in the asymptotic regime of convergence.

Figure 3. Maximum pointwise error using the second-order Kapur–Rokhlin rule to compute $\boldsymbol {B}_{V,R}^{{\rm pol}}$ and $\boldsymbol {B}_{V,Z}^{{\rm pol}}$ via principal value.

Figure 4. Maximum pointwise error using the sixth-order Kapur–Rokhlin rule to compute $\boldsymbol {B}_{V,R}^{{\rm pol}}$ and $\boldsymbol {B}_{V,Z}^{{\rm pol}}$ via principal value. We compare the Kapur–Rokhlin (KR) error with the alternating trapezoidal (AT) error of figure 2.

Figure 5. Maximum pointwise error using the tenth-order Kapur–Rokhlin rule to compute $\boldsymbol {B}_{V,R}^{{\rm pol}}$ and $\boldsymbol {B}_{V,Z}^{{\rm pol}}$ via principal value. We compare the Kapur–Rokhlin (KR) error with the alternating trapezoidal (AT) error of figure 2.

Notably, the tenth-order scheme converges slightly slower than predicted. We offer two potential explanations for this phenomenon. First, the smoothness assumptions of Kapur–Rokhlin are violated, as we explain presently. Similar to the alternating trapezoidal rule, the Kapur–Rokhlin rule requires that all of $\{p_R, p_Z, q_R, q_Z\}$ are sufficiently smooth to guarantee the theoretical convergence rate (Kapur & Rokhlin 1997). However, both $q_R$ and $q_Z$ include the elliptic integral $E$. The mapping $t \mapsto E(k(t)^2)$ is continuous, but its first derivative is not. This is readily seen by the identity of Gradshteyn & Ryzhik (2014) that

(4.6)\begin{equation} \frac{\partial E(k^2)}{\partial k^2} = \frac{E(k^2) - K(k^2)}{2 k^2}. \end{equation}

The first derivative of $t \mapsto E(k(t)^2)$, and hence the first derivatives of $q_R$ and $q_Z$, are logarithmically singular at $t_0$ due to the appearance of $K$. Second, very high orders of the Kapur–Rokhlin rule have weight corrections $\{\gamma _j\}$ which are large in magnitude and are sign-indefinite. This introduces numerical instability as a price to pay for the higher convergence order (Martinsson 2019).

We hypothesize that the instability from the weight corrections also explains the observation that the quadrature error saturates at approximately nine digits of accuracy in figures 4 and 5. It is numerically common for algorithmic instabilities to raise the error floor from machine precision. The instability is known to worsen as the theoretical convergence order increases, so we see this phenomenon earlier and more markedly for the tenth-order rule compared with the sixth-order rule. In contemporary plasma physics applications relevant to this work, however, accuracy at the level of eight or nine digits is typically not a limiting factor in overall calculations, so we do not explore the error saturation further.

Though the second- and sixth-order rules have agreed with or exceeded the theoretical baselines in our experiment, we emphasize that we cannot guarantee such performance in general. The behaviour of all of the Kapur–Rokhlin schemes in this setting are less predictable, but nevertheless promising in practice.

When we directly compare the performance of the Kapur–Rokhlin scheme and the alternating trapezoidal scheme in figures 4 and 5, we see that both methods are competitive for computing $\boldsymbol {B}^{{\rm pol}}_{V,R}$, whose virtual-casing integrand contains a less severe singularity. The Kapur–Rokhlin rule is asymptotically faster to converge, but the alternating trapezoidal rule has better accuracy when the number of quadrature points is small. However, the Kapur–Rokhlin rule more clearly outperforms the alternating trapezoidal rule for computing $\boldsymbol {B}^{{\rm pol}}_{V,Z}$, whose virtual-casing integrand contains the more severe singularity. In this case, the alternating trapezoidal rule only has higher accuracy when the number of quadrature points is very limited: less than approximately 50.

Finally, we consider a third quadrature scheme which directly combines the alternating trapezoidal rule and the Kapur–Rokhlin rule. Returning to the decomposition (4.3), we integrate the logarithmically singular terms with $p_R$ and $p_Z$ of (4.4) and (4.5) via the periodic Kapur–Rokhlin rule, and we integrate the remaining terms $q_R$ and $q_Z$ via the alternating trapezoidal rule. One might imagine that this scheme would outperform both the pure Kapur–Rokhlin rule and the pure alternating trapezoidal rule, since this combined scheme uses each constituent quadrature on the term for which it is specialized. However, figure 6 reveals that the lack of smoothness in the integrand again stalls convergence. For few quadrature nodes, the Kapur–Rokhlin corrections considerably improve convergence, but for more than approximately 100 nodes, the slower convergence of the alternating trapezoidal rule dominates. We observe that the limiting asymptotic convergence rate appears to be third-order convergence in this experiment for both the radial and vertical components of $\boldsymbol {B}_{V}^{{\rm pol}}$. This improves the asymptotic rate from the pure alternating trapezoidal rule for $\boldsymbol {B}_{V,Z}^{{\rm pol}}$, indicating that the alternating trapezoidal rule originally achieved only $O(h)$ convergence because it ignored the logarithmic singularity.

Figure 6. Maximum pointwise error using a combination of the tenth-order Kapur–Rokhlin rule and the alternating trapezoidal rule to compute $\boldsymbol {B}_{V,R}^{{\rm pol}}$ and $\boldsymbol {B}_{V,Z}^{{\rm pol}}$. We compare the error with the alternating trapezoidal (AT) error of figure 2.

We additionally note that the asymptotic rates for the second- and tenth-order Kapur–Rokhlin rules in figures 2 and 4 experimentally appear to agree with rates of $O(h^3 \log h)$ and $O(h^9 \log h)$, respectively, though this is our hypothesis and we are not aware of quadrature theory that would imply the presence of a $\log h$ factor in the convergence rate in this context. However, convergence orders with $\log h$ factors are known to appear when ill-suited quadrature rules, like the alternating trapezoidal rule, are used to integrate a function with a purely logarithmic singularity (Martinsson 2019). We display this feature for the tenth-order rule in figure 6.

The pure Kapur–Rokhlin method remains the most efficient quadrature scheme discussed for this principle value application. It empirically displays high-order convergence and accuracy close to the theorized rates. The tenth-order rule achieves a maximal accuracy of approximately nine digits with approximately 400 quadrature nodes.

5. Conclusion

We have shown that both the alternating trapezoidal rule and the Kapur–Rokhlin rule are convergent schemes for the computation of all the components of $\boldsymbol {B}_{V}^{{\rm pol}}$ expressed as one-dimensional Cauchy principal value integrals of the virtual-casing principle for axisymmetric equilibria. However, the convergence of the alternating trapezoidal rule is limited to very low order, which can be as low as one, because of the lack of smoothness of the integrand. In contrast, high-order Kapur–Rokhlin schemes, which are nearly as easy to implement as the alternating trapezoidal rule, provide high-order convergence for this calculation because the schemes are specifically designed for the singularities encountered in this application.

We therefore recommend the direct implementation of high-order Kapur–Rokhlin rules for the computation of the axisymmetric virtual-casing principle as presented in this letter, since it is easier to use than the method originally proposed by us in Toler et al. (2023), it provides all the components of the magnetic field and the loss of accuracy is negligible as compared to the results obtained by us in Toler et al. (2023) for the normal component of the magnetic field only.