In this question for example, Divergence-free vector field on a 2-sphere., divergence free vector field on a 2-sphere is shown to be infinite-dimenaional as each $X_h$ is associated with a smooth scalar function $h$.
However, using the direct formula of divergence I got contradiction. Let me specialize to axisymmetric case, i.e. the divergence-free vector field $X$ and everything is only function of $\theta$. $$div(X) = \frac{1}{\sqrt{\sigma}}\partial_\theta (\sqrt{\sigma}X^\theta)=\frac{1}{\sqrt{\sigma}}\partial_\theta (PQ\sin\theta X^\theta)= $$ where I used that axisymmetric metric should be of form $\sigma=P(\theta)^2d\theta^2 + Q(\theta)^2 \sin^2\theta d\phi^2$. Note $\sin\theta=0$ at $\theta=0,\pi$ and $P,Q$ should not vanish anywhere; thus $PQ\sin\theta X^\theta\equiv 0$ and hence $X^\theta\equiv0$! Only trivial divergence vector field is obtained.
Can anyone point our where I could got wrong?
I think the problem arises because you use the polar angle $\theta$ as a dummy index for the summation in your divergence formula. Another problem I see is that in polar coordinates neither a vector field, nor its divergence can be defined at the poles because the azimuth $\varphi$ is not defined there.
From scratch:
In polar coordinates the natural metric on the two-sphere is, as we know, $$ g_{ij}=\begin{pmatrix}1&0\\0&\sin^2\theta\end{pmatrix} $$ where $\theta\in[0,\pi]$ is the polar angle. It sounds like you have in mind a more general axisymmetric metric $$ \sigma_{ij}=\begin{pmatrix}p^2(\theta)&0\\0&q^2(\theta)\sin^2\theta\end{pmatrix}\,. $$ The divergence of a vector field $$ X=\begin{pmatrix}X^\theta\\ X^\varphi\end{pmatrix} $$ is then \begin{align} \operatorname{div}X&=\frac{\partial_\theta(\sqrt{\det\sigma}X^\theta)+\color{red}{\partial_\varphi(\sqrt{\det\sigma}X^\varphi)}}{\sqrt{\det\sigma}}\\[2mm] &=\frac{1}{p(\theta)q(\theta)\sin\theta}\Big\{\partial_\theta\big(p(\theta)q(\theta)\sin\theta X^\theta\big)+\partial_\varphi\big(p(\theta)q(\theta)\sin\theta X^\varphi\big)\Big\}\\[2mm] &=\frac{\partial_\theta(p(\theta)q(\theta)\sin\theta)}{p(\theta)q(\theta)\sin\theta}X^\theta+\partial_\theta X^\theta+\partial_\varphi X^\varphi\\[2mm] &=\partial_\theta\Big(\log (p(\theta)q(\theta)\sin\theta)\Big)X^\theta+\partial_\theta X^\theta+\partial_\varphi X^\varphi\,. \end{align} When $p\equiv q\equiv 1$ (natural metric) this simplifies to $$ \operatorname{div}X=\cot\theta\,X^\theta+\partial_\theta X^\theta+\partial_\varphi X^\varphi\,. $$ An obvious example of a non-zero vector field whose divergence vanishes away from the poles is $$ X=\begin{pmatrix}1\\ -\varphi\cot\theta\end{pmatrix}\,. $$ This $X$ is however not continuous at $\varphi=0$ because $\varphi\in[0,2\pi)$ and the limit from the "left" of $X^\varphi=-\varphi\cot\theta$ at $\varphi=0$ is $-2\pi\cot\theta\,.$
A continuous example is $$ X=\begin{pmatrix}\frac{1}{\sin\theta}\\ c\end{pmatrix} $$ where $c$ is an arbitrary constant. The divergence is zero because $$ \partial_\theta X^\theta=-\frac{\cos\theta}{\sin^2\theta}=-\cot\theta X^\theta\,. $$ This $X$ however does not have removable coordinate singularities at the poles because there it blows up.