(Vishik's Normal Form) Behavior of a vector field near the boundary of a manifold I'm trying to prove a special case of Vishik's Normal Form.
Consider $T^2 := \frac{1}{\sqrt{2}}\cdot\mathbb{T}^2 \subset \mathbb{S}^3$, let $h: \mathbb{\mathbb{S}^3}\to \mathbb{R}$ be the function $h(x_1,x_2,x_3,x_4) = x_1^2+x_2^2 -\frac{1}{2}$, and $\mathbb{S}^+ :=\{x \in \mathbb{S}^3; h(x)\geq 0\}$.
We say that $X$ is a vector field on $\mathbb{S}^+$ if
$$X:\mathbb{S}^+ \to \mathbb{R^4} $$
is a smooth function and $X(p)$ $\in$ $T_p\mathbb{S^3} \subset \mathbb{R}^4$, $\forall p \in \mathbb{S}^+$.
For our purposes we will always consider $X(p) \neq 0$, $\ \forall \ p \in T^2$.

Notations: $Xh(p)= \nabla h(p) \cdot X(p)$ and $X^2 h(p)= \nabla Xh(p) \cdot X(p).$
Definition: Let $X$ be a vector field on $\mathbb{S}^+$ and $p$ $\in$ $T^2$. If $Xh(p) = 0$ and $X^2 h(p) \neq 0$, $p$ is called a fold point.

The theorem that I am trying to demonstrate is as follows:

Theorem (Vishik's Normal Form) Let $X$ be a vector field on $\mathbb{S}^+$ such that, all the points on  $T^2$ that satisfy $Xh(q) =0$ are fold points. Then, if $p$ $\in T^2$ is a fold point, there exists an open set $V_p \subset \mathbb{S}^3$ and a local chart $\varphi: V_p \to U_0 \subset \mathbb{R}^3$ $(\varphi(p) = 0)$  such that $\varphi_*X|_{\varphi(V_p)}$ is a germ at $\{0\}\times \mathbb{R}^2 \cap \varphi (V_p)$  of the vector field given by:
$$\left\{\begin{array}{l}
\dot{x}_1=\varepsilon x_2,\\
\dot{x}_2=1,\\
\dot{x}_3=0. \\
\end{array}\right.$$
Where $\varepsilon =\text{sgn}(X^2 h(p))$ and $\varphi^{-1}(\{0\}\times \mathbb{R}^2 \cap \varphi (V_p) )= V_p\cap T^2$.

I do not have many ideas of how to prove this theorem and I do not know any good reference. Does anyone know how to prove this theorem or can give me some hints (or references)?

My ideas
Note that $h$ is a function such that $0$ is regular value, then using local submersion theorem there exists a local coordinate system $\phi: V_0 \subset \mathbb{R}^3 \to U_p \subset \mathbb{S}^3$ ($\phi(0) = p$), such that $h \circ \phi (x,y,z) = x$.
Defining $Y(q) = D\phi^{-1}(\phi(q)) \cdot X(\phi(q))$, we are able to study the problem in an open neighborhood of $0$ in the topological space $\mathbb{H}^3 = \{(x,y,z); 0\leq x\}$, with some calculation is possible to show that if $f=h\circ \phi$, then $Y f(q) = Y_1(q)$  (where $Y=(Y_1,Y_2,Y_3))$, implying by the definition of $X$, that $Yf(0) =0$  and $Y^2 f(0) = Y(0) \cdot \nabla Y_1(0) \neq 0$.
Once $Yf(0,0,0) =0$  and, by hypotesis, $0\neq Y^2 f(0,0,0) = Y(0,0,0) \cdot \nabla Y_1(0,0,0)$, we are able to conclude that either $$\frac{\partial}{\partial y} Y_1(0,0,0) \neq 0\  \text{or} \ \frac{\partial}{\partial z} Y_1(0,0,0) \neq 0,$$
assuming without loss of generality $\frac{\partial}{\partial y} Y_1(0,0,0)\neq 0$, and using implicit function theorem, there exists $\tau: (-\varepsilon,\varepsilon)^2 \to (-\delta,\delta)$, such that $Y(x,\tau(x,z),z) =0$. Defining the local diffeomorphism $\psi (x,y,z) := (x,y+\tau(x,z),z)$ we can (by changing the coordinates and and shrinking the definition domain of $Y$) assume that $Y = (Y_1,Y_2,Y_3)$, with $Y_1(x,0,z) =0$ (redefining $Y$ as $\psi^{-1}_* Y$).
It is clear that $$\frac{\partial Y_1}{\partial y}(0)\neq 0 \ \text{and}\ \frac{\partial Y_1}{\partial x}(0)  = \frac{\partial Y_1}{\partial x}(0) = 0. $$
Using this question A special change of coordiantes of a Vector Field, we can make a change of coordinates $\zeta (x,y,z) = (x,\pi_{(2,3)} \circ \varphi (y,(x,0,z) )^{-1}$, where $\varphi(t,z)$ is the solution of the ODE
$$
\left\{\begin{array}{l}
\dot{x} = Y(x)\\
x(0)=z,\\
\end{array}\right.
$$
(we are able to extend $Y(x)$ to a open neighboorhood of $(0,0,0)$), implying that
$$\zeta_* Y(x,y,z) = \text{d}\zeta (\zeta^{-1} (x)) Y(\zeta^{-1}(x))=  \left(\begin{array}{c}
Y_1\circ \zeta^{-1}(x,y,z)\\
1\\
0 \\
\end{array}\right)$$
Note that we still have $Y_1(\zeta^{-1}(x,0,z)) = Y_1(x,\pi_{(2,3)}\circ \varphi(0,(x,0,z)) = Y_1(x,0,z) =0,$ and consequently
$$\frac{\partial Y_1(\zeta^{-1})}{\partial x}(0) = \frac{\partial Y_1(\zeta^{-1})}{\partial z}(0) =0 \ \text{and}\ \frac{\partial Y_1(\zeta^{-1})}{\partial y}(0)\neq 0.$$
The vector field is almost in the desired form, but I don't know how to proceed, any hints?

Another thing that I noticed but didn't lead me to anywhere, is that we can define the local diffeomorphism $\theta (x,y,z) = (x, Y_1\circ \zeta^{-1}(x,y,z),z)$ and we can easily check that
\begin{align*}
\theta_* ( \zeta_* Y)(x,y,z) &=  \left(\begin{array}{c}
y\\
\frac{\partial Y_1\circ \zeta^{-1} }{\partial x}\left(\theta^{-1}(x,y,z)\right) \cdot  Y_1\circ \zeta^{-1} \circ \theta^{-1} (x,y,z) + \frac{\partial Y_1\circ \zeta^{-1}}{\partial y}(\theta^{-1}(x,y,z)) \\
0 \\
\end{array}\right)\\
&=  \left(\begin{array}{c}
y\\
\nabla Y_1(\zeta^{-1})(\theta^{-1} (x,y,z)) \cdot \zeta_* Y (\theta^{-1}(x,y,z)) \\
0 \\
\end{array}\right)\\
&=  \left(\begin{array}{c}
y\\
  (\zeta_* Y)^2 \pi_1 \  (\theta^{-1}(x,y,z)) \\
0 \\
\end{array}\right)\\
&=  \left(\begin{array}{c}
y\\
  (X^2 h) \  (\theta \circ \zeta \circ \phi )^{-1}(x,y,z)) \\
0 \\
\end{array}\right)\\
\end{align*}
however $\frac{\partial Y_1\circ \zeta^{-1} }{\partial x}\left(\theta^{-1}(x,y,z)\right) \cdot  Y_1\circ \zeta^{-1} \circ \theta^{-1} (x,y,z) + \frac{\partial Y_1\circ \zeta^{-1}}{\partial y}(\theta^{-1}(x,y,z))$ does not seem like a constant.
 A: I think you are almost there. So, just a hint:
Think on the coordinate system $\tilde\varphi:V_p\to U$ satisfying $d\tilde\varphi(X)=(0,1,0)$. If you replace the first coordinate by the function $h$, you get something whose first jets at $p$ is $(ax_2,1,0)$, $a=X^2h(p)$. It shouldn't be hard to proceed from here.
The geometric meaning of the situation is that $h$ is the coordinate normal to the Torus. So, if you go in the direction of $X$, you start escaping slowly, with starting escaping velocity 0, but accelerating. Note that $Xh(p)=0$ just mean that $X$ is tangent to the torus; $X^2h(p)\neq 0$ means it is trying to escape.
ERRATA v.21.10.2018 (changing the construction from $\{0\}\times \mathbb{R}^2$ to $\{x\}\times \mathbb{R}^2$)
Actually, the diffeomorphism I suggested does not send $T^2\cap V_p$ to $\{0\}{\times} \mathbb R^2$, but the greatest issue is still the first coordinate. Let us think on the other way around:
Note that $V_p$ can be thought as a a neighborhood $U'\subset \mathbb R^3$ (up to diffeomorphism), where $T^2\cap V_p$ corresponds to $U'\cap \{0\}\times\mathbb R^2$ and the first coordinate is $h$. Let us transpose the whole problem to $U'$, still denoting by $X$ and $h$ their corresponding objects in $ U'$. Decompose $X=(X_1,X_2,X_3)$ and denote $X^T=(0,X_2,X_3)$. 
The function $h$ satisfies $h(T^2)=0$, $Xh(p)=0$ and $X^2h(p)\neq 0$ (let us assume $X^2h(p)> 0$, for simplicity). The last inequality guarantees that $p$ is a regular point for the function $h_1=Xh$, therefore,  $C'=h_1^{-1}(0)\cap U'$ is a regular surface transversal to $X$ and  $X^T$ (you can shrink $U'$ for that). In particular, $C_x=C'\cap h^{-1}(x)$ is a smooth family of curves on $U'$ for small $x$. 
Shrinking $U'$ again, we assume both $X^Th_1$ and $X^2h$ positive on $C'$. 
Let $\psi:U'\to U''\subset \mathbb R^3$ be a diffeomorphism where $h$ is the  first coordinate (thus preserving the first coordinate), $d\psi(X^T)=(0,1,0)$ and $z\mapsto\psi(x,0,z)$ parametrizes $C_x$ (it exists since the flow of $X^T$ preserves the first coordinate and $C_x$ is transverse to $X^T$). Now let us see what we need to  change in the first coordinate. Given $(x_0,0,z_0)\in C_x$,  $h_2=h\circ\psi^{-1}=first~coordinate$ satisfies 
$$\psi(X)h_2(x_0,y,z+z_0)=Xh(\psi^{-1}(x_0,y,z+z_0))\\=Xh(\psi^{-1}(x_0,0,z_0))+yX^Th_1(\psi^{-1}(x_0,0,z+z_0))+z\frac{d}{dz}Xh(\psi^{-1}(x_0,0,z+z_0))+o(x_0,y,z+z_0)\\=yX^Th(\psi^{-1}(x_0,0,z+z_0))+o(x_0,y,z+z_0)
$$
where $\lim_{(y,z)\to(0,0)}o(y,z)/|(y,z)|=0$ and $Xh(\psi^{-1}(0,0,z_0))=\frac{d}{dz}h(\psi^{-1}(0,y,z+z_0))=0$ since $Xh$ vanishes identicaly on $C'$.
Now we get rid of $o(x_0,y,z+z_0)$.
We just integrate $o$  on the $\psi(X)$-direction and subtract the result from $h_2$. That would give you the desired  first coordinate, and if you know how to do this integration, I highly recommend to not proceed and do it yourself, since the text gets messy.
To make this integration, consider the chart  $\varphi:U'''\to U''$ where $\varphi(0,0,z)\in C$ and $\varphi(t,y,z)=c(t)$, $ \dot c=\psi(X)$. Take the function $O:U''\to \mathbb R$:
$$O(x,y,z)=\int_0^{\varphi^{-1}_1(x,y,z)} o(\varphi(t,\varphi^{-1}_2(x,y,z),\varphi^{-1}_3(x,y,z)))dt.$$
If nothing went wrong, $\psi(X)(h_2-O)=\psi(X)h_2-o$, as desired.
Completing the proof:
Consider the coordinate chart $\tilde \psi:U'\to V'$ that you get by replacing the first coordinate by $h_3=h_2-O$. We have:
$$\tilde\psi(X)h_3(x_0,y,z+z_0)=Xh(\tilde\psi^{-1}(x_0,y,z+z_0))-\tilde\psi(X)O(\tilde\psi^{-1}(x_0,y,z+z_0))\\=yX^Th_1(\tilde\psi^{-1}(x_0,0,z_0))+o(\tilde\psi^{-1}(x_0,y,z+z_0))-XO(\tilde\psi^{-1}(x_0,y,z+z_0))\\=yX^Th(\tilde\psi^{-1}(x_0,0,z+z_0)).
$$
(I understand that you be not happy with this term multiplying $y$.  You can try to replace $O$ by another function to really finish the proof. I will try something along the week.)
A: See the following paper on  ArXiv.

