The integral $2 \iint r^2 \cos\theta\ dxdy $ was correctly evaluated to $4\pi R^4/3$. But it represents only the contribution of the $z$-component of the vector field to the surface integral. By symmetry, the $x$- and $y$- components contribute as much. The total is $4\pi R^4$, as expected.
The issue here is that even if $\mathbf{c}$ is an arbitrary fixed vector, $c_r$ and $c_\theta$ are not constant with respect to $r$ and $\theta$. So, you're missing some other terms in your computation.
Adding the missing terms, we get:
\begin{align}
\nabla \cdot (\mathbf{v}\times\mathbf{c}) &= \frac{\partial (v_{\theta}c_{z}-v_{z}c_{\theta})}{\partial r} + \frac{(v_{\theta}c_{z}-v_{z}c_{\theta})}{r} + \frac{1}{r}\frac{\partial (v_{z}c_{r}-v_{r}c_{z})}{\partial \theta} + \frac{\partial (v_{r}c_{\theta}-v_{\theta}c_{r})}{\partial z}\\
&= c_{r}\left(\frac{1}{r}\frac{\partial v_{z}}{\partial \theta}-\frac{\partial v_{\theta}}{\partial z}\right) + c_{\theta}\left(\frac{\partial v_{r}}{\partial z}-\frac{\partial v_z}{\partial r}-\frac{v_z}{r}\right) + c_z\left(\frac{\partial v_{\theta}}{\partial r} +\frac{v_{\theta}}{r}-\frac{1}{r}\frac{\partial v_{r}}{\partial\theta}\right) + v_z \left(\frac{1}{r}\frac{\partial c_r}{\partial \theta} - \frac{\partial c_\theta}{\partial r} \right).
\end{align}
Now, $\displaystyle \frac{\partial c_\theta}{\partial r} = 0$ and $\displaystyle \frac{\partial c_r}{\partial \theta} = c_\theta$. Thus, this cancels with the extra term that you found, and proves your formula.
Why is $\displaystyle \frac{\partial c_\theta}{\partial r} = 0$?
For each point $P$ in $\mathbb{R}^3$, we have an orthonormal basis $\{ \hat{r}(P), \hat{\theta}(P), \hat{z}(P) \}$ at that point. Now, consider two points $P$ and $Q$ in $\mathbb{R}^3$ such that $Q - P = (\Delta r) \hat{r}(P)$ for some small $\Delta r > 0$. That is, $Q$ is obtained from $P$ by a small shift along the $\hat{r}(P)$ direction. Notice that this does not change the directions of the corresponding unit vectors at $Q$! Hence, the projection of $\mathbf{c}$ along each of the basis vectors $\hat{r}(P), \hat{\theta}(P), \hat{z}(P)$ at $P$ is identical to the projection along the corresponding basis vectors $\hat{r}(Q), \hat{\theta}(Q), \hat{z}(Q)$ at $Q$. Thus, taking the limit as $\Delta r$ tends to zero, we have $\displaystyle \frac{\partial c_\theta}{\partial r} = 0$.
Why is $\displaystyle \frac{\partial c_r}{\partial \theta} = c_\theta$?
Consider two points $P$ and $Q$ in $\mathbb{R}^3$ such that $\angle(\hat{r}(Q),\hat{r}(P)) = \Delta \theta$ for some small $\Delta \theta > 0$. That is, $Q$ is obtained from $P$ by a small rotation in the anti-clockwise direction about the origin. Since the projection of $\mathbf{c}$ along the $\hat{z}$ direction never changes, there is no loss of generality in assuming that $\mathbf{c}$ is an arbitrary fixed vector lying in the plane. Now, a little bit of trigonometry shows that if $\mathbf{c} = c_r \hat{r}(P) + c_\theta \hat{\theta}(P) = \tilde{c_r} \hat{r}(Q) + \tilde{c_\theta} \hat{\theta}(Q)$, then
$$
\frac{\tilde{c_r} - c_r}{\Delta \theta} = c_r \left( \frac{\cos(\Delta \theta) - 1}{\Delta \theta} \right) + c_\theta \frac{\sin (\Delta \theta)}{\Delta \theta}.
$$
Thus, in the limit as $\Delta \theta$ tends to zero, we have $\displaystyle \frac{\partial c_r}{\partial \theta} = c_\theta$.
Best Answer
If one integrates the line integral of the normal component of $\frac{\hat \rho}{\rho}$ along the circle $\rho=\epsilon$, then one finds that
$$\oint_{\rho=\epsilon}\left(\frac{\hat \rho}{\rho}\right)\cdot d\vec\ell=\int_0^{2\pi}\left(\frac{\hat \rho}{\epsilon}\right)\cdot \hat \rho\,\epsilon\,d\phi=2\pi$$
Heuristically, we apply naively the divergence theorem (the conditions of the divergence theorem are not satisfied due to the singularity at $\rho=0$) and find for a smooth function $\phi(\vec \rho)$
$$\begin{align} \lim_{\epsilon\to 0}\int_0^{2\pi}\int_0^\epsilon \nabla\cdot\left(\frac{\hat\rho}{\rho}\right)\phi(\vec\rho)\,dS&\approx \phi(0) \lim_{\epsilon\to 0}\int_0^{2\pi}\int_0^\epsilon \nabla\cdot\left(\frac{\hat\rho}{\rho}\right)\,dS\\\\ &=\phi(0)\lim_{\epsilon\to0}\oint_{\rho=\epsilon}\left(\frac{\hat \rho}{\rho}\right)\cdot d\vec\ell\\\\ &=2\pi\phi(0) \end{align}$$
Hence, in distribution, we can write
$$\nabla\cdot\left(\frac{\hat \rho}{\rho}\right)\sim 2\pi \delta(\vec\rho)$$
To be precise, we begin with the regularization
$$\begin{align} \delta_a(\vec \rho)&=\nabla \cdot \left(\frac{\vec \rho}{\rho^2+a^2}\right)\\\\ &=\frac{2a^2}{(\rho^2+a^2)^2} \end{align}$$
Then, for any test function $\phi(\vec \rho)$ and domain $S\in \mathbb{R}^2$ we have
$$\lim_{a\to0}\int_S \phi(\vec \rho)\delta_a(\vec\rho)\,dS=\begin{cases}2\pi \phi(0)&,\vec\rho\in S\\\\0&,\vec\rho\notin S\end{cases}\tag1$$
To show that $(1)$ holds, we used an analogous approach of the analysis in THIS ANSWER for the $3$-D Dirac Delta.
Thus, it is in this sense that
$$\lim_{a\to0}\delta_a(\vec\rho)=2\pi\delta(\vec\rho)$$