[proofplan]
The argument is the local Moser trick, with one extra adjustment needed to guarantee that the relevant local flow exists for the whole time interval. Near the chosen point $x_0$, we subtract from $\sigma_t$ a closed $1$-form $\alpha_t$ with the same value at $x_0$; this does not change $d\sigma_t$, but it makes the Moser vector field vanish at $x_0$. A compactness and Gronwall argument then gives a neighbourhood $V$ on which the time-dependent flow exists for all $t\in[0,1]$. Finally, Cartan's formula shows that $\psi_t^*\omega_t$ has zero $t$-derivative, hence remains equal to $\omega_0|_V$.
[/proofplan]
[step:Choose local closed forms that agree with $\sigma_t$ at $x_0$]
Let $m=\dim M$. Choose a coordinate chart $(W,\varphi)$ with $x_0\in W$, $\overline{W}\subset U$ after replacing $W$ by a smaller coordinate neighbourhood if necessary, and write
\begin{align*}
\varphi=(y_1,\dots,y_m):W\to \varphi(W)\subset \mathbb R^m.
\end{align*}
For each $i\in\{1,\dots,m\}$, define the smooth function $a_i:[0,1]\to\mathbb R$ by
\begin{align*}
a_i(t)=\sigma_t(x_0)\bigl((d\varphi_{x_0})^{-1}e_i\bigr),
\end{align*}
where $(e_1,\dots,e_m)$ is the standard ordered basis of $\mathbb R^m$.
Define a smooth family of $1$-forms $\alpha_t\in\Omega^1(W)$ by
\begin{align*}
\alpha_t=\sum_{i=1}^m a_i(t)\,dy_i.
\end{align*}
Since the coefficients $a_i(t)$ are constant as functions on $W$, each $\alpha_t$ is closed:
\begin{align*}
d\alpha_t=0.
\end{align*}
Moreover, by construction,
\begin{align*}
(\alpha_t)_{x_0}=(\sigma_t)_{x_0}.
\end{align*}
Set
\begin{align*}
\beta_t=\sigma_t|_W-\alpha_t\in\Omega^1(W).
\end{align*}
Then $\beta_t$ is a smooth family of $1$-forms on $W$ satisfying
\begin{align*}
d\beta_t=d\sigma_t|_W
\end{align*}
and
\begin{align*}
(\beta_t)_{x_0}=0.
\end{align*}
[guided]
The point of this preliminary step is to use the freedom in the equation
\begin{align*}
\frac{d}{dt}\omega_t=d\sigma_t.
\end{align*}
Only $d\sigma_t$ appears in the hypothesis and in the Moser calculation, so replacing $\sigma_t$ by $\sigma_t$ minus a closed $1$-form changes neither $d\sigma_t$ nor the desired derivative computation.
Choose a coordinate chart $(W,\varphi)$ around $x_0$, with
\begin{align*}
\varphi=(y_1,\dots,y_m):W\to \varphi(W)\subset \mathbb R^m,
\end{align*}
where $m=\dim M$. We shrink $W$ so that $\overline{W}\subset U$; this gives room for the local flow to remain inside $U$ later.
For each coordinate direction $e_i\in\mathbb R^m$, define
\begin{align*}
a_i(t)=\sigma_t(x_0)\bigl((d\varphi_{x_0})^{-1}e_i\bigr).
\end{align*}
Thus $a_i(t)$ is the $dy_i$-coordinate of the covector $(\sigma_t)_{x_0}$. Since $(\sigma_t)_{t\in[0,1]}$ is a smooth family of $1$-forms, each function $a_i:[0,1]\to\mathbb R$ is smooth.
Now define
\begin{align*}
\alpha_t=\sum_{i=1}^m a_i(t)\,dy_i.
\end{align*}
For fixed $t$, the coefficients $a_i(t)$ do not depend on the point of $W$. Therefore
\begin{align*}
d\alpha_t=\sum_{i=1}^m d(a_i(t))\wedge dy_i+\sum_{i=1}^m a_i(t)\,d(dy_i)=0.
\end{align*}
The equality $(\alpha_t)_{x_0}=(\sigma_t)_{x_0}$ follows from how the coefficients $a_i(t)$ were chosen in the coordinate coframe at $x_0$.
Finally set
\begin{align*}
\beta_t=\sigma_t|_W-\alpha_t.
\end{align*}
Then
\begin{align*}
d\beta_t=d\sigma_t|_W-d\alpha_t=d\sigma_t|_W
\end{align*}
and
\begin{align*}
(\beta_t)_{x_0}=(\sigma_t)_{x_0}-(\alpha_t)_{x_0}=0.
\end{align*}
This vanishing at $x_0$ is the repair that makes the local-in-space flow exist for the full time interval after shrinking the initial neighbourhood.
[/guided]
[/step]
[step:Define the Moser vector field that vanishes at $x_0$]
For each $t\in[0,1]$, the form $\omega_t|_W$ is symplectic, hence nondegenerate at every point of $W$. Therefore, for each $x\in W$, the [linear map](/page/Linear%20Map)
\begin{align*}
T_xW\to T_x^*W,\qquad v\mapsto \iota_v(\omega_t)_x
\end{align*}
is an isomorphism. Define the time-dependent vector field $X:[0,1]\times W\to TW$ by requiring
\begin{align*}
\iota_{X_t}\omega_t|_W=-\beta_t
\end{align*}
for every $t\in[0,1]$, where $X_t\in\mathfrak X(W)$ denotes the vector field $x\mapsto X(t,x)$.
The inverse of the nondegenerate bundle map depends smoothly on $(t,x)$ because $(\omega_t)$ and $(\beta_t)$ are smooth. Hence $X$ is a smooth time-dependent vector field on $W$. Since $(\beta_t)_{x_0}=0$, the defining equation gives
\begin{align*}
X_t(x_0)=0
\end{align*}
for every $t\in[0,1]$.
[/step]
[step:Shrink the initial neighbourhood so the flow exists for all $t\in[0,1]$]
Choose $r>0$ such that the closed Euclidean ball $\overline{B}(\varphi(x_0),r)$ is contained in $\varphi(W)$. Let
\begin{align*}
Y:[0,1]\times B(\varphi(x_0),r)\to\mathbb R^m
\end{align*}
be the coordinate expression of $X$, defined by
\begin{align*}
Y(t,z)=d\varphi_{\varphi^{-1}(z)}\bigl(X_t(\varphi^{-1}(z))\bigr).
\end{align*}
Since $X_t(x_0)=0$, one has
\begin{align*}
Y(t,\varphi(x_0))=0
\end{align*}
for every $t\in[0,1]$.
By smoothness of $Y$ on the compact set $[0,1]\times\overline{B}(\varphi(x_0),r)$, there exists a constant $L>0$ such that
\begin{align*}
|Y(t,z)-Y(t,z')|\le L|z-z'|
\end{align*}
for all $t\in[0,1]$ and all $z,z'\in \overline{B}(\varphi(x_0),r)$. Choose $\rho>0$ such that
\begin{align*}
0<\rho<r e^{-L}.
\end{align*}
Define
\begin{align*}
V=\varphi^{-1}\bigl(B(\varphi(x_0),\rho)\bigr).
\end{align*}
For $x\in V$, let $z_x:[0,T_x)\to \varphi(W)$ be the maximal solution of the time-dependent [ordinary differential equation](/page/Ordinary%20Differential%20Equation)
\begin{align*}
\frac{d}{dt}z_x(t)=Y(t,z_x(t)),\qquad z_x(0)=\varphi(x).
\end{align*}
Local existence and smooth dependence for time-dependent ordinary differential equations gives such a maximal smooth solution. As long as $z_x(t)\in B(\varphi(x_0),r)$, the Lipschitz estimate and $Y(t,\varphi(x_0))=0$ give
\begin{align*}
\left|\frac{d}{dt}|z_x(t)-\varphi(x_0)|\right|\le L|z_x(t)-\varphi(x_0)|
\end{align*}
at times where $z_x(t)\ne \varphi(x_0)$, and the same estimate follows by the upper Dini derivative at times where equality holds. Gronwall's inequality gives
\begin{align*}
|z_x(t)-\varphi(x_0)|\le e^{Lt}|z_x(0)-\varphi(x_0)|<e^L\rho<r
\end{align*}
for all $t$ for which the solution remains in $B(\varphi(x_0),r)$.
Therefore the solution cannot reach the boundary of $B(\varphi(x_0),r)$ before time $1$. The standard [continuation criterion for ordinary differential equations](/theorems/8314) then implies that $z_x(t)$ exists for every $t\in[0,1]$ and remains in $B(\varphi(x_0),r)$. Define
\begin{align*}
\psi:[0,1]\times V\to U,\qquad \psi(t,x)=\varphi^{-1}(z_x(t)).
\end{align*}
Then $\psi$ is smooth by smooth dependence on initial conditions, and $\psi_0=\operatorname{id}_V$.
[guided]
We now explain why the vector field has a flow for the full interval $[0,1]$ after shrinking the starting neighbourhood. This is the only delicate local point: a general vector field on $U$ may flow a chosen point out of $U$ before time $1$, and shrinking the initial neighbourhood would not fix the trajectory of that point. Here the repair is that $X_t(x_0)=0$ for every $t$.
Choose $r>0$ so that
\begin{align*}
\overline{B}(\varphi(x_0),r)\subset \varphi(W).
\end{align*}
Let
\begin{align*}
Y:[0,1]\times B(\varphi(x_0),r)\to\mathbb R^m
\end{align*}
be the coordinate representative of $X$, namely
\begin{align*}
Y(t,z)=d\varphi_{\varphi^{-1}(z)}\bigl(X_t(\varphi^{-1}(z))\bigr).
\end{align*}
Because $X_t(x_0)=0$, we have
\begin{align*}
Y(t,\varphi(x_0))=0
\end{align*}
for every $t\in[0,1]$.
The map $Y$ is smooth, so its first derivatives in the spatial variables are bounded on the compact set $[0,1]\times\overline{B}(\varphi(x_0),r)$. Hence there is a constant $L>0$ such that
\begin{align*}
|Y(t,z)-Y(t,z')|\le L|z-z'|
\end{align*}
for every $t\in[0,1]$ and all $z,z'\in\overline{B}(\varphi(x_0),r)$. This is the uniform Lipschitz estimate that controls how fast trajectories can separate from the constant trajectory $z(t)=\varphi(x_0)$.
Choose $\rho>0$ with
\begin{align*}
0<\rho<r e^{-L}
\end{align*}
and define
\begin{align*}
V=\varphi^{-1}\bigl(B(\varphi(x_0),\rho)\bigr).
\end{align*}
For $x\in V$, let $z_x:[0,T_x)\to\varphi(W)$ be the maximal solution of
\begin{align*}
\frac{d}{dt}z_x(t)=Y(t,z_x(t)),\qquad z_x(0)=\varphi(x).
\end{align*}
Local existence and smooth dependence for time-dependent ordinary differential equations gives this solution for at least a short time, and it is smooth in the initial point $x$.
As long as the solution lies in $B(\varphi(x_0),r)$, the Lipschitz estimate with $z'=\varphi(x_0)$ gives
\begin{align*}
|Y(t,z_x(t))|\le L|z_x(t)-\varphi(x_0)|.
\end{align*}
Thus the distance from $z_x(t)$ to the fixed point $\varphi(x_0)$ satisfies the differential inequality
\begin{align*}
\left|\frac{d}{dt}|z_x(t)-\varphi(x_0)|\right|\le L|z_x(t)-\varphi(x_0)|
\end{align*}
at all times where $z_x(t)\ne\varphi(x_0)$. At times where the distance is zero, the same estimate is interpreted through the upper Dini derivative, which is enough for Gronwall's inequality. Therefore
\begin{align*}
|z_x(t)-\varphi(x_0)|\le e^{Lt}|z_x(0)-\varphi(x_0)|.
\end{align*}
Since $t\in[0,1]$ and $x\in V$, this yields
\begin{align*}
|z_x(t)-\varphi(x_0)|<e^L\rho<r.
\end{align*}
So the trajectory cannot hit the boundary of the coordinate ball before time $1$.
The ordinary differential equation continuation criterion says that a maximal solution can fail to extend only by leaving every compact subset of the domain of the vector field. Here the trajectory remains in the compact set $\overline{B}(\varphi(x_0),e^L\rho)$, which is contained in $B(\varphi(x_0),r)\subset\varphi(W)$. Hence the solution exists on all of $[0,1]$. Defining
\begin{align*}
\psi(t,x)=\varphi^{-1}(z_x(t))
\end{align*}
gives a smooth map $\psi:[0,1]\times V\to U$, and the initial condition gives $\psi_0=\operatorname{id}_V$.
[/guided]
[/step]
[step:Differentiate the pulled-back forms along the flow]
For each $t\in[0,1]$, write $\psi_t:V\to U$ for $x\mapsto\psi(t,x)$. Since $\psi_t$ is the time-dependent flow of $X_t$, the standard differentiation formula for pullbacks along a time-dependent flow gives
\begin{align*}
\frac{d}{dt}\bigl(\psi_t^*\omega_t\bigr)=\psi_t^*\left(\frac{d}{dt}\omega_t+\mathcal L_{X_t}\omega_t\right).
\end{align*}
The family $(\omega_t)$ is symplectic, so each $\omega_t$ is closed:
\begin{align*}
d\omega_t=0.
\end{align*}
By Cartan's formula for the Lie derivative of differential forms,
\begin{align*}
\mathcal L_{X_t}\omega_t=d(\iota_{X_t}\omega_t)+\iota_{X_t}(d\omega_t).
\end{align*}
Using $d\omega_t=0$ and $\iota_{X_t}\omega_t=-\beta_t$, we obtain
\begin{align*}
\mathcal L_{X_t}\omega_t=-d\beta_t.
\end{align*}
Since $d\beta_t=d\sigma_t|_W$ and $\frac{d}{dt}\omega_t=d\sigma_t$, it follows on $W$ that
\begin{align*}
\frac{d}{dt}\omega_t+\mathcal L_{X_t}\omega_t=d\sigma_t-d\beta_t=0.
\end{align*}
Therefore
\begin{align*}
\frac{d}{dt}\bigl(\psi_t^*\omega_t\bigr)=0
\end{align*}
as a $2$-form on $V$ for every $t\in[0,1]$.
[/step]
[step:Evaluate at $t=0$ to obtain the local symplectic isotopy]
Since the derivative of $t\mapsto \psi_t^*\omega_t$ is zero on $[0,1]$, the family $\psi_t^*\omega_t$ is constant in $t$. Evaluating at $t=0$ and using $\psi_0=\operatorname{id}_V$, we get
\begin{align*}
\psi_t^*\omega_t=\psi_0^*\omega_0=\omega_0|_V
\end{align*}
for every $t\in[0,1]$. Thus $V\subset U$ and $\psi:[0,1]\times V\to U$ satisfy all required conclusions.
[/step]