divergence in polar coordinates
$\DeclareMathOperator\div{div}$The formula for $\nabla\cdot X$ is incorrect. The notation with the 'usual' dot product is misleading. Properly it is: $$\div F = \frac 1\rho\frac{\partial(\rho F^i)}{\partial x^i}$$ where $\rho=\sqrt{\det g}$ is the coefficient of the differential volume element $dV=\rho\, dx^1\wedge\ldots \wedge dx^n$, meaning $\rho$ is also the Jacobian determinant, and where $F^i$ are the components of $F$ with respect to an unnormalized basis.
In polar coordinates we have $\rho=\sqrt{\det g}=r$, and: $$\div X = \frac 1r \frac{\partial(r X^r)}{\partial r} + \frac 1r\frac{\partial(r X^\theta)}{\partial \theta}$$
In the usual normalized coordinates $X=\hat X^{r}\frac{\partial}{\partial r} + \hat X^{\theta}\frac 1r\frac{\partial}{\partial\theta}$ this becomes: $$\div X = \frac 1r \frac{\partial(r \hat X^{r})}{\partial r} + \frac 1r\frac{\partial \hat X^{\theta}}{\partial \theta}$$ which agrees with the usual formula given in calculus books.
If you look at equation (114) here https://mathworld.wolfram.com/CylindricalCoordinates.html
you can see that you said is essentially right but you also need to add on the "connection" terms with the Christoffel symbols. Then you simplify and you get the usual formula.