Solution (source code)

= Solution

Axisymmetry makes all cylindrical components independent of $\phi$. The <divergence-free> condition is therefore $\partial_R(RB_R)+\partial_z(RB_z)=0$. On a simply connected meridional patch introduce a <streamfunction> $\psi$ so that
$$
\boxed{B_R=-\frac{\partial_z\psi}{R},\qquad B_z=\frac{\partial_R\psi}{R}.}
$$
These definitions satisfy the solenoidal condition identically and represent any such poloidal field locally. Since $\nabla\phi=e_\phi/R$, they give $B_p=\nabla\psi\times\nabla\phi$. The toroidal field is independent of this constraint, yielding the full representation $B=B_p+B_\phi e_\phi$. Global use of a single flux function assumes the corresponding meridional domain has no obstructing flux topology; a regular field extends to the axis by the usual limiting construction.

The <axisymmetric magnetic flux function> has a direct physical normalization. The flux through an annulus in a plane of constant $z$ is $\int2\pi RB_z\,dR=2\pi\Delta\psi$. Thus $\psi$ is poloidal <magnetic flux> divided by $2\pi$, up to an arbitrary additive constant. Also $B\cdot\nabla\psi=0$, including the toroidal part, so its level surfaces are <axisymmetric magnetic flux surfaces>. Field lines remain on those surfaces even when they wind azimuthally.