Is there any sufficient/necessary condition for a vector field $X$ on some compact Riemannian manifold $M$ (for simplicity just imagine the two-sphere $S^2$) that ensures the existence of a function $h\in C^\infty(M)$ such that
\begin{equation}\nabla h=X,\end{equation} where $\nabla$ is the gradient with respect to the Riemannian metric.
Best Answer
The equation $\nabla h = X$ is equivalent to $dh = X^\flat$, where $X^\flat$ is the 1-form dual to $X$ (this is basically the definition of $\nabla$ on a Riemannian manifold). Thus, your question transforms to the following: for which 1-forms $\omega$ there exists a function $h$ with $dh=\omega$?
By the properties of differentiation of forms, $d\omega=d^2h=0$, hence $d\omega=0$ is a necessary condition for the existence of such $h$. For a form $\omega = \sum \omega_i dx_i$ this condition means
$$\frac{\partial\omega_i}{\partial x_j}=\frac{\partial\omega_j}{\partial x_i}$$
So, there are 1-forms with $\omega = dh$ (called exact forms), which we are interested in, and there are forms with $d\omega=0$ (called closed forms). Any exact form is closed, and we are interested in whether the converse is true. This is precisely what the first de Rham cohomology $H^1_{dR}(M)$ measures: it is defined as the space of closed forms modulo the space of exact forms.
For the case of a 2-sphere $M=S^2$, it is known that $H^1_{dR}(S^2)=0$, meaning that closed and exact forms coincide, so for any closed 1-form $\omega$ there is a function $h$ with $dh=\omega$.
For a 2-torus $M=T^2$ we have $H^1_{dR}(T^2)=\mathbb R^2$, meaning that any closed 1-form decomposes as an exact form plus a linear combination of 2 linearly independent non-exact forms. In fact, using a coordinate system $(\theta, \phi)$ with $\theta, \phi \in [0; 2\pi)$ we can construct a basis of the space of non-exact forms as $\{d\theta, d\phi\}$.