In the paper they include also the following:
Since $\nabla_k X_j$ is symmetric in $k$ and $j$ we can integrate $X_i$ locally, i.e. on any simply connected domain $\Omega \subset N$ we can find a smooth function $f$
such that $X_i = \nabla_i f$.
This is indeed how one proves the Poincare lemma for 1 one form. Indeed you are always dealing with the one form $X = X_i dx^i$. Since
$$ (dX)_{ij} = \partial _i X_j - \partial _j X_i = \nabla _j X_i - \nabla _i X_j = 0,$$
$X$ is a closed one form and thus $X = df$ for some smooth function $f$.
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\}$.
Best Answer
Since the Riemannian metric $g$ is non-degenerate, it induces an isomorphism $X\in\mathfrak{X}(M) \mapsto X^\flat\in\Omega(M)$ given by $$ X^\flat(Y) := g(X,Y). $$ The potential condition can be written as $X^\flat = d\psi$. This is true locally if $dX^\flat = 0$. It's true globally if in addition $H^1(M) = \lbrace 0 \rbrace$.