Theory and numerical implementation

Theory

The code is based on an unpublished 1995 note by Steve Hirshman, "Transformation from VMEC to Boozer coordinates". (The method is not specific to VMEC, and can be applied to other equilibrium representations.) The calculation has two steps. In the first step, the difference is computed between the original toroidal and poloidal angles and the Boozer angles. In the second step, this information is used to transform other quantities such as \(R\), \(Z\), and \(B\) that are known as functions of the original angles to functions of the new angles.

1. Determining the toroidal angle difference

We assume that at the start, each flux surface is parameterized by a poloidal angle \(\theta_0\) and a toroidal angle \(\zeta_0\), which need not be straight-field-line angles. The toroidal angle need not be the standard toroidal angle (although it typically is). We assume that we know the quantity \(\lambda\) by which the poloidal angle can be shifted to obtain the straight-field-line angle \(\theta^*=\theta_0 + \lambda\) compatible with \(\zeta_0\). Thus, the magnetic field \(\vec{B}\) can be written as

(1)\[\vec{B} = \nabla\psi\times\nabla\theta^* + \iota(\psi) \nabla\zeta_0\times\nabla\psi\]

where \(2 \pi \psi\) is the toroidal flux, and \(\iota\) is the rotational transform. The Boozer poloidal angle \(\theta_B\) and toroidal angle \(\zeta_B\) are also straight field line angles, so

(2)\[\vec{B} = \nabla\psi\times\nabla\theta_B + \iota(\psi) \nabla\zeta_B\times\nabla\psi.\]

Moreover, Boozer coordinates satisfy

(3)\[\vec{B} = \beta(\psi,\theta_B,\zeta_B)\nabla\psi + I(\psi)\nabla\theta_B + G(\psi)\nabla\zeta_B.\]

We define \(\nu\) to be the difference between the two toroidal angles:

(4)\[\zeta_B = \zeta_0 + \nu.\]

Plugging this expression into (2) and equating the result with (1) gives \(\nabla\psi\times\nabla\left( \theta_B - \theta_0^* - \iota \nu \right) = 0\). Therefore \(\theta_B - \theta_0^* - \iota \nu = y(\psi)\) for some flux function \(y(\psi)\). We are free to choose \(y(\psi)=0\), which amounts to specifying the origin of \(\theta_B\) on each surface. Thus,

(5)\[\theta_B = \theta_0 + \lambda + \iota \nu.\]

Note

In Hirshman’s note “Transformation from VMEC to Boozer coordinates”, the angle difference is defined as \(p = \zeta_B - \zeta_0\). However in the fortran booz_xform code, a minus sign appears on line 83 of boozer.f that reverses the sign, so the p quantity saved in boozmn_*.nc files is in fact \(\zeta_0 - \zeta_B\). In the C++/python booz_xform code, the sign convention \(\nu = \zeta_B - \zeta_0\) is used consistently throughout the code and output files.

Plugging (4)-(5) into the covariant Boozer representation (3) gives

(6)\[\vec{B} = G \nabla \zeta_0 + (\nabla \theta_0 + \nabla \lambda) I + (G + \iota I) \nabla \nu + \left(\beta + I \nu \frac{d \iota}{d\psi}\right) \nabla \psi.\]

Note that the covariant components of \(\vec{B}\) in the original coordinates are

(7)\[ \begin{align}\begin{aligned}\begin{split}B_{\theta_0} = \vec{B} \cdot \frac{\partial\vec{r}}{\partial\theta_0} = \vec{B}\cdot\frac{\nabla\zeta_0\times\nabla\psi}{\nabla\psi\cdot\nabla\theta_0\cdot\nabla\zeta_0}, \\\end{split}\\B_{\zeta_0} = \vec{B} \cdot \frac{\partial\vec{r}}{\partial\zeta_0} = \vec{B}\cdot\frac{\nabla\psi\times\nabla\theta_0}{\nabla\psi\cdot\nabla\theta_0\cdot\nabla\zeta_0},\end{aligned}\end{align} \]

where \(\vec{r}\) is the position vector, and the dual relations have been used to get the right expressions from the central ones. Plugging (6) into these expressions gives

