Solution (source code)

= Solution

The solenoidal constraint on the <magnetic field> reduces to $\partial_xB_x+\partial_zB_z=0$. Locally, or globally in a simply connected cross-section, it permits a <Cartesian magnetic flux function> satisfying
$$
B_x=-\psi_z,\qquad B_z=\psi_x,\qquad
\boxed{\mathbf B=\nabla\times(\psi\mathbf e_y)+B_y\mathbf e_y.}
$$
The transverse <magnetic field lines> are contours of $\psi$, because $\mathbf B\cdot\nabla\psi=0$. The independent component $B_y$ supplies the twist of the <flux tube>; it is not restricted by the solenoidal constraint.

For the <Lorentz force>, take the <curl> explicitly:
$$
\nabla\times\mathbf B
=(-B_{y,z},-\nabla^2\psi,B_{y,x}).
$$
Its <cross product> with $\mathbf B$ has transverse components $-(\nabla^2\psi)\nabla\psi-B_y\nabla B_y$ and $y$ component $\psi_xB_{y,z}-\psi_zB_{y,x}$. Consequently
$$
\boxed{\mathbf f_L=
-\frac1{\mu_0}\left[
(\nabla^2\psi)\nabla\psi+B_y\nabla B_y+
\nabla\psi\times\nabla B_y\right].}
$$
The term $-B_y\nabla B_y/\mu_0$ is the transverse gradient of the <magnetic pressure> associated with the axial field, while the other terms include <magnetic tension>.