I am experiencing a slight dilemma here. Starting from the Divergence theorem, I am trying to show how a divergence-free field $\vec F$ implies that $\vec F$ can be written as the curl of another vector, say $\vec G$. Below I have written some steps which outline my current chain of logic.
Note that I'm assuming we don't already know that $\vec F =\vec\nabla\times \vec G$, since this is what I am trying to prove. Also, I am not using anything like the Helmholtz theorem either.
Here is goes:
Step $0$: Start with the divergence theorem in $\mathbb R^3$:
$$\iint_S\vec F\cdot d\vec S=\iiint _V\text{div}(\vec F)\;dV,$$
where $S$ is of course a closed surface.
Step $1$: If $\vec F$ is a $C^1$ vector field on an open region containing the volume V and $\text{div}(\vec F)=0$ everywhere, then
$$\iint_S\vec F\cdot d\vec S=0,$$
which says the flux through any closed surface must be zero.
Step $2$: Suppose we break the closed surface $S$ into two orientable surfaces $S_1$ and $S_2$ (note that $S_1$ and $S_2$ share the same boundary curve). Let one surface assume the orientation of the other. E.g, if $S_2$ assumes the orientation of $S_1$, then $S_2$ starts off with "negative" orientation, which we will denote by $-S_2$. Since $\iint_{-S_2}\vec F\cdot d\vec S =-\iint_{S_2}\vec F\cdot d\vec S$, hence
$$0=\iint_S\vec F\cdot d\vec S=\iint_{S_1} \vec F\cdot d\vec S + \iint_{-S_2} \vec F\cdot d\vec S=\iint_{S_1} \vec F\cdot d\vec S - \iint_{S_2} \vec F\cdot d\vec S.$$
Therefore,
$$\iint_{S_1}\vec F \cdot d\vec S=\iint_{S_2}\vec F \cdot d\vec S.$$
This says that the flux of a divergence-free vector field through any two open surfaces is surface independent and only depends on the surface boundary. In other words, for any fixed non-empty boundary, we may deform the surface $S$.
You may view this thread for more details regarding this step: Let $F$ be a vector field in $\mathbb{R}^3$. If $F$ is divergence free, we may deform the surface. Why?
Step $3$: This is where I am stuck. I've seen several articles somehow relate this surface independence property back to Stokes' theorem. In particular, I've seen: If $\text{div}(\vec F)=0$ everywhere, then
$$\tag{$\star$}\iint_{\text{any open surface}}\!\!\!\!\vec F\cdot d\vec S=\text{constant}=\oint_{\partial S}\vec G\cdot d\vec r=\iint_S\text{curl}(\vec G)\cdot d\vec S,$$ which implies $\vec F=\vec\nabla\times\vec G$. Here's the link to the article where that's from (on page 3): https://physics56.files.wordpress.com/2016/01/tutorial-7-scalar-and-vector-potential3.pdf
My confusion lies with equation $(\star)$. I just don't see how one simply relates the LHS to Stokes' theorem. I know this is analogous to the case when $\text{curl}(\vec F)=0$, which by Stokes' theorem implies path independence and moreover, the existence of a scalar function $\psi$, such that $\vec F=\vec\nabla \psi$. We can show that this is indeed the case for curl-free fields by using the fundament theorem of calculus. Again, I'm trying to come to this conclusion using vector calculus but not Helmholtz theorem or calculus on manifolds. Just following my outlined chain of logic starting from the Divergence theorem.
I appreciate any input I can get. :)