(8)\[B_{\theta_0} = \left( 1 + \frac{\partial\lambda}{\partial\theta_0}\right) I + (G + \iota I) \frac{\partial\nu}{\partial\theta_0}\]

and

(9)\[B_{\zeta_0} = G + I \frac{\partial\lambda}{\partial\zeta_0} + (G + \iota I) \frac{\partial\nu}{\partial\zeta}.\]

Equations (8) and (9) determine \(\nu\) up to a flux function (i.e. a constant on each surface). From (4), it can be seen that this constant effectively determines the origin of the \(\zeta_B\) coordinate. We are free to fix this constant by requiring that the average of \(\nu\) over \(\theta_0\) and \(\zeta_0\) be zero. We can write a double Fourier series for \(\nu\),

(10)\[\nu(\psi,\theta_0,\zeta_0) = \sum_{m,n} \left[ \hat{\nu}_{m,n}^s(\psi) \sin(m \theta_0 - n \zeta_0) + \hat{\nu}_{m,n}^c(\psi) \cos(m \theta_0 - n \zeta_0) \right],\]

with analogous series for \(\lambda\), \(B_{\theta_0}\), and \(B_{\zeta_0}\). The \(m=n=0\) modes of (8)-(9) then give

(11)\[I = \hat{B}_{\theta_0, 0, 0}^c, \;\; G = \hat{B}_{\zeta_0, 0, 0}^c,\]

allowing the flux functions \(I\) and \(G\) to be computed from the input data. When at least one of \(m\) or \(n\) is nonzero, the \(\cos(m \theta_0 - n \zeta_0)\) modes of (8)-(9) are

(12)\[\hat{\nu}_{m,n}^s = \frac{1}{G+\iota I} \left( \frac{\hat{B}_{\theta_0, m, n}^c}{m} - I \hat{\lambda}_{m,n}^s \right)\]

and

(13)\[\hat{\nu}_{m,n}^s = \frac{1}{G+\iota I} \left( \frac{-\hat{B}_{\zeta_0, m, n}^c}{n} - I \hat{\lambda}_{m,n}^s \right),\]

and the \(\sin(m \theta_0 - n \zeta_0)\) modes of (8)-(9) are

(14)\[\hat{\nu}_{m,n}^c = \frac{1}{G+\iota I} \left( \frac{-\hat{B}_{\theta_0, m, n}^s}{m} - I \hat{\lambda}_{m,n}^c \right)\]

and

(15)\[\hat{\nu}_{m,n}^c = \frac{1}{G+\iota I} \left( \frac{\hat{B}_{\zeta_0, m, n}^s}{n} - I \hat{\lambda}_{m,n}^c \right).\]

Equations (12)-(15) enable \(\nu\) to be calculated from the input data. To see that (12)-(13) are consistent with each other, the curl of \(\vec{B} = B_\psi\nabla\psi + B_{\theta_0}\nabla\theta_0 + B_{\zeta_0}\nabla\zeta_0\) can be plugged into the MHD equilibrium property \(\vec{J}\cdot\nabla\psi=0\), yielding \(\partial B_{\zeta_0}/\partial\theta_0 - \partial B_{\theta_0}/\partial\zeta_0=0\). In Fourier space, then,

(16)\[m \hat{B}_{\zeta_0,m,n}^s = -n \hat{B}_{\zeta_0,m,n}^s, \;\; m \hat{B}_{\zeta_0,m,n}^c = -n \hat{B}_{\zeta_0,m,n}^c.\]

The modes of \(\nu\) with \(m \ne 0\) can be computed from (12) and (14). The modes of \(\nu\) with \(n \ne 0\) can be computed from (13) and (15). For modes of \(\nu\) with both \(n \ne 0\) and \(m \ne 0\), we can use either (12) or (13), since the expressions are equivalent by (16); for the same reason we can use either (14) or (15). The mode of \(\nu\) with \(m=n=0\) can be set to zero, for as described above, this choice amounts to specifying the origin of \(\zeta_B\).

The results (12)-(15) can be summarized as

(17)\[\nu = \frac{w - I \lambda}{G + \iota I}\]

