Field-aligned coordinates
In order to efficiently simulate (predominantly) field-aligned
structures, grid-points are placed in a field-aligned coordinate system.
We define
\(\sigma_{B\theta} \equiv {B_{\text{pol}}}/ \left|{B_{\text{pol}}}\right|\)
i.e. the sign of the poloidal field. The new coordinates
\(\left(x,y,z\right)\) are defined by:
(11)\[\begin{aligned}
x = {\sigma_{B\theta}}\left(\psi - \psi_0\right) \qquad y = \theta \qquad z = \sigma_{B\theta}
\left(\zeta - \int_{\theta_0}^{\theta}\nu\left(\psi,\theta\right)d\theta\right)
\end{aligned}\]
Where \(\nu\) is the local field-line pitch given by
\[\begin{aligned}
\nu\left(\psi, \theta\right) = \frac{{\boldsymbol{B}}\cdot\nabla\zeta}{{\boldsymbol{B}}\cdot\nabla\theta} =
\frac{{B_{\text{tor}}}{h_\theta}}{{B_{\text{pol}}}R} = \frac{\left(F/R\right)h_\theta}{{B_{\text{pol}}}R} = FJ/R^2\end{aligned}\]
where \(F={B_{\text{tor}}}R\) is a function only of \(\psi\)
(sometimes called the poloidal current function).
The coordinate system is chosen so that \(x\) increases radially
outwards, from plasma to the wall. The sign of the toroidal field
\({B_{\text{tor}}}\) can then be either + or -.
The contravariant basis vectors are therefore
\[\begin{aligned}
\nabla x = {\sigma_{B\theta}}\nabla \psi \qquad \nabla y = \nabla \theta \qquad \nabla z =
{\sigma_{B\theta}}\left(\nabla\zeta - \left[\int_{\theta_0}^\theta{\frac{\partial \nu\left(\psi,
\theta\right)}{\partial \psi}} d\theta\right] \nabla\psi - \nu\left(\psi, \theta\right)\nabla\theta\right)\end{aligned}\]
The term in square brackets is the integrated local shear:
\[\begin{aligned}
I = \int_{y_0}^y\frac{\partial\nu\left(x, y\right)}{\partial\psi}dy\end{aligned}\]
Magnetic field
Magnetic field is given in Clebsch form by:
\[\begin{aligned}
{\boldsymbol{B}}= \nabla z\times \nabla x = \frac{1}{J}{\boldsymbol{e}}_y\end{aligned}\]
The contravariant components of this are then
\[\begin{aligned}
B^y = \frac{{B_{\text{pol}}}}{{h_\theta}} \qquad B^x = B^z = 0\end{aligned}\]
i.e. \({\boldsymbol{B}}\) can be written as
\[\begin{aligned}
{\boldsymbol{B}}= \frac{{B_{\text{pol}}}}{{h_\theta}}{\boldsymbol{e}}_y\end{aligned}\]
and the covariant components calculated using \(g_{ij}\) as
\[\begin{aligned}
B_x = {\sigma_{B\theta}}{B_{\text{tor}}}I R \qquad B_y = \frac{B^2 {h_\theta}}{{B_{\text{pol}}}} \qquad B_z = {\sigma_{B\theta}}{B_{\text{tor}}}R\end{aligned}\]
The unit vector in the direction of equilibrium \({\boldsymbol{B}}\) is
therefore
\[\begin{aligned}
{\boldsymbol{b}} = \frac{1}{JB}{\boldsymbol{e}}_y = \frac{1}{JB}\left[g_{xy}\nabla x + g_{yy}\nabla y
+ g_{yz}\nabla z\right]\end{aligned}\]
Jacobian and metric tensors
The Jacobian of this coordinate system is
\[\begin{aligned}
J^{-1} \equiv \left(\nabla x\times\nabla y\right)\cdot\nabla z = {B_{\text{pol}}}/ {h_\theta}\end{aligned}\]
which can be either positive or negative, depending on the sign of
\({B_{\text{pol}}}\). The contravariant metric tensor is
given by:
\[\begin{split}\begin{aligned}
g^{ij} \equiv {\boldsymbol{e}}^i \cdot{\boldsymbol{e}}^j \equiv \nabla u^i \cdot \nabla u^j = \left(%
\begin{array}{ccc}
\left(R{B_{\text{pol}}}\right)^2 & 0 & -I\left(R{B_{\text{pol}}}\right)^2 \\
0 & 1 / {h_\theta}^2 & -{\sigma_{B\theta}}\nu / {h_\theta}^2 \\
-I\left(R{B_{\text{pol}}}\right)^2 & -{\sigma_{B\theta}}\nu / {h_\theta}^2 & I^2\left(R{B_{\text{pol}}}\right)^2 + B^2 /
\left(R{B_{\text{pol}}}\right)^2
\end{array}
%
\right)\end{aligned}\end{split}\]
and the covariant metric tensor:
\[\begin{split}\begin{aligned}
g_{ij} \equiv {\boldsymbol{e}}_i \cdot{\boldsymbol{e}}_j = \left(%
\begin{array}{ccc}
I^2 R^2 + 1 / {\left({R{B_{\text{pol}}}}\right)^2}& {\sigma_{B\theta}}{B_{\text{tor}}}{h_\theta}I R / {B_{\text{pol}}}& I R^2 \\
{\sigma_{B\theta}}{B_{\text{tor}}}{h_\theta}I R / {B_{\text{pol}}}& B^2{h_\theta}^2 / {B_{\text{pol}}}^2 & {\sigma_{B\theta}}{B_{\text{tor}}}{h_\theta}R / {B_{\text{pol}}}\\
I R^2 & {\sigma_{B\theta}}{B_{\text{tor}}}{h_\theta}R / {B_{\text{pol}}}& R^2
\end{array}
%
\right)\end{aligned}\end{split}\]
Differential operators
The derivative of a scalar field \(f\) along the unperturbed
magnetic field \({\boldsymbol{b}}_0\) is given by
\[\begin{aligned}
\partial^0_{||}f \equiv {\boldsymbol{b}}_0 \cdot\nabla f =
\frac{1}{\sqrt{g_{yy}}}{\frac{\partial f}{\partial y}} = \frac{{B_{\text{pol}}}}{B{h_\theta}}{\frac{\partial f}{\partial y}}\end{aligned}\]
whilst the parallel divergence is given by
\[\begin{aligned}
\nabla^0_{||}f = B_0\partial^0_{||}\left(\frac{f}{B_0}\right)\end{aligned}\]
Using equation (25),
the Laplacian operator is given by
\[\begin{split}\begin{aligned}
\nabla^2 = &\frac{\partial^2}{\partial x^2}\left|\nabla x\right|^2 +
\frac{\partial^2}{\partial y^2}\left|\nabla y\right|^2 +
\frac{\partial^2}{\partial z^2}\left|\nabla z\right|^2 \nonumber \\
&-2\frac{\partial^2}{\partial x\partial z}I\left(R{B_{\text{pol}}}\right)^2 -
2\frac{\partial^2}{\partial y\partial z}\frac{\nu}{h_\theta^2}\\
&+\frac{\partial}{\partial x}\nabla^2x + \frac{\partial}{\partial
y}\nabla^2y + \frac{\partial}{\partial z}\nabla^2z \nonumber\end{aligned}\end{split}\]
Using equation (24) for
\(\nabla^2x = G^x\) etc, the values are
\[\begin{aligned}
\nabla^2x = \frac{{B_{\text{pol}}}}{h_\theta}\frac{\partial}{\partial x}\left(h_\theta
R^2{B_{\text{pol}}}\right) \qquad \nabla^2y = \frac{{B_{\text{pol}}}}{h_\theta}\frac{\partial}{\partial
y}\left(\frac{1}{{B_{\text{pol}}}h_\theta}\right)\end{aligned}\]
\[\begin{aligned}
\nabla^2z = -\frac{{B_{\text{pol}}}}{h_\theta}\left[\frac{\partial}{\partial x}\left(IR^2{B_{\text{pol}}}
h_\theta\right) + \frac{\partial}{\partial y}\left(\frac{\nu}{{B_{\text{pol}}}h_\theta}\right)\right]\end{aligned}\]
Neglecting some parallel derivative terms, the perpendicular Laplacian
can be written:
\[\begin{aligned}
\nabla_\perp^2= {\left({R{B_{\text{pol}}}}\right)^2}\left[{\frac{\partial^2 }{\partial {x}^2}} - 2I\frac{\partial^2}{\partial z\partial x} +
\left(I^2 + \frac{B^2}{\left({R{B_{\text{pol}}}}\right)^4}\right){\frac{\partial^2 }{\partial {z}^2}}\right] + \nabla^2 x {\frac{\partial }{\partial x}} +
\nabla^2 z{\frac{\partial }{\partial z}}\end{aligned}\]
The second derivative along the equilibrium field
\[\begin{aligned}
\partial^2_{||}\phi = \partial^0_{||}\left(\partial^0_{||}\phi\right) =
\frac{1}{\sqrt{g_{yy}}}{\frac{\partial }{\partial y}}\left(\frac{1}{\sqrt{g_{yy}}}\right){\frac{\partial \phi}{\partial y}}
+ \frac{1}{g_{yy}}\frac{\partial^2\phi}{\partial y^2}\end{aligned}\]
A common expression (the Poisson bracket in reduced MHD) is (from
equation (29))):
\[\begin{aligned}
{\boldsymbol{b}}_0\cdot\nabla\phi\times\nabla A =
\frac{1}{J\sqrt{g_{yy}}}\left[\left(g_{yy}{\frac{\partial \phi}{\partial z}} -
g_{yz}{\frac{\partial \phi}{\partial y}}\right){\frac{\partial A}{\partial x}} + \left(g_{yz}{\frac{\partial \phi}{\partial x}} -
g_{xy}{\frac{\partial \phi}{\partial z}}\right){\frac{\partial A}{\partial y}} + \left(g_{xy}{\frac{\partial \phi}{\partial y}} -
g_{yy}{\frac{\partial \phi}{\partial x}}\right){\frac{\partial A}{\partial z}}\right]\end{aligned}\]
The perpendicular nabla operator:
\[\begin{split}\begin{aligned}
\nabla_\perp \equiv& \nabla - {\boldsymbol{b}}\left({\boldsymbol{b}}\cdot\nabla\right) \\ =& \nabla
x\left({\frac{\partial }{\partial x}} - \frac{g_{xy}}{\left(JB\right)^2}{\frac{\partial }{\partial y}}\right) + \nabla
z\left({\frac{\partial }{\partial z}} - \frac{g_{yz}}{\left(JB\right)^2}{\frac{\partial }{\partial y}}\right)\end{aligned}\end{split}\]
J x B in field-aligned coordinates
Components of the magnetic field in field-aligned coordinates:
\[\begin{aligned}
B^y = \frac{{B_{\text{pol}}}}{{h_\theta}} \qquad B^x = B^z = 0\end{aligned}\]
and
\[\begin{aligned}
B_x = {\sigma_{B\theta}}{B_{\text{tor}}}I R \qquad B_y = \frac{B^2{h_\theta}}{{B_{\text{pol}}}} \qquad B_z = {\sigma_{B\theta}}{B_{\text{tor}}}R\end{aligned}\]
Calculate current \({\boldsymbol{J}}= \frac{1}{\mu}{\nabla\times
{\boldsymbol{B}} }\)
\[\begin{aligned}
\left({\nabla\times {\boldsymbol{B}} }\right)^x = \frac{1}{J}\left({\frac{\partial B_z}{\partial y}} - {\frac{\partial B_y}{\partial z}}\right) = 0\end{aligned}\]
since \({B_{\text{tor}}}R\) is a flux-surface quantity, and
\({\boldsymbol{B}}\) is axisymmetric.
\[\begin{split}\begin{aligned}
\left({\nabla\times {\boldsymbol{B}} }\right)^y =& -{\sigma_{B\theta}}\frac{{B_{\text{pol}}}}{{h_\theta}}{\frac{\partial }{\partial x}}\left({B_{\text{tor}}}R\right) \\
\left({\nabla\times {\boldsymbol{B}} }\right)^z =&
\frac{{B_{\text{pol}}}}{{h_\theta}}\left[{\frac{\partial }{\partial x}}\left(\frac{B^2{h_\theta}}{{B_{\text{pol}}}}\right) -
{\sigma_{B\theta}}{\frac{\partial }{\partial y}}\left({B_{\text{tor}}}I R\right)\right]\end{aligned}\end{split}\]
The second term can be simplified, again using
\({B_{\text{tor}}}R\) constant on flux-surfaces:
\[\begin{aligned}
{\frac{\partial }{\partial y}}\left({B_{\text{tor}}}I R\right) = {\sigma_{B\theta}}{B_{\text{tor}}}R{\frac{\partial \nu}{\partial x}} \qquad \nu =
\frac{{h_\theta}{B_{\text{tor}}}}{R{B_{\text{pol}}}}\end{aligned}\]
From these, calculate covariant components:
(12)\[\begin{split}\begin{aligned}
\left({\nabla\times {\boldsymbol{B}} }\right)_x =& -{B_{\text{tor}}}I R {\frac{\partial }{\partial x}}\left({B_{\text{tor}}}R\right) +
\frac{IR^2{B_{\text{pol}}}}{{h_\theta}}\left[{\frac{\partial }{\partial x}}\left(\frac{B^2{h_\theta}}{{B_{\text{pol}}}}\right) - {B_{\text{tor}}}
R{\frac{\partial \nu}{\partial x}}\right] \nonumber\\
%
\left({\nabla\times {\boldsymbol{B}} }\right)_y =& -{\sigma_{B\theta}}\frac{B^2{h_\theta}}{{B_{\text{pol}}}}{\frac{\partial }{\partial x}}\left({B_{\text{tor}}}R\right) +
{\sigma_{B\theta}}{B_{\text{tor}}}R\left[{\frac{\partial }{\partial x}}\left(\frac{B^2{h_\theta}}{{B_{\text{pol}}}}\right) - {B_{\text{tor}}}R{\frac{\partial \nu}{\partial x}}\right]
\\
%
\left({\nabla\times {\boldsymbol{B}} }\right)_z =& -{B_{\text{tor}}}R{\frac{\partial }{\partial x}}\left({B_{\text{tor}}}R\right) +
\frac{R^2{B_{\text{pol}}}}{{h_\theta}}\left[{\frac{\partial }{\partial x}}\left(\frac{B^2{h_\theta}}{{B_{\text{pol}}}}\right) - {B_{\text{tor}}}
R{\frac{\partial \nu}{\partial x}}\right] \nonumber\end{aligned}\end{split}\]
Calculate \({\boldsymbol{J}}\times{\boldsymbol{B}}\) using
\[\begin{aligned}
{\boldsymbol{e}}^i = \frac{1}{J}\left({\boldsymbol{e}}_j \times {\boldsymbol{e}}_k\right) \qquad {\boldsymbol{e}}_i =
J\left({\boldsymbol{e}}^j \times {\boldsymbol{e}}^k\right) \qquad i,j,k \texttt{ cyc } 1,2,3\end{aligned}\]
gives
\[\begin{split}\begin{aligned}
\mu_0 \left({\boldsymbol{J}}\times{\boldsymbol{B}}\right)^x =& \frac{1}{J}\left[\left({\nabla\times {\boldsymbol{B}} }\right)_y B_z -
\left({\nabla\times {\boldsymbol{B}} }\right)_z B_y \right]\\ =& -\frac{{B_{\text{pol}}}^3
R^2}{{h_\theta}}\left[{\frac{\partial }{\partial x}}\left(\frac{B^2{h_\theta}}{{B_{\text{pol}}}}\right) - {B_{\text{tor}}}R{\frac{\partial \nu}{\partial x}}\right]\end{aligned}\end{split}\]
Covariant components of \(\nabla P\):
\[\begin{aligned}
\left(\nabla P\right)_x = {\frac{\partial P}{\partial x}} \qquad \left(\nabla P\right)_y = \left(\nabla P\right)_z = 0\end{aligned}\]
and contravariant:
\[\begin{aligned}
\left(\nabla P\right)^x = {\left({R{B_{\text{pol}}}}\right)^2}{\frac{\partial P}{\partial x}} \qquad \left(\nabla P\right)^y = 0 \qquad
\left(\nabla P\right)^z = -I{\left({R{B_{\text{pol}}}}\right)^2}{\frac{\partial P}{\partial x}}\end{aligned}\]
Hence equating contravariant x components of
\({\boldsymbol{J}}\times{\boldsymbol{B}}= \nabla P\),
(13)\[\begin{aligned}
{\frac{\partial }{\partial x}}\left(\frac{B^2{h_\theta}}{{B_{\text{pol}}}}\right) - {B_{\text{tor}}}
R{\frac{\partial }{\partial x}}\left(\frac{{B_{\text{tor}}}{h_\theta}}{R{B_{\text{pol}}}}\right) + \frac{\mu_0{h_\theta}}{{B_{\text{pol}}}}{\frac{\partial P}{\partial x}} =
0
\end{aligned}\]
Use this to calculate \({h_\theta}\) profiles (need to fix
\({h_\theta}\) at one radial location).
Close to x-points, the above expression becomes singular, so a better
way to write it is:
\[\begin{aligned}
{\frac{\partial }{\partial x}}\left(B^2{h_\theta}\right) - {h_\theta}{B_{\text{pol}}}{\frac{\partial {B_{\text{pol}}}}{\partial x}} - {B_{\text{tor}}}
R{\frac{\partial }{\partial x}}\left(\frac{{B_{\text{tor}}}{h_\theta}}{R}\right) + \mu_0{h_\theta}{\frac{\partial P}{\partial x}} = 0\end{aligned}\]
For solving force-balance by adjusting \(P\) and \(f\)
profiles, the form used is
\[\begin{aligned}
{B_{\text{tor}}}{h_\theta}{\frac{\partial {B_{\text{tor}}}}{\partial x}} + \frac{{B_{\text{tor}}}^2{h_\theta}}{R}{\frac{\partial R}{\partial x}} +
\mu_0{h_\theta}{\frac{\partial P}{\partial x}} = -{B_{\text{pol}}}{\frac{\partial }{\partial x}}\left({B_{\text{pol}}}{h_\theta}\right)\end{aligned}\]
A quick way to calculate f is to rearrange this to:
\[\begin{aligned}
{\frac{\partial {B_{\text{tor}}}}{\partial x}} = {B_{\text{tor}}}\left[-\frac{1}{R}{\frac{\partial R}{\partial x}}\right] +
\frac{1}{{B_{\text{tor}}}}\left[-\mu_0{\frac{\partial P}{\partial x}} -
{\frac{\partial {B_{\text{pol}}}}{\partial {h_\theta}}}{\frac{\partial }{\partial x}}\left({B_{\text{pol}}}{h_\theta}\right)\right]\end{aligned}\]
and then integrate this using LSODE.
Parallel current
\[\begin{aligned}
J_{||} = {\boldsymbol{b}}\cdot{\boldsymbol{J}}\qquad b^y = \frac{{B_{\text{pol}}}}{B{h_\theta}}\end{aligned}\]
and from equation (12):
\[\begin{aligned}
J_y = \frac{{\sigma_{B\theta}}}{\mu_0}\left\{-\frac{B^2{h_\theta}}{{B_{\text{pol}}}}{\frac{\partial }{\partial x}}\left({B_{\text{tor}}}R\right) + {B_{\text{tor}}}
R\left[{\frac{\partial }{\partial x}}\left(\frac{B^2{h_\theta}}{{B_{\text{pol}}}}\right) - {\sigma_{B\theta}}{B_{\text{tor}}}R{\frac{\partial \nu}{\partial x}}\right]\right\}\end{aligned}\]
since \(J_{||} = b^yJ_y\),
\[\begin{aligned}
\mu_0 J_{||} ={\sigma_{B\theta}}\frac{{B_{\text{pol}}}{B_{\text{tor}}}
R}{B{h_\theta}}\left[{\frac{\partial }{\partial x}}\left(\frac{B^2{h_\theta}}{{B_{\text{pol}}}}\right) - {B_{\text{tor}}}R{\frac{\partial \nu}{\partial x}}\right] -
{\sigma_{B\theta}}B{\frac{\partial }{\partial x}}\left({B_{\text{tor}}}R\right)\end{aligned}\]
Curvature
For reduced MHD, need to calculate curvature term
\({\boldsymbol{b}}\times{\boldsymbol{\kappa}}\), where
\({\boldsymbol{\kappa}} =
\left({\boldsymbol{b}}\cdot\nabla\right){\boldsymbol{b}}=
-{\boldsymbol{b}}\times\left(\nabla\times{\boldsymbol{b}}\right)\).
Re-arranging, this becomes:
\[\begin{aligned}
{\boldsymbol{b}}\times{\boldsymbol{\kappa}} = \nabla\times{\boldsymbol{b}}-
{\boldsymbol{b}}\left({\boldsymbol{b}}\cdot\left(\nabla\times{\boldsymbol{b}}\right)\right)\end{aligned}\]
Components of \(\nabla\times{\boldsymbol{b}}\) are:
\[\begin{split}\begin{aligned}
\left(\nabla\times{\boldsymbol{b}}\right)^x =& {\sigma_{B\theta}}\frac{{B_{\text{pol}}}}{{h_\theta}}{\frac{\partial }{\partial y}}\left(\frac{{B_{\text{tor}}}
R}{B}\right) \\ \left(\nabla\times{\boldsymbol{b}}\right)^y =&
-{\sigma_{B\theta}}\frac{{B_{\text{pol}}}}{{h_\theta}}{\frac{\partial }{\partial x}}\left(\frac{{B_{\text{tor}}}R}{B}\right) \\
\left(\nabla\times{\boldsymbol{b}}\right)^z =&
\frac{{B_{\text{pol}}}}{{h_\theta}}{\frac{\partial }{\partial x}}\left(\frac{B{h_\theta}}{{B_{\text{pol}}}}\right) - {\sigma_{B\theta}}\frac{{B_{\text{pol}}}{B_{\text{tor}}}
R}{{h_\theta}B}{\frac{\partial \nu}{\partial x}} - {\sigma_{B\theta}}I\frac{{B_{\text{pol}}}}{{h_\theta}}{\frac{\partial }{\partial y}}\left(\frac{{B_{\text{tor}}}
R}{B}\right) \\\end{aligned}\end{split}\]
giving:
(14)\[\begin{split}\begin{aligned}
{\boldsymbol{\kappa}} =& -\frac{{B_{\text{pol}}}}{B h_\theta}\left[{\frac{\partial }{\partial x}}\left(\frac{B
h_\theta}{{B_{\text{pol}}}}\right) - {\sigma_{B\theta}}{\frac{\partial }{\partial y}}\left(\frac{{B_{\text{tor}}}I R}{B}\right)\right]\nabla x \nonumber
\\ &+ {\sigma_{B\theta}}\frac{{B_{\text{pol}}}}{B h_\theta}{\frac{\partial }{\partial y}}\left(\frac{{B_{\text{tor}}}R}{B}\right)\nabla z
\end{aligned}\end{split}\]
\[\begin{aligned}
{\boldsymbol{b}}\cdot\left(\nabla\times{\boldsymbol{b}}\right) = -{\sigma_{B\theta}}B{\frac{\partial }{\partial x}}\left(\frac{{B_{\text{tor}}}R}{B}\right) +
{\sigma_{B\theta}}\frac{{B_{\text{tor}}}{B_{\text{pol}}}R}{B{h_\theta}}{\frac{\partial }{\partial x}}\left(\frac{B{h_\theta}}{{B_{\text{pol}}}}\right) -
\frac{{B_{\text{pol}}}{B_{\text{tor}}}^2R^2}{{h_\theta}B^2}{\frac{\partial \nu}{\partial x}}\end{aligned}\]
therefore,
\[\begin{split}\begin{aligned}
\left({\boldsymbol{b}}\times{\boldsymbol{\kappa}}\right)^x =& {\sigma_{B\theta}}\frac{{B_{\text{pol}}}}{{h_\theta}}{\frac{\partial }{\partial y}}\left(\frac{{B_{\text{tor}}}
R}{B}\right) = -{\sigma_{B\theta}}\frac{{B_{\text{pol}}}{B_{\text{tor}}}R}{{h_\theta}B^2}{\frac{\partial B}{\partial y}} \\
\left({\boldsymbol{b}}\times{\boldsymbol{\kappa}}\right)^y =& \frac{{B_{\text{pol}}}^2{B_{\text{tor}}}^2
R^2}{B^3{h_\theta}^2}{\frac{\partial \nu}{\partial x}} - {\sigma_{B\theta}}\frac{{B_{\text{pol}}}^2{B_{\text{tor}}}
R}{B^2{h_\theta}^2}{\frac{\partial }{\partial x}}\left(\frac{B{h_\theta}}{{B_{\text{pol}}}}\right) \\
\left({\boldsymbol{b}}\times{\boldsymbol{\kappa}}\right)^z =&
\frac{{B_{\text{pol}}}}{{h_\theta}}{\frac{\partial }{\partial x}}\left(\frac{B{h_\theta}}{{B_{\text{pol}}}}\right) - {\sigma_{B\theta}}\frac{{B_{\text{pol}}}{B_{\text{tor}}}
R}{{h_\theta}B}{\frac{\partial \nu}{\partial x}} - I\left({\boldsymbol{b}}\times{\boldsymbol{\kappa}}\right)^x\end{aligned}\end{split}\]
Using equation (13):
\[\begin{aligned}
B{\frac{\partial }{\partial x}}\left(\frac{B{h_\theta}}{{B_{\text{pol}}}}\right) + \frac{B{h_\theta}}{{B_{\text{pol}}}}{\frac{\partial B}{\partial x}} - {\sigma_{B\theta}}{B_{\text{tor}}}
R{\frac{\partial }{\partial x}}\left(\frac{{B_{\text{tor}}}{h_\theta}}{R{B_{\text{pol}}}}\right) + \frac{\mu_0{h_\theta}}{{B_{\text{pol}}}}{\frac{\partial P}{\partial x}} =
0\end{aligned}\]
we can re-write the above components as:
\[\begin{split}\begin{aligned}
\left({\boldsymbol{b}}\times{\boldsymbol{\kappa}}\right)^y =& {\sigma_{B\theta}}\frac{{B_{\text{pol}}}{B_{\text{tor}}}
R}{B^2{h_\theta}}\left[\frac{\mu_0}{B}{\frac{\partial P}{\partial x}} + {\frac{\partial B}{\partial x}}\right] \\
\left({\boldsymbol{b}}\times{\boldsymbol{\kappa}}\right)^z =& -\frac{\mu_0}{B}{\frac{\partial P}{\partial x}} - {\frac{\partial B}{\partial x}} -
I\left({\boldsymbol{b}}\times{\boldsymbol{\kappa}}\right)^x\end{aligned}\end{split}\]
Curvature from div (b/B)
The vector \({\boldsymbol{b}}\times{\boldsymbol{\kappa}}\) is an
approximation of
\[\begin{aligned}
\frac{B}{2}\nabla\times\left(\frac{{\boldsymbol{b}}}{B}\right) \simeq {\boldsymbol{b}}\times{\boldsymbol{\kappa}}\end{aligned}\]
so can just derive from the original expression. Using the
contravariant components of \({\boldsymbol{b}}\), and the curl
operator in curvilinear coordinates (see appendix):
\[\begin{split}\begin{aligned}
\nabla\times\left(\frac{{\boldsymbol{b}}}{B}\right) =&
\frac{{B_{\text{pol}}}}{{h_\theta}}\left[\left({\frac{\partial }{\partial x}}\left(\frac{{h_\theta}}{{B_{\text{pol}}}}\right) -
{\frac{\partial }{\partial y}}\left(\frac{{\sigma_{B\theta}}{B_{\text{tor}}}IR}{B^2}\right)\right){\boldsymbol{e}}_z \right. \\ &+
{\frac{\partial }{\partial y}}\left(\frac{{\sigma_{B\theta}}{B_{\text{tor}}}R}{B^2}\right){\boldsymbol{e}}_x \\ &+
\left.{\frac{\partial }{\partial x}}\left(\frac{{\sigma_{B\theta}}{B_{\text{tor}}}R}{B^2}\right){\boldsymbol{e}}_y\right]\end{aligned}\end{split}\]
This can be simplified using
\[\begin{aligned}
{\frac{\partial }{\partial y}}\left(\frac{{\sigma_{B\theta}}{B_{\text{tor}}}IR}{B^2}\right) = I{\sigma_{B\theta}}{B_{\text{tor}}}
R{\frac{\partial }{\partial y}}\left(\frac{1}{B^2}\right) + \frac{{B_{\text{tor}}}R}{B^2}{\frac{\partial \nu}{\partial x}}\end{aligned}\]
to give
\[\begin{split}\begin{aligned}
\left({\boldsymbol{b}}\times{\boldsymbol{\kappa}}\right)^x =& -{\sigma_{B\theta}}\frac{{B_{\text{pol}}}{B_{\text{tor}}}R}{{h_\theta}B^2}{\frac{\partial B}{\partial y}} \\
\left({\boldsymbol{b}}\times{\boldsymbol{\kappa}}\right)^y =& -{\sigma_{B\theta}}\frac{B{B_{\text{pol}}}}{2{h_\theta}}{\frac{\partial }{\partial x}}\left(\frac{{B_{\text{tor}}}
R}{B^2}\right) \\ \left({\boldsymbol{b}}\times{\boldsymbol{\kappa}}\right)^z =&
\frac{B{B_{\text{pol}}}}{2{h_\theta}}{\frac{\partial }{\partial x}}\left(\frac{{h_\theta}}{{B_{\text{pol}}}}\right) - \frac{{B_{\text{pol}}}{B_{\text{tor}}}
R}{2{h_\theta}B}{\frac{\partial \nu}{\partial x}} - I\left({\boldsymbol{b}}\times{\boldsymbol{\kappa}}\cdot\nabla\right)^x\end{aligned}\end{split}\]
The first and second terms in
\(\left({\boldsymbol{b}}\times{\boldsymbol{\kappa}}\cdot\nabla\right)^z\)
almost cancel, so by expanding out \(\nu\) a better expression is
\[\begin{aligned}
\left({\boldsymbol{b}}\times{\boldsymbol{\kappa}}\right)^z = \frac{{B_{\text{pol}}}^3}{2{h_\theta}
B}{\frac{\partial }{\partial x}}\left(\frac{{h_\theta}}{{B_{\text{pol}}}}\right) - \frac{{B_{\text{tor}}}
R}{2B}{\frac{\partial }{\partial x}}\left(\frac{{h_\theta}}{{B_{\text{pol}}}}\right)\end{aligned}\]
Curvature of a single line
The curvature vector can be calculated from the field-line toroidal
coordinates \(\left(R,Z,\phi\right)\) as follows. The line element
is given by
\[\begin{aligned}
d{\boldsymbol{r}} = dR{\hat{{\boldsymbol{R}}}}+ dZ{\hat{{\boldsymbol{Z}}}}+ Rd\phi{\hat{{\boldsymbol{\phi}}}}\end{aligned}\]
Hence the tangent vector is
\[\begin{aligned}
\hat{{\boldsymbol{T}}} \equiv {\frac{d {\boldsymbol{r}}}{d s}} = {\frac{d R}{d s}}{\hat{{\boldsymbol{R}}}}+ {\frac{d Z}{d s}}{\hat{{\boldsymbol{Z}}}}+
R{\frac{d \phi}{d s}}{\hat{{\boldsymbol{\phi}}}}\end{aligned}\]
where \(s\) is the distance along the field-line. From this, the
curvature vector is given by
\[\begin{split}\begin{aligned}
{\boldsymbol{\kappa}}\equiv {\frac{d {\boldsymbol{T}}}{d s}} =& {\frac{d^2 R}{d s^2}}{\hat{{\boldsymbol{R}}}}+ {\frac{d R}{d s}}{\frac{d \phi}{d s}}{\hat{{\boldsymbol{\phi}}}}
\\ &+ {\frac{d^2 Z}{d s^2}}{\hat{{\boldsymbol{Z}}}}\\ &+ {\frac{d R}{d s}}{\frac{d \phi}{d s}}{\hat{{\boldsymbol{\phi}}}}+
R{\frac{d^2 \phi}{d s^2}}{\hat{{\boldsymbol{\phi}}}}- R\left({\frac{d \phi}{d s}}\right)^2 {\hat{{\boldsymbol{R}}}}\end{aligned}\end{split}\]
i.e.
(15)\[\begin{aligned}
{\boldsymbol{\kappa}}= \left[{\frac{d^2 R}{d s^2}} - R\left({\frac{d \phi}{d s}}\right)^2\right]{\hat{{\boldsymbol{R}}}}+ {\frac{d^2 Z}{d s^2}}{\hat{{\boldsymbol{Z}}}}+
\left[2{\frac{d R}{d s}}{\frac{d \phi}{d s}} + R{\frac{d^2 \phi}{d s^2}}\right]{\hat{{\boldsymbol{\phi}}}}
\end{aligned}\]
Want the components of
\({\boldsymbol{b}}\times{\boldsymbol{\kappa}}\),
and since the vector \({\boldsymbol{b}}\) is just the
tangent vector \({\boldsymbol{T}}\) above, this can be
written using the cross-products
\[\begin{aligned}
{\hat{{\boldsymbol{R}}}}\times{\hat{{\boldsymbol{Z}}}}= -{\hat{{\boldsymbol{\phi}}}}\qquad {\hat{{\boldsymbol{\phi}}}}\times{\hat{{\boldsymbol{Z}}}}= {\hat{{\boldsymbol{R}}}}\qquad
{\hat{{\boldsymbol{R}}}}\times{\hat{{\boldsymbol{\phi}}}}= {\hat{{\boldsymbol{Z}}}}\end{aligned}\]
This vector must then be dotted with \(\nabla\psi\),
\(\nabla\theta\), and \(\nabla\phi\). This is done by writing
these vectors in cylindrical coordinates:
\[\begin{split}\begin{aligned}
\nabla\psi =& {\frac{\partial \psi}{\partial R}}\hat{{\boldsymbol{R}}} + {\frac{\partial \psi}{\partial Z}}\hat{{\boldsymbol{Z}}} \\ \nabla\theta =&
\frac{1}{{B_{\text{pol}}}{h_\theta}}\nabla\phi\times\nabla\psi =
\frac{1}{R{B_{\text{pol}}}{h_\theta}}\left({\frac{\partial \psi}{\partial Z}}\hat{{\boldsymbol{R}}} - {\frac{\partial \psi}{\partial R}}\hat{{\boldsymbol{Z}}}\right) \\\end{aligned}\end{split}\]
An alternative is to use
\[\begin{aligned}
{\boldsymbol{b}}\times \nabla\phi = \frac{{\sigma_{B\theta}}}{BR^2}\nabla\psi\end{aligned}\]
and that the tangent vector \({\boldsymbol{T}} =
{\boldsymbol{b}}\). This gives
(16)\[\begin{aligned}
\nabla\psi = {\sigma_{B\theta}}BR\left[\frac{dR}{ds}{\boldsymbol{Z}} - \frac{dZ}{ds}{\boldsymbol{R}}\right]
\end{aligned}\]
and so because
\(d\phi / ds = {B_{\text{tor}}}/ \left(RB\right)\)
(17)\[\begin{aligned}
{\boldsymbol{\kappa}}\cdot\nabla\psi = {\sigma_{B\theta}}BR\left[ \left( \frac{{B_{\text{tor}}}^2}{RB^2} -
{\frac{d^2 R}{d s^2}}\right){\frac{d Z}{d s}} + {\frac{d^2 Z}{d s^2}}\frac{dR}{ds} \right]
\end{aligned}\]
Taking the cross-product of the tangent vector with the curvature in
equation (15) above gives
\[\begin{split}\begin{aligned}
{\boldsymbol{b}}\times{\boldsymbol{\kappa}}=& \left[\frac{{B_{\text{tor}}}}{B}{\frac{d^2 Z}{d s^2}} -
{\frac{d Z}{d s}}\left(2{\frac{d R}{d s}}{\frac{d \phi}{d s}} + R{\frac{d^2 \phi}{d s^2}}\right)\right]{\boldsymbol{R}} \\ &+
\left[{\frac{d R}{d s}}\left(2{\frac{d R}{d s}}{\frac{d \phi}{d s}} + R{\frac{d^2 \phi}{d s^2}}\right) -
\frac{{B_{\text{tor}}}}{B}\left({\frac{d^2 R}{d s^2}} - R\left({\frac{d \phi}{d s}}\right)^2\right)\right]{\boldsymbol{Z}} \\ &+
\left[{\frac{d Z}{d s}}\left({\frac{d^2 R}{d s^2}} - R\left({\frac{d \phi}{d s}}\right)^2\right) -
{\frac{d R}{d s}}{\frac{d^2 Z}{d s^2}}\right]{\hat{{\boldsymbol{\phi}}}}\end{aligned}\end{split}\]
The components in field-aligned coordinates can then be calculated:
\[\begin{split}\begin{aligned}
\left({\boldsymbol{b}}\times{\boldsymbol{\kappa}}\right)^x =& {\sigma_{B\theta}}\left({\boldsymbol{b}}\times{\boldsymbol{\kappa}}\right)\cdot\nabla\psi \\ =&
\frac{R{B_{\text{pol}}}^2}{B}\left(2{\frac{d R}{d s}}{\frac{d \phi}{d s}} + R{\frac{d^2 \phi}{d s^2}}\right) -
R{B_{\text{tor}}}\left({\frac{d R}{d s}}{\frac{d^2 R}{d s^2}} + {\frac{d Z}{d s}}{\frac{d^2 Z}{d s^2}}\right) +
\frac{{B_{\text{tor}}}^3}{B^2}{\frac{d R}{d s}}\end{aligned}\end{split}\]
Curvature in toroidal coordinates
In toroidal coordinates \(\left(\psi,\theta,\phi\right)\), the
\({\boldsymbol{b}}\) vector is
\[\begin{split}\begin{aligned}
{\boldsymbol{b}}=& \frac{{B_{\text{pol}}}}{B}{\hat{{\boldsymbol{e}}}}_\theta + \frac{{B_{\text{tor}}}}{B}{\hat{{\boldsymbol{e}}}}_\phi \\ =&
\frac{{B_{\text{pol}}}{h_\theta}}{B}\nabla\theta + \frac{R{B_{\text{tor}}}}{B}\nabla\phi\end{aligned}\end{split}\]
The curl of this vector is
\[\begin{split}\begin{aligned}
\left(\nabla\times{\boldsymbol{b}}\right)^\psi =& \frac{1}{\sqrt{g}}\left({\frac{\partial b_\phi}{\partial \theta}} -
{\frac{\partial b_\theta}{\partial \phi}}\right) \\ \left(\nabla\times{\boldsymbol{b}}\right)^\theta =&
\frac{1}{\sqrt{g}}\left({\frac{\partial b_\psi}{\partial \phi}} - {\frac{\partial b_\phi}{\partial \psi}}\right) \\
\left(\nabla\times{\boldsymbol{b}}\right)^\phi =& \frac{1}{\sqrt{g}}\left({\frac{\partial b_\theta}{\partial \psi}}
- {\frac{\partial b_\psi}{\partial \theta}}\right)\end{aligned}\end{split}\]
where
\(1/\sqrt{g} = {B_{\text{pol}}}/{h_\theta}\).
Therefore, in terms of unit vectors:
\[\begin{aligned}
\nabla\times{\boldsymbol{b}}=
\frac{1}{R{h_\theta}}{\frac{\partial }{\partial \theta}}\left(\frac{R{B_{\text{tor}}}}{B}\right){\hat{{\boldsymbol{e}}}}_\psi -
{B_{\text{pol}}}{\frac{\partial }{\partial \psi}}\left(\frac{R{B_{\text{tor}}}}{B}\right){\hat{{\boldsymbol{e}}}}_\theta + \frac{{B_{\text{pol}}}
R}{{h_\theta}}{\frac{\partial }{\partial \psi}}\left(\frac{{h_\theta}{B_{\text{pol}}}}{B}\right){\hat{{\boldsymbol{e}}}}_\phi\end{aligned}\]
psi derivative of the B field
Needed to calculate magnetic shear, and one way to get the curvature.
The simplest way is to use finite differencing, but there is another way
using local derivatives (implemented using DCT).
\[\begin{aligned}
{B_{\text{pol}}}= \frac{\left|\nabla\psi\right|}{R} = \frac{1}{R}\sqrt{\left({\frac{\partial \psi}{\partial R}}\right)^2 +
\left({\frac{\partial \psi}{\partial R}}\right)^2}\end{aligned}\]
Using
\[\begin{aligned}
\nabla{B_{\text{pol}}}= {\frac{\partial {B_{\text{pol}}}}{\partial \psi}}\nabla\psi + {\frac{\partial {B_{\text{pol}}}}{\partial \theta}}\nabla\theta +
{\frac{\partial {B_{\text{pol}}}}{\partial \phi}}\nabla\phi\end{aligned}\]
we get
\[\begin{aligned}
\nabla{B_{\text{pol}}}\cdot\nabla\psi = {\frac{\partial {B_{\text{pol}}}}{\partial \psi}}\left|\nabla\psi\right|^2\end{aligned}\]
and so
\[\begin{aligned}
{\frac{\partial {B_{\text{pol}}}}{\partial \psi}} = \nabla{B_{\text{pol}}}\cdot\nabla\psi / \left(R{B_{\text{pol}}}\right)^2\end{aligned}\]
The derivatives of \({B_{\text{pol}}}\) in \(R\) and
\(Z\) are:
\[\begin{split}\begin{aligned}
{\frac{\partial {B_{\text{pol}}}}{\partial R}} =& -\frac{{B_{\text{pol}}}}{R} + \frac{1}{{B_{\text{pol}}}
R^2}\left[{\frac{\partial \psi}{\partial R}}{\frac{\partial^2 \psi}{\partial {R}^2}} +
{\frac{\partial \psi}{\partial Z}}\frac{\partial^2\psi}{\partial R\partial Z}\right] \\ {\frac{\partial {B_{\text{pol}}}}{\partial Z}}
=& \frac{1}{{B_{\text{pol}}}R^2}\left[{\frac{\partial \psi}{\partial Z}}{\frac{\partial^2 \psi}{\partial {Z}^2}} +
{\frac{\partial \psi}{\partial R}}\frac{\partial^2\psi}{\partial R\partial Z}\right]\end{aligned}\end{split}\]
For the toroidal field, \({B_{\text{tor}}}= f/R\)
\[\begin{aligned}
{\frac{\partial {B_{\text{tor}}}}{\partial \psi}} = \frac{1}{R}{\frac{\partial f}{\partial \psi}} - \frac{f}{R^2}{\frac{\partial R}{\partial \psi}}\end{aligned}\]
As above,
\({\frac{\partial R}{\partial \psi}} = \nabla R \cdot\nabla\psi / \left(R{B_{\text{pol}}}\right)^2\),
and since \(\nabla R\cdot\nabla R = 1\),
\[\begin{aligned}
{\frac{\partial R}{\partial \psi}} = {\frac{\partial \psi}{\partial R}} / \left(R{B_{\text{pol}}}\right)^2\end{aligned}\]
similarly,
\[\begin{aligned}
{\frac{\partial Z}{\partial \psi}} = {\frac{\partial \psi}{\partial Z}} / \left(R{B_{\text{pol}}}\right)^2\end{aligned}\]
and so the variation of toroidal field with \(\psi\) is
\[\begin{aligned}
{\frac{\partial {B_{\text{tor}}}}{\partial \psi}} = \frac{1}{R}{\frac{\partial f}{\partial \psi}} -
\frac{{B_{\text{tor}}}}{R^3{B_{\text{pol}}}^2}{\frac{\partial \psi}{\partial R}}\end{aligned}\]
From the definition
\(B=\sqrt{{B_{\text{tor}}}^2 + {B_{\text{pol}}}^2}\),
\[\begin{aligned}
{\frac{\partial B}{\partial \psi}} = \frac{1}{B}\left({B_{\text{tor}}}{\frac{\partial {B_{\text{tor}}}}{\partial \psi}} + {B_{\text{pol}}}{\frac{\partial {B_{\text{pol}}}}{\partial \psi}}\right)\end{aligned}\]
Parallel derivative of the B field
To get the parallel nablaients of the \(B\) field components, start
with
\[\begin{aligned}
{\frac{\partial }{\partial s}}\left(B^2\right) = {\frac{\partial }{\partial s}}\left({B_{\text{tor}}}^2\right) + {\frac{\partial }{\partial s}}\left({B_{\text{pol}}}^2\right)\end{aligned}\]
Using the fact that \(R{B_{\text{tor}}}\) is constant
along \(s\),
\[\begin{aligned}
{\frac{\partial }{\partial s}}\left(R^2{B_{\text{tor}}}^2\right) = R^2{\frac{\partial }{\partial s}}\left({B_{\text{tor}}}^2\right) +
{B_{\text{tor}}}^2{\frac{\partial }{\partial s}}\left(R^2\right) = 0\end{aligned}\]
which gives
\[\begin{aligned}
{\frac{\partial }{\partial s}}\left({B_{\text{tor}}}^2\right) = -\frac{{B_{\text{tor}}}^2}{R^2}{\frac{\partial }{\partial s}}\left(R^2\right)\end{aligned}\]
The poloidal field can be calculated from
\[\begin{aligned}
{\frac{\partial }{\partial s}}\left(\nabla\psi \cdot \nabla\psi\right) = {\frac{\partial }{\partial s}}\left(R^2{B_{\text{pol}}}^2\right) =
R^2{\frac{\partial }{\partial s}}\left({B_{\text{pol}}}^2\right) + {B_{\text{pol}}}^2{\frac{\partial }{\partial s}}\left(R^2\right)\end{aligned}\]
Using equation (16),
\(\nabla\psi \cdot \nabla\psi\) can also be written as
\[\begin{aligned}
\nabla\psi \cdot \nabla\psi = B^2R^2\left[\left({\frac{\partial R}{\partial s}}\right)^2 +
\left({\frac{\partial Z}{\partial s}}\right)^2\right]\end{aligned}\]
and so (unsurprisingly)
\[\begin{aligned}
\frac{{B_{\text{pol}}}^2}{B^2} = \left[\left({\frac{\partial R}{\partial s}}\right)^2 + \left({\frac{\partial Z}{\partial s}}\right)^2\right]\end{aligned}\]
Hence
\[\begin{aligned}
{\frac{\partial }{\partial s}}\left({B_{\text{pol}}}^2\right) = B^2{\frac{\partial }{\partial s}}\left[\left({\frac{\partial R}{\partial s}}\right)^2 +
\left({\frac{\partial Z}{\partial s}}\right)^2\right] + \frac{{B_{\text{pol}}}^2}{B^2}{\frac{\partial }{\partial s}}\left(B^2\right)\end{aligned}\]
Which gives
\[\begin{aligned}
{\frac{\partial }{\partial s}}\left(B^2\right) = -\frac{B^2}{R^2}{\frac{\partial }{\partial s}}\left(R^2\right) +
\frac{B^4}{{B_{\text{tor}}}^2}{\frac{\partial }{\partial s}}\left[\left({\frac{\partial R}{\partial s}}\right)^2 + \left({\frac{\partial Z}{\partial s}}\right)^2\right]\end{aligned}\]
\[\begin{aligned}
{\frac{\partial }{\partial s}}\left({B_{\text{pol}}}^2\right) = \left(1 +
\frac{{B_{\text{pol}}}^2}{{B_{\text{tor}}}^2}\right)B^2{\frac{\partial }{\partial s}}\left[\left({\frac{\partial R}{\partial s}}\right)^2 +
\left({\frac{\partial Z}{\partial s}}\right)^2\right] - \frac{{B_{\text{pol}}}^2}{R^2}{\frac{\partial }{\partial s}}\left(R^2\right)\end{aligned}\]
Magnetic shear from J x B
Re-arranging the radial force balance
equation (13) gives
\[\begin{aligned}
\frac{{B_{\text{pol}}}^2R}{{B_{\text{tor}}}}{\frac{\partial \nu}{\partial \psi}} + \nu\left(\frac{2RB}{{B_{\text{tor}}}}{\frac{\partial B}{\partial \psi}} +
\frac{B^2}{{B_{\text{tor}}}}{\frac{\partial R}{\partial \psi}} - \frac{B^2R}{{B_{\text{tor}}}^2}{\frac{\partial {B_{\text{tor}}}}{\partial \psi}}\right) +
\frac{\mu_0{h_\theta}}{{B_{\text{pol}}}}{\frac{\partial P}{\partial \psi}} = 0\end{aligned}\]
Magnetic shear
The field-line pitch is given by
\[\begin{aligned}
\nu = \frac{{h_\theta}{B_{\text{tor}}}}{{B_{\text{pol}}}R}\end{aligned}\]
and so
\[\begin{aligned}
{\frac{\partial \nu}{\partial \psi}} = \frac{\nu}{{h_\theta}}{\frac{\partial {h_\theta}}{\partial \psi}} +
\frac{\nu}{{B_{\text{tor}}}}{\frac{\partial {B_{\text{tor}}}}{\partial \psi}} - \frac{\nu}{{B_{\text{pol}}}}{\frac{\partial {B_{\text{pol}}}}{\partial \psi}} -
\frac{\nu}{R}{\frac{\partial R}{\partial \psi}}\end{aligned}\]
The last three terms are given in the previous section, but
\(\partial{h_\theta}/\partial\psi\) needs to be evaluated
psi derivative of h
From the expression for curvature (equation (14)),
and using
\(\nabla x \cdot \nabla \psi = {\sigma_{B\theta}}\left(R{B_{\text{pol}}}\right)^2\)
and
\(\nabla z\cdot\nabla \psi = -{\sigma_{B\theta}}I \left(R{B_{\text{pol}}}\right)^2\)
\[\begin{split}\begin{aligned}
{\boldsymbol{\kappa}}\cdot\nabla\psi =& -{\sigma_{B\theta}}
\frac{{B_{\text{pol}}}}{B{h_\theta}}{\left({R{B_{\text{pol}}}}\right)^2}\left[{\frac{\partial }{\partial x}}\left(\frac{B{h_\theta}}{{B_{\text{pol}}}}\right) -
{\sigma_{B\theta}}{\frac{\partial }{\partial y}}\left(\frac{{B_{\text{tor}}}IR}{B}\right)\right] \\ &- I{\left({R{B_{\text{pol}}}}\right)^2}
\frac{{B_{\text{pol}}}}{B{h_\theta}}{\frac{\partial }{\partial y}}\left(\frac{{B_{\text{tor}}}R}{B}\right)\end{aligned}\end{split}\]
The second and third terms partly cancel, and using
\({\frac{\partial I}{\partial y}} = {\sigma_{B\theta}}
{\frac{\partial \nu}{\partial x}}\)
\[\begin{split}\begin{aligned}
\frac{{\boldsymbol{\kappa}}\cdot\nabla\psi}{{\left({R{B_{\text{pol}}}}\right)^2}} =&
-{\sigma_{B\theta}}\frac{{B_{\text{pol}}}}{B{h_\theta}}{\frac{\partial }{\partial x}}\left(\frac{B{h_\theta}}{{B_{\text{pol}}}}\right) +
{\sigma_{B\theta}}\frac{{B_{\text{pol}}}}{B{h_\theta}}\frac{{B_{\text{tor}}}R}{B}{\frac{\partial \nu}{\partial x}} \\ =&
-{\sigma_{B\theta}}\frac{{B_{\text{pol}}}}{B{h_\theta}}\left[{\frac{\partial }{\partial x}}\left(\frac{B{h_\theta}}{{B_{\text{pol}}}}\right) - \frac{{B_{\text{tor}}}
R}{B}{\frac{\partial }{\partial x}}\left(\frac{{B_{\text{tor}}}{h_\theta}}{{B_{\text{pol}}}R}\right)\right] \\ =&
-{\sigma_{B\theta}}\frac{{B_{\text{pol}}}}{B{h_\theta}}\left[{h_\theta}{\frac{\partial }{\partial x}}\left(\frac{B}{{B_{\text{pol}}}}\right) -
{h_\theta}\frac{{B_{\text{tor}}}R}{B}{\frac{\partial }{\partial x}}\left(\frac{{B_{\text{tor}}}}{{B_{\text{pol}}}R}\right) +
\frac{B^2}{B{B_{\text{pol}}}}{\frac{\partial {h_\theta}}{\partial x}} -
\frac{{B_{\text{tor}}}^2}{B{B_{\text{pol}}}}{\frac{\partial {h_\theta}}{\partial x}}\right] \\ =& -{\sigma_{B\theta}}
\frac{{B_{\text{pol}}}}{B^2{h_\theta}}{\frac{\partial {h_\theta}}{\partial x}} -
{\sigma_{B\theta}}\frac{{B_{\text{pol}}}}{B^2}\left[B{\frac{\partial }{\partial x}}\left(\frac{B}{{B_{\text{pol}}}}\right) - {B_{\text{tor}}}
R{\frac{\partial }{\partial x}}\left(\frac{{B_{\text{tor}}}}{{B_{\text{pol}}}R}\right)\right]\end{aligned}\end{split}\]
Writing
\[\begin{split}\begin{aligned}
B{\frac{\partial }{\partial x}}\left(\frac{B}{{B_{\text{pol}}}}\right) =& {\frac{\partial }{\partial x}}\left(\frac{B^2}{{B_{\text{pol}}}}\right) -
\frac{B}{{B_{\text{pol}}}}{\frac{\partial B}{\partial x}} \\ {B_{\text{tor}}}R{\frac{\partial }{\partial x}}\left(\frac{{B_{\text{tor}}}}{{B_{\text{pol}}}R}\right) =&
{\frac{\partial }{\partial x}}\left(\frac{{B_{\text{tor}}}^2}{{B_{\text{pol}}}}\right) - \frac{{B_{\text{tor}}}}{{B_{\text{pol}}}R}{\frac{\partial }{\partial x}}\left({B_{\text{tor}}}
R\right)\end{aligned}\end{split}\]
and using
\(B{\frac{\partial B}{\partial x}} = {B_{\text{tor}}}{\frac{\partial {B_{\text{tor}}}}{\partial x}} + {B_{\text{pol}}}{\frac{\partial {B_{\text{pol}}}}{\partial x}}\),
this simplifies to give
(18)\[\begin{aligned}
\frac{{\boldsymbol{\kappa}}\cdot\nabla\psi}{{\left({R{B_{\text{pol}}}}\right)^2}} =
-{\sigma_{B\theta}}\frac{{B_{\text{pol}}}^2}{B^2{h_\theta}}{\frac{\partial {h_\theta}}{\partial x}} - {\sigma_{B\theta}}\frac{{B_{\text{tor}}}^2}{B^2
R}{\frac{\partial R}{\partial x}}
\end{aligned}\]
This can be transformed into an expression for
\({\frac{\partial {h_\theta}}{\partial x}}\)
involving only derivatives along field-lines. Writing \(\nabla R =
{\frac{\partial R}{\partial \psi}}\nabla\psi + {\frac{\partial R}{\partial \theta}}\nabla\theta\),
\[\begin{aligned}
\nabla R \cdot \nabla\psi = {\frac{\partial R}{\partial \psi}}{\left({R{B_{\text{pol}}}}\right)^2}\end{aligned}\]
Using (16),
\[\begin{aligned}
\nabla\psi \cdot \nabla R = -{\sigma_{B\theta}}B R\frac{dZ}{ds}\end{aligned}\]
and so
\[\begin{aligned}
{\frac{\partial R}{\partial x}} = -\frac{BR}{{\left({R{B_{\text{pol}}}}\right)^2}}\frac{dZ}{ds}\end{aligned}\]
Substituting this and equation (17)
for \({\boldsymbol{\kappa}}\cdot\nabla\psi\) into
equation (18) the
\({\frac{\partial R}{\partial x}}\) term cancels with
part of the \({\boldsymbol{\kappa}}\cdot\nabla\psi\)
term, simplifying to
\[\begin{aligned}
{\frac{\partial {h_\theta}}{\partial x}} =
-{h_\theta}\frac{B^3R}{{B_{\text{pol}}}^2{\left({R{B_{\text{pol}}}}\right)^2}}\left[\frac{d^2Z}{ds^2}\frac{dR}{ds} -
\frac{d^2R}{ds^2}\frac{dZ}{ds}\right]\end{aligned}\]
Differential geometry
Warning
Several mistakes have been found (and is now corrected)
in this section, so it should be proof read before removing this
warning! The following are notes from [haeseler].
Sets of vectors \(\left\{\mathbf{A, B, C}\right\}\) and
\(\left\{\mathbf{a, b, c}\right\}\) are reciprocal if
\[\begin{split}\begin{aligned}
\mathbf{A\cdot a} = \mathbf{B\cdot b} = \mathbf{C\cdot c} = 1\\ \mathbf{A\cdot
b} = \mathbf{A\cdot c} = \mathbf{B\cdot a} = \mathbf{B\cdot c} = \mathbf{C\cdot
a} = \mathbf{C\cdot b} = 0 \\\end{aligned}\end{split}\]
which implies that \(\left\{\mathbf{A, B, C}\right\}\) and
\(\left\{\mathbf{a, b, c}\right\}\) are each linearly independent.
Equivalently,
\[\begin{aligned}
\mathbf{a} = \frac{\mathbf{B\times C}}{\mathbf{A\cdot\left(B\times C\right)}}\qquad
{\boldsymbol{b}}= \frac{\mathbf{C\times A}}{\mathbf{B\cdot\left(C\times A\right)}}\qquad
\mathbf{c} = \frac{\mathbf{A\times B}}{\mathbf{C\cdot\left(A\times B\right)}}\end{aligned}\]
Either of these sets can be used as a basis, and any vector
\(\mathbf{w}\) can be represented as
\(\mathbf{w} = \left(\mathbf{w\cdot a}\right)\mathbf{A} +
\left(\mathbf{w\cdot b}\right){\boldsymbol{B}}+ \left(\mathbf{w\cdot c}\right)\mathbf{C}\)
or as
\(\mathbf{w} = \left(\mathbf{w\cdot A}\right)\mathbf{a} + \left(\mathbf{w\cdot B}\right){\boldsymbol{b}}
+ \left(\mathbf{w\cdot C}\right)\mathbf{c}\). In the Cartesian coordinate
system, the basis vectors are reciprocal to themselves so this
distinction is not needed. For a general set of coordinates
\(\left\{u^1, u^2, u^3\right\}\), tangent basis vectors can be
defined. If the Cartesian coordinates of a point are given by
\(\left(x, y, z\right) = \mathbf{R}\left(u^1, u^2, u^3\right)\) then
the tangent basis vectors are:
\[\begin{aligned}
{\boldsymbol{e}}_i = \frac{\partial\mathbf{R}}{\partial u^i}\end{aligned}\]
and in general these will vary from point to point. The scale factor or
metric coefficient
\(h_i =\left|{\boldsymbol{e}}_i\right|\) is the distance
moved for a unit change in \(u^i\). The unit vector
\(\hat{{\boldsymbol{e}}}_i = {\boldsymbol{e}}_i/h_i\).
Definition of nabla operator:
From the chain rule,
\(d\mathbf{R} = \frac{\partial\mathbf{R}}{\partial u^i}du^i
= {\boldsymbol{e}}_idu^i\) and substituting \(\Phi = u^i\)
\[\begin{aligned}
du^i = \nabla u^i\cdot{\boldsymbol{e}}_jdu^j\end{aligned}\]
which can only be true if
\(\nabla u^i\cdot{\boldsymbol{e}}_j = \delta^i_j\) i.e.
if
\[\text{Sets of vectors $\boldsymbol{e}^i\equiv\nabla u^i$ and
$\boldsymbol{e}_j$ are reciprocal}\]
Since the sets of vectors
\(\left\{{\boldsymbol{e}}^i\right\}\) and
\(\left\{{\boldsymbol{e}}_i\right\}\) are reciprocal, any
vector \(\mathbf{D}\) can be written as
\(\mathbf{D} = D_i{\boldsymbol{e}}^i
= D^i{\boldsymbol{e}}_i\) where
\(D_i = \mathbf{D\cdot e}_i\) are the covariant components and
\(D^i = \mathbf{D\cdot e}^i\) are the contravariant components. To
convert between covariant and contravariant components, define the
metric coefficients \(g_{ij} = \mathbf{e_i\cdot e_j}\) and
\(g^{ij} =
\mathbf{e^i\cdot e^j}\) so that
\({\boldsymbol{e}}_i = g_{ij}{\boldsymbol{e}}^j\).
\(g_{ij}\) and \(g^{ij}\) are symmetric and if the basis is
orthogonal then \(g_{ij}=g^{ij} = 0\) for \(i\neq j\) i.e. the
metric is diagonal.
\[\text{$g_{ij} = h_ih_j\boldsymbol{e}_i\cdot\boldsymbol{e}_j$ and so $g_{ii} = h_i^2$}\]
For a general set of coordinates, the nabla operator can be expressed as
\[\begin{aligned}
\nabla = \nabla u^i\frac{\partial}{\partial u^i} =
{\boldsymbol{e}}^i\frac{\partial}{\partial u^i}\end{aligned}\]
and for a general set of (differentiable) coordinates
\(\left\{u^i\right\}\), the Laplacian is given by
(23)\[\begin{aligned}
\nabla^2\phi = \frac{1}{J}\frac{\partial}{\partial
u^i}\left(Jg^{ij}\frac{\partial\phi}{\partial u^j}\right)
\end{aligned}\]
which can be expanded as
(24)\[\begin{aligned}
\nabla^2\phi = g^{ij}\frac{\partial^2\phi}{\partial u^i\partial u^j} +
\underbrace{\frac{1}{J}\frac{\partial}{\partial
u^i}\left(Jg^{ij}\right)}_{G^j}\frac{\partial\phi}{\partial u^j}
\end{aligned}\]
where \(G^j\) must not be mistaken as the so called connection
coefficients (i.e. the Christoffel symbols of second kind). Setting
\(\phi =
u^k\) in equation (23) gives
\(\nabla^2u^k = G^k\). Expanding
(23) and setting
\(\left\{u^i\right\} = \left\{x, y, z\right\}\) gives
(25)\[\begin{split}\begin{aligned}
\nabla^2f &= \nabla\cdot\nabla f = \nabla\cdot\left(\frac{\partial}{\partial
x}\nabla x + \frac{\partial}{\partial y}\nabla y + \frac{\partial}{\partial
z}\nabla z\right) \nonumber \\
&= \frac{\partial^2 f}{\partial x^2}\left|\nabla x\right|^2 + \frac{\partial^2
f}{\partial y^2}\left|\nabla y\right|^2 + \frac{\partial^2 f}{\partial z^2}\left|\nabla
z\right|^2 \\ +2\frac{\partial^2 f}{\partial x\partial y}\left(\nabla x\cdot\nabla
y\right) +2\frac{\partial^2 f}{\partial x\partial z}\left(\nabla x\cdot\nabla z\right)
+2\frac{\partial^2 f}{\partial y\partial z}\left(\nabla y\cdot\nabla z\right)
\nonumber \\ +\nabla^2x\frac{\partial f}{\partial x} +\nabla^2y\frac{\partial
f}{\partial y} + \nabla^2z\frac{\partial f}{\partial z} \nonumber
\end{aligned}\end{split}\]
Curl defined as:
(26)\[\begin{aligned}
\nabla\times\mathbf{A} = \frac{1}{\sqrt{g}}\sum_k\left(\frac{\partial
A_j}{\partial u_i} - \frac{\partial A_i}{\partial u_j}\right){\boldsymbol{e}}_k \qquad i,j,k
\texttt{ cyc } 1,2,3 \end{aligned}\]
Cross-product relation between contravariant and covariant vectors:
\[\begin{aligned}
{\boldsymbol{e}}^i = \frac{1}{J}\left({\boldsymbol{e}}_j \times {\boldsymbol{e}}_k\right) \qquad {\boldsymbol{e}}_i =
J\left({\boldsymbol{e}}^j \times {\boldsymbol{e}}^k\right) \qquad i,j,k \texttt{ cyc } 1,2,3\end{aligned}\]