Why the divernce of this magnetic field is not zero?

Question

I am working on a project on which I need to calculate the geomagnetic field in different coordinates. When I use the conventional form of the dipole field in spherical coordinates: $$\vec{B}_{r,\phi}=\frac{B_0 R^3}{r^3}(-2 \cos \phi\, \hat{e}_r- \sin \phi\, \hat{e}_\phi),$$ with $B_0$ as the field intensity in the equator, $R$ the radius of the Earth, $r$ is an arbitrary distance from the center from the center, $\phi$ is the polar angle. The divergence of this field is zero. $$\nabla \cdot \vec{B}_{r,\phi}=0.$$ When I transform it to cartesian coordinates $\hat{I}$,$\hat{J}$,$\hat{K}$ $$\vec{B}_{X,Y,Z}=\frac{B_0 R^3}{(X^2+Y^2+Z^2)^{5/2}}[-3 X Z\, \hat{I} -3 YZ\, \hat{J}+(X^2+Y^2-2Z^2)\hat{K}],$$ the divergence of this field is also zero $$\nabla \cdot \vec{B}_{X,Y,Z}=0$$ The transformation on which I am having trouble is happening when I convert the $\vec{B}_{X,Y,Z}$ to the cyllindrical orbital coordinates through and ordinary rotation matrix 3-1-3, wehre the first rotation is about the $\hat{K}$-axis, the longitude of ascending node $\Omega$, the second around the new $\hat{X}^\prime$-axis, the inclination $i$, and the final rotation around the $\hat{K}^\prime$-axis, is the sum of the periapsos argument $\omega$ and the true anomaly $f$, so the total angle of this final rotation is $\omega+f$, and the variables $X$, $Y$ and $Z$ are given by: $$X=r (\sin (f) (\cos (i) \cos (\omega ) \sin (\Omega )+\sin (\omega ) \cos (\Omega ))+\cos (f) (\cos (i) \sin (\omega ) \sin (\Omega )-\cos (\omega ) \cos (\Omega )));$$ $$Y=-r(\cos (i) \cos (\Omega ) \sin (f+\omega )+\sin (\Omega ) \cos (f+\omega ));$$ $$Z=-r \sin (i) \sin (f+\omega ).$$ When I do these transformations, I obtain a magnetic field in this form: $$\vec{B}_{r,f}=\frac{B_0 R^3}{r^3}[-2 \sin i \sin (\omega+f)\, \hat{r} + \sin i \cos(\omega+f)\, \hat{f} + \cos i\, \hat{k}],$$ since is a cyllindrical frame, I choosed to name the the basis as $\hat{r}$, $\hat{f}$ and $\hat{k}$, wehre the variable $r$ has the same interpretation of distance and the variable $f$, the true anomaly is the angle analogous to an angle of cyllindrical coordinates.

However when I take the divergence of $\vec{B}_{r,f}$, it is not equals zero, actually: $$\nabla \cdot \vec{B}_{r,f}=\frac{3 B_0 R^3 \sin i \sin (\omega+f)}{r^4}.$$ I know that Maxwell's equations states that the divergence of any magnetic field is zero so I am just wondering why this is happening, I am pretty sure that the transformations are correct, but I can't actually understand why the divergence of this field is non-zero.

Obs: I also made the transformation using the electromagnetic field tensor formalism for the rotation and it gives the same result as the commom vector transformation by a matrix of rotation.

Answer 1

I confirm that the divergence of $$\tag{1} \vec{B}_{X,Y,Z}=\frac{B_0 R^3}{(X^2+Y^2+Z^2)^{5/2}}[-3X Z\, \hat{I} -3YZ\, \hat{J}+(X^2+Y^2-2Z^2)\hat{K}], $$ which I prefer to write as $$\tag{2} B(x,y,z)=\frac{1}{(x^2+y^2+z^2)^{5/2}}\begin{pmatrix}-3xz\\-3yz\\x^2+y^2-2z^2\end{pmatrix} $$ is zero: \begin{align} \nabla\cdot B&=\frac{-3z-3z-4z}{(x^2+y^2+z^2)^{5/2}}-\frac{5}{2}\frac{-6x^2z-6y^2z+2z(x^2+y^2-2z^2)}{(x^2+y^2+z^2)^{7/2}}\\[2mm] &=\frac{-10z(x^2+y^2+z^2)+15x^2z+15y^2z-5x^2z-5y^2z+10z^3}{(x^2+y^2+z^2)^{7/2}}=0\,.\tag{3} \end{align} The conversion of the components of any vector field from Cartesian coordinates $(x,y,z)$ into cylindrical coordinates $(\rho,\varphi,z)$ is done as follows: \begin{align}\tag{4} \begin{pmatrix}B_\rho\\B_\varphi\\B_z\end{pmatrix}=\begin{pmatrix}\cos\varphi & \sin\varphi& 0\\-\sin\varphi &\cos\varphi&0\\0&0&1\end{pmatrix} \begin{pmatrix}B_x\\B_y\\B_z\end{pmatrix}\,. \end{align} With (2) we get \begin{align} \begin{pmatrix}B_\rho\\B_\varphi\\B_z\end{pmatrix}&= \frac{1}{(x^2+y^2+z^2)^{5/2}}\begin{pmatrix}-3xz\cos\varphi-3yz\sin\varphi\\3xz\sin\varphi-3yz\cos\varphi\\x^2+y^2-2z^2 \end{pmatrix}\\ &=\frac{1}{(\rho^2+z^2)^{5/2}}\begin{pmatrix}-3\rho z\cos^2\varphi-3\rho z\sin^2\varphi\\ 3\rho z\cos\varphi\sin\varphi-3\rho z\cos\varphi\sin\varphi\\ \rho^2-2z^2 \end{pmatrix}\\ &=\frac{1}{(\rho^2+z^2)^{5/2}}\begin{pmatrix}-3\rho z\\0\\\rho^2-2z^2\end{pmatrix}\,.\tag{5} \end{align} The divergence in cylindrical coordinates is $$ \nabla\cdot B=\frac{\partial_\rho(\rho B_\rho)}{\rho}+\frac{\partial_\varphi B_\varphi}{\rho}+\partial_z B_z\,.\tag{6} $$ Using (5) this is zero according to Wolfram Alpha as it must.