Is the above $\mathbf {B}$ divergence-free?

Let us check this. Write $\mathbf {B}$ as

$$\mathbf{B}= B^{(1)} \mathcal{J} \nabla \theta \times \nabla \phi ... ...nabla \phi \times \nabla r + B^{(3)} \mathcal{J} \nabla r \times \nabla \theta,$$ (345)

where $\mathcal{J}$ is the Jacobian of $(r, \theta, \phi)$ coordinates; $B^{(1)}$, $B^{(2)}$, and $B^{(3)}$ are given by

$$B^{(1)} =\mathbf{B} \cdot \nabla r$$ (346)

$$B^{(2)} =\mathbf{B} \cdot \nabla \theta$$ (347)

$$B^{(3)} =\mathbf{B} \cdot \nabla \phi$$ (348)

Use $\mathbf{B}_p$ given by (342), then $B^{(1)}$, $B^{(2)}$, and $B^{(3)}$ are written as

$$B^{(1)} = 0,$$ (349)

$$B^{(2)} = - \frac{1}{\mathcal{J}} \frac{1}{q (r)} \frac{g_0 r}{\sqrt{R_0^2 - r^2}},$$ (350)

and

$$B^{(3)} = \frac{B_{\phi}}{R},$$ (351)

respectively. Then, by using the divergence formula in $(r, \theta, \phi)$ coordinates, $\nabla \cdot \mathbf{B}$ is written as

$$\nabla \cdot \mathbf{B} = \frac{1}{\mathcal{J}} \left( \frac{\partial B^{(1)} \mathcal{J}}{... ...}{\partial \theta} + \frac{\partial B^{(3)} \mathcal{J}}{\partial \phi} \right)$$

$$= \frac{1}{\mathcal{J}} \left( 0 + \frac{\partial B^{(2)} \mathcal{J}}{\partial \theta} + 0 \right)$$

$$= \frac{1}{\mathcal{J}} \frac{\partial}{\partial \theta} \left( - \frac{1}{q (r)} \frac{g_0 r}{\sqrt{R_0^2 - r^2}} \right)$$

$$=$$ 0 (352)

i.e., $\mathbf {B}$ in this case is indeed divergence-free.