where \(w\) has a Fourier series analogous to (10), with \(\hat{w}_{0,0}^c=0\) (assuming \(\lambda_{0,0}^c=0\)) and

(18)\[\begin{split}\hat{w}_{m,n}^s = \hat{B}_{\theta_0,m,n}^c / m = -\hat{B}_{\zeta_0,m,n}^c / n, \\ \hat{w}_{m,n}^c = -\hat{B}_{\theta_0,m,n}^s / m = \hat{B}_{\zeta_0,m,n}^s / n.\end{split}\]

2. Transforming other quantities

Now that \(\nu\) is determined, the remaining task is to express other scalar quantities like \(B\) as functions of the new angles instead of the old angles. We let \(\Omega\) denote any scalar quantity that we wish to transform in this way, including the cylindrical coordinates \(R\) and \(Z\), and \(\nu\) itself. We write a double Fourier expansion, using a bar instead of hat to denote that the independent variables are now the Boozer angles:

(19)\[\Omega(\psi,\theta_B,\zeta_B) = \sum_{m,n} \left[ \bar{\Omega}_{m,n}^s(\psi) \sin(m \theta_B - n \zeta_B) + \bar{\Omega}_{m,n}^c(\psi) \cos(m \theta_B - n \zeta_B) \right].\]

The amplitudes satisfy

(20)\[\begin{split}\bar{\Omega}_{m,n}^s = \frac{1}{4\pi^2} \int_0^{2\pi}d\theta_B \int_0^{2\pi}d\zeta_B \; \Omega \sin(m \theta_B - n \zeta_B), \\ \bar{\Omega}_{m,n}^c = \frac{1}{4\pi^2} \int_0^{2\pi}d\theta_B \int_0^{2\pi}d\zeta_B \; \Omega \cos(m \theta_B - n \zeta_B).\end{split}\]

Changing the integration variables to the old angles,

(21)\[\begin{split}\bar{\Omega}_{m,n}^s = \frac{1}{4\pi^2} \int_0^{2\pi}d\theta_0 \int_0^{2\pi}d\zeta_0 \frac{\partial(\theta_B,\zeta_B)}{\partial(\theta_0,\zeta_0)} \Omega \sin( m[\theta_0 + \lambda + \iota \nu] - n[\zeta_0 + \nu]), \\ \bar{\Omega}_{m,n}^c = \frac{1}{4\pi^2} \int_0^{2\pi}d\theta_0 \int_0^{2\pi}d\zeta_0 \frac{\partial(\theta_B,\zeta_B)}{\partial(\theta_0,\zeta_0)} \Omega \cos( m[\theta_0 + \lambda + \iota \nu] - n[\zeta_0 + \nu]),\end{split}\]

where (4)-(5) have been applied. The Jacobian appearing in (21) is

(22)\[\begin{split}\frac{\partial(\theta_B,\zeta_B)}{\partial(\theta_0,\zeta_0)} &=\frac{\partial\theta_B}{\partial\theta_0} \frac{\partial\zeta_B}{\partial\zeta_0} -\frac{\partial\theta_B}{\partial\zeta_0} \frac{\partial\zeta_B}{\partial\theta_0} \\ &=\left(1 + \frac{\partial\lambda}{\partial\theta_0}\right) \left(1 + \frac{\partial\nu}{\partial\zeta_0}\right) +\left(\iota - \frac{\partial\lambda}{\partial\zeta_0}\right) \frac{\partial\nu}{\partial\theta_0}.\end{split}\]

It can be seen that (21) and the last expression in (22) express the new Fourier amplitudes we seek in terms of the old angles and known quantities on these angles.

Implementation details

Consistent with the discussion following (16), the modes of \(\nu\) are computed via \(w\) as follows. If \(m \ne 0\), (12) and (14) are used, equivalent to the first equality on each line of (18). If \(m = 0\) but \(n \ne 0\), (13) and (15) are used, equivalent to the last expression on each line of (18). This logic appears in surface_solve.cpp.

The integrals in (21) are computed using a uniform tensor-product grid in \(\theta_0\) and \(\zeta_0\).