Solution (source code)

= Solution

Setting $I=0$ in <cylindrical magnetostatic pressure balance> gives $p+B_z^2/(8\pi)=\mathrm{constant}$. The axial boundary value sets this constant to $p_0$. Thus the <magnetic pressure support with vanishing axial field> solution is
$$
\boxed{B_z(R)=\sigma\sqrt{8\pi p_0}\sqrt{1-e^{-R^2/a^2}},\qquad \sigma=\pm1.}
$$
The azimuthal component of the magnetostatic <Ampère-Maxwell equation> gives
$$
\boxed{j_\phi(R)=-\sigma\frac{c\sqrt{8\pi p_0}}{4\pi a^2}
\frac{R e^{-R^2/a^2}}{\sqrt{1-e^{-R^2/a^2}}}.}
$$
Since $1-e^{-R^2/a^2}=R^2/a^2+O(R^4)$, the radial limit is
$$
\boxed{\lim_{R\downarrow0}j_\phi(R)=-\sigma\frac{c\sqrt{8\pi p_0}}{4\pi a}.}
$$
There is a regularity subtlety: $B_z\sim\sigma\sqrt{8\pi p_0}R/a$ is not a <smooth function> of Cartesian position at the axis, and the nonzero limiting $j_\phi$ multiplies an azimuthal unit vector with no unique direction there. These formulas solve the radial problem for $R>0$ and give the requested scalar limit, but the data in this part do not admit the fully smooth vector-field regularity assumed in part (a).