Edit: I would like to reopen this question since the linked potential duplicate question is not useful in showing that we get a tempered distribution as the potential, and I have found no easy way to do so.
Edit 2: If $S(\mathbb R^n)$ below were replaced by the space of compactly supported smooth functions $C_{\text c}^\infty(\mathbb R^n)$ everywhere, the question is known to have the answer "yes". See for instance MathOverflow question 230997.
Let $S(\mathbb R^n)$ denote the space of all Schwartz functions on $\mathbb R^n$ equipped with the topology induced by the usual Schwartz semi-norms. Let $S(\mathbb R^n)^*$ denote its dual. For a distribution $T\in S(\mathbb R^n)^*$ and $\phi\in S(\mathbb R^n)$, I will use the notation $$\langle T, \phi\rangle \overset{\text{Def.}}= T(\phi).$$
My question. Suppose that I have a "vector-distribution" $\vec T=(T^1, T^2)\in \big(S(\mathbb R^2)^*\big)^2$ that weakly satisfies $\operatorname{div} T=0$. Explicitly, $$\langle T^1,\partial_1 \phi\rangle+\langle T^2,\partial_2\phi\rangle=0$$ for all $\phi\in S(\mathbb R^2)$. Does it follow that $\vec T$ has a weak potential $\psi\in S(\mathbb R^2)^*$, in the sense that $\nabla^\bot \psi=\vec T$ ? Or, put more explicitly, such that
$$\left\langle\psi,\frac{\partial\phi}{\partial x^2}\right\rangle=\langle T^1,\phi\rangle$$ and $$\left\langle\psi,\frac{\partial\phi}{\partial x^1}\right\rangle=-\langle T^2,\phi\rangle$$ for all $\phi \in S(\mathbb R^2)$ ?
Motivation (the smooth case). Suppose that we are given a smooth vector field $$v=(v^1, v^2)\in\mathcal C^\infty(\mathbb R^2; \mathbb R^2)$$ satisfying $$\operatorname{div} v \overset{\text{Def.}}=\frac{\partial v^1}{\partial x^1}+\frac{\partial v^2}{\partial x^2}=0.$$
We can look at the differential form $$\omega\overset{\text{Def.}}=v^2\,\mathrm dx^1-v^1\,\mathrm dx^2.$$
Formally, $\omega$ is a differential $1$-form over $\mathbb R^2$, i.e. using traditional notation from differential geometry, $$\omega\in\Omega(\mathbb R^2)\overset{\text{Def.}}=\Gamma(\Lambda^2(T^*\mathbb R^2)).$$
We can compute the exterior derivative $$\mathrm d\omega = -\operatorname{div}(v) \,\mathrm dx^1\land\mathrm dx^2=0.$$
Poincaré's Lemma therefore implies that there exists a differential $0$-form, i.e. a function $f\in\mathcal C^\infty(\mathbb R^2)$, such that $$\mathrm df=\omega.$$
In other words,
$$\boxed{-\frac{\partial f}{\partial x^2}=v^1, \frac{\partial f}{\partial x^1}=v^2, \quad\text{i.e.} \quad\nabla^\bot f=v}$$
