[proofplan]
We prove the backward direction: if the flows $g_1^t$ and $g_2^s$ of two smooth vector fields $V_1$ and $V_2$ commute (i.e., $g_1^t \circ g_2^s = g_2^s \circ g_1^t$ for all $t, s$), then $[V_1, V_2] = 0$. The method is to expand the difference $h(t, s, x_0) = g_1^t(g_2^s(x_0)) - g_2^s(g_1^t(x_0))$ in a Taylor series in $t$ and $s$. The zeroth-order and first-order terms cancel, and the leading cross term proportional to $ts$ is exactly $ts\,[V_1, V_2](x_0)$. Since commutativity forces $h \equiv 0$, this coefficient must vanish at every point, giving $[V_1, V_2] = 0$.
[/proofplan]
[step:Set up the Taylor expansion of the composed flows]
Fix $x_0 \in \mathbb{R}^n$. Since $V_1$ and $V_2$ are smooth, the flow maps $g_1^t$ and $g_2^s$ are smooth in all variables. We expand $g_1^t(g_2^s(x_0))$ as a Taylor series in $t$ and $s$ about $(t, s) = (0, 0)$.
By definition of the flow, $g_i^0 = \operatorname{id}$ and $\frac{\partial}{\partial t}\big|_{t=0} g_1^t(y) = V_1(y)$ for any $y$. Similarly, $\frac{\partial}{\partial s}\big|_{s=0} g_2^s(x_0) = V_2(x_0)$. Therefore, expanding $g_2^s(x_0)$ to first order in $s$:
\begin{align*}
g_2^s(x_0) = x_0 + s\, V_2(x_0) + O(s^2).
\end{align*}
Now expand $g_1^t$ applied to this expression. Writing $y = g_2^s(x_0)$ and expanding $g_1^t(y)$ in $t$:
\begin{align*}
g_1^t(y) = y + t\, V_1(y) + \tfrac{1}{2} t^2 \frac{\partial}{\partial t} V_1(g_1^t(y))\big|_{t=0} + O(t^3).
\end{align*}
We need only the terms through order $ts$. Substituting $y = x_0 + s\, V_2(x_0) + O(s^2)$ and expanding $V_1(y)$ about $x_0$:
\begin{align*}
V_1(y) = V_1(x_0) + s\, (DV_1)_{x_0}(V_2(x_0)) + O(s^2),
\end{align*}
where $(DV_1)_{x_0}$ denotes the Jacobian matrix of $V_1$ evaluated at $x_0$, acting on $V_2(x_0)$. In component form, $((DV_1)_{x_0}(V_2(x_0)))_i = \sum_j \frac{\partial (V_1)_i}{\partial x_j}(x_0) (V_2)_j(x_0)$.
Combining, the expansion of $g_1^t(g_2^s(x_0))$ through order $ts$ is
\begin{align*}
g_1^t(g_2^s(x_0)) = x_0 + t\, V_1(x_0) + s\, V_2(x_0) + ts\, (DV_1)_{x_0}(V_2(x_0)) + O(t^2) + O(s^2).
\end{align*}
[/step]
[step:Expand the flows in the opposite order and extract the $ts$-coefficient of the difference]
By the same argument with the roles of $V_1, t$ and $V_2, s$ exchanged, we obtain
\begin{align*}
g_2^s(g_1^t(x_0)) = x_0 + t\, V_1(x_0) + s\, V_2(x_0) + ts\, (DV_2)_{x_0}(V_1(x_0)) + O(t^2) + O(s^2).
\end{align*}
Subtracting:
\begin{align*}
h(t, s, x_0) &:= g_1^t(g_2^s(x_0)) - g_2^s(g_1^t(x_0)) \\
&= ts \Big[ (DV_1)_{x_0}(V_2(x_0)) - (DV_2)_{x_0}(V_1(x_0)) \Big] + O(t^2 s) + O(t s^2).
\end{align*}
The expression in brackets is the Lie bracket (commutator) of $V_1$ and $V_2$ evaluated at $x_0$. In components:
\begin{align*}
[V_1, V_2]_i(x_0) = \sum_{j=1}^n \left[ (V_2)_j(x_0)\, \frac{\partial (V_1)_i}{\partial x_j}(x_0) - (V_1)_j(x_0)\, \frac{\partial (V_2)_i}{\partial x_j}(x_0) \right].
\end{align*}
(Note: the sign convention here is $[V_1, V_2] = (DV_1)(V_2) - (DV_2)(V_1)$, which is the standard Lie bracket of vector fields.) Therefore
\begin{align*}
h(t, s, x_0) = ts\, [V_1, V_2](x_0) + O(t^2 s) + O(t s^2).
\end{align*}
[/step]
[step:Conclude that commutativity of flows implies vanishing of the Lie bracket]
If the flows commute, then $g_1^t \circ g_2^s = g_2^s \circ g_1^t$ for all $t, s$, so $h(t, s, x_0) = 0$ for all $t, s$ and all $x_0$. In particular, the coefficient of $ts$ in the Taylor expansion of $h$ must vanish. From the previous step, this coefficient is $[V_1, V_2](x_0)$. Since $x_0 \in \mathbb{R}^n$ was arbitrary, we conclude
\begin{align*}
[V_1, V_2] = 0
\end{align*}
identically.
[guided]
More precisely, dividing both sides of $h(t, s, x_0) = 0$ by $ts$ (for $t, s \neq 0$) and taking the limit as $(t, s) \to (0, 0)$ gives
\begin{align*}
[V_1, V_2](x_0) = \lim_{(t,s) \to (0,0)} \frac{h(t, s, x_0)}{ts} = 0.
\end{align*}
This uses the fact that the higher-order remainder terms $O(t^2 s) + O(ts^2)$, when divided by $ts$, tend to zero as $(t, s) \to (0, 0)$, which follows from the smoothness of the flow maps.
[/guided]
[/step]
