[proofplan] The argument has two steps. First, assuming we already have a commuting local frame $\hat E_1, \ldots, \hat E_k$ for $E$ (i.e., $[\hat E_i, \hat E_j] = 0$), the composition of their flows provides a local chart whose first $k$ coordinate vector fields coincide with the $\hat E_i$. Second, we build such a commuting frame from involutivity: starting with a chart in which $E|_p$ aligns with $\partial_{x_1}|_p, \ldots, \partial_{x_k}|_p$, we find a local basis of $E$ of the specific form $\partial_{x_i} + \sum_{l>k} a_{il}\partial_{x_l}$, and the structure of this basis forces $[\hat E_i, \hat E_j]$ to lie outside $E$ unless it is zero. Involutivity then compels it to be zero. [/proofplan]

[step:Commuting flows commute as diffeomorphisms when the fields commute] We record the following ingredient for the construction. Let $V, W$ be smooth vector fields on an open set with $[V, W] = 0$ on a neighbourhood of $p$. Then by the [commuting flows theorem](/theorems/???), there exists $\varepsilon > 0$ such that the flows $\varphi^V_t, \varphi^W_s$ are defined on a common neighbourhood $N$ of $p$ for $|t|, |s| < \varepsilon$ and satisfy \begin{align*} \varphi^V_t \circ \varphi^W_s = \varphi^W_s \circ \varphi^V_t \qquad \text{on } N. \end{align*} [/step]

[step:Build a chart from a commuting frame by composing flows] Assume for now that we have smooth sections $\hat E_1, \ldots, \hat E_k \in \Gamma(E|_U)$ on some open $U \ni p$ that (i) form a pointwise basis of $E_q$ for every $q \in U$, and (ii) satisfy $[\hat E_i, \hat E_j] = 0$ on $U$ for every $1 \le i, j \le k$. We construct the Frobenius chart under this assumption. Extend $\hat E_1(p), \ldots, \hat E_k(p) \in T_pM$ to a basis $\hat E_1(p), \ldots, \hat E_k(p), v_{k+1}, \ldots, v_n$ of $T_pM$ by adjoining vectors $v_{k+1}, \ldots, v_n \in T_pM$. Choose a smooth chart $(\tilde U, \tilde\varphi)$ at $p$ with $\tilde\varphi(p) = 0$ and with coordinate directions $\partial_{\tilde x_l}|_p = v_l$ for $l = k+1, \ldots, n$ (obtained from any initial chart by a linear change of coordinates on $\mathbb{R}^n$). Let $\varphi^{\hat E_i}_{t}$ denote the local flow of $\hat E_i$. Define \begin{align*} F: (-\varepsilon, \varepsilon)^n &\to M \\ (y_1, \ldots, y_n) &\mapsto \bigl(\varphi^{\hat E_1}_{y_1} \circ \cdots \circ \varphi^{\hat E_k}_{y_k}\bigr)\bigl(\tilde\varphi^{-1}(0, \ldots, 0, y_{k+1}, \ldots, y_n)\bigr), \end{align*} for some $\varepsilon > 0$ small enough that all flows are defined and the image stays inside $U \cap \tilde U$. Then $F(0) = p$. [/step]

[step:Compute the partial derivatives of $F$ using commutativity of the flows] For $j \in \{1, \ldots, k\}$, we show $\partial_{y_j} F = \hat E_j \circ F$. Because the flows $\varphi^{\hat E_i}_{y_i}$ pairwise commute by $[\hat E_i, \hat E_j] = 0$ and the previous step, we may reorder the composition to place $\varphi^{\hat E_j}_{y_j}$ on the outside: \begin{align*} F(y) = \bigl(\varphi^{\hat E_j}_{y_j} \circ \varphi^{\hat E_1}_{y_1} \circ \cdots \widehat{\varphi^{\hat E_j}_{y_j}} \cdots \circ \varphi^{\hat E_k}_{y_k}\bigr)\bigl(\tilde\varphi^{-1}(0, \ldots, 0, y_{k+1}, \ldots, y_n)\bigr), \end{align*} where the hat indicates omission. Differentiating in $y_j$ while the other variables are held fixed gives, by the defining property of the flow of $\hat E_j$, \begin{align*} \partial_{y_j} F(y) = \hat E_j\bigl(F(y)\bigr). \end{align*} For $l \in \{k+1, \ldots, n\}$, only the final term $\tilde\varphi^{-1}(0, \ldots, 0, y_{k+1}, \ldots, y_n)$ depends on $y_l$. Evaluating at $y = 0$, \begin{align*} \partial_{y_l} F(0) = (d\tilde\varphi^{-1})_0(e_l) = \partial_{\tilde x_l}\big|_p = v_l. \end{align*} [/step]

[step:Apply the inverse function theorem and pull back coordinates to obtain the Frobenius chart] The differential of $F$ at $0$ has columns $\hat E_1(p), \ldots, \hat E_k(p), v_{k+1}, \ldots, v_n$, a basis of $T_pM$, so $dF_0: \mathbb{R}^n \to T_pM$ is a linear isomorphism. By the [inverse function theorem](/theorems/???), $F$ is a local diffeomorphism onto an open neighbourhood $W$ of $p$. Set $\psi := F^{-1}: W \to (-\varepsilon, \varepsilon)^n$ with coordinates $(x_1, \ldots, x_n) := (y_1, \ldots, y_n) \circ \psi$. From $\partial_{y_j} F = \hat E_j \circ F$ for $1 \le j \le k$, we get $dF_y(\partial_{y_j}|_y) = \hat E_j(F(y))$, so $\partial_{x_j}|_q = \hat E_j(q)$ for every $q \in W$ and $1 \le j \le k$. In particular, \begin{align*} E_q = \operatorname{Span}\langle \hat E_1(q), \ldots, \hat E_k(q)\rangle = \operatorname{Span}\langle \partial_{x_1}|_q, \ldots, \partial_{x_k}|_q \rangle \qquad \text{for every } q \in W. \end{align*} This is the desired chart, provided we exhibit a commuting frame. [/step]

[step:Pick a chart in which $E|_p$ aligns with the first $k$ coordinate directions] It remains to construct a commuting local basis of $E$. Since $E_p \subset T_pM$ is a $k$-dimensional subspace, choose a linear isomorphism $T_pM \cong \mathbb{R}^n$ sending $E_p$ to $\mathbb{R}^k \times \{0\}$. Pulling back a chart through this linear identification, we obtain a smooth chart $(V, \rho)$ at $p$ with coordinates $(x_1, \ldots, x_n)$, $\rho(p) = 0$, such that \begin{align*} E_p = \operatorname{Span}\langle \partial_{x_1}|_p, \ldots, \partial_{x_k}|_p \rangle. \end{align*} [/step]

[step:Construct a local basis of $E$ of adapted form $\hat E_i = \partial_{x_i} + \sum_{l>k} a_{il}\partial_{x_l}$] On $V$, the distribution $E$ is a smooth rank-$k$ subbundle, so it has a smooth local basis $Z_1, \ldots, Z_k \in \Gamma(E|_V)$, each expanded in the coordinate basis as $Z_i = \sum_{m=1}^n Z_{i,m}\, \partial_{x_m}$. At $p$, the matrix $(Z_{i,j}(p))_{1 \le i, j \le k} \in \mathbb{R}^{k \times k}$ has columns giving $Z_i(p)$ projected to $\mathrm{Span}\langle \partial_{x_1}|_p, \ldots, \partial_{x_k}|_p\rangle$. Since $E_p = \mathrm{Span}\langle \partial_{x_j}|_p \rangle_{j \le k}$ and the $Z_i(p)$ form a basis of $E_p$, this matrix is invertible at $p$, hence invertible on a smaller neighbourhood $V' \subset V$ of $p$ by continuity of $\det$. On $V'$, let $B(q) := (Z_{i,j}(q))_{1 \le i, j \le k}^{-1}$ and define \begin{align*} \hat E_i: V' &\to TM|_{V'} \\ q &\mapsto \sum_{j=1}^k B(q)_{ij}\, Z_j(q). \end{align*} Each $\hat E_i$ is a smooth section of $E|_{V'}$, and by construction the first $k$ components of $\hat E_i$ in the coordinate basis are the Kronecker symbols $\delta_{ij}$. Thus \begin{align*} \hat E_i = \partial_{x_i} + \sum_{l = k+1}^n a_{il}(x)\, \partial_{x_l} \end{align*} with smooth coefficients $a_{il} \in C^\infty(V')$. The $\hat E_1, \ldots, \hat E_k$ form a local basis of $E|_{V'}$. [/step]

[step:Show $[\hat E_i, \hat E_j]$ has no $\partial_{x_l}$-component for $l \le k$] Compute the bracket in the coordinate basis. Using the general formula for the Lie bracket of coordinate-expanded vector fields from [theorem 1520](/theorems/1520): if $X = \sum_m X_m \partial_{x_m}$ and $Y = \sum_m Y_m \partial_{x_m}$, then \begin{align*} [X, Y] = \sum_{m=1}^n \left( \sum_{l=1}^n X_l\, \partial_{x_l} Y_m - Y_l\, \partial_{x_l} X_m \right) \partial_{x_m}. \end{align*} Applying this with $X = \hat E_i$, $Y = \hat E_j$, the components $\hat E_{i,m} = \delta_{im}$ for $m \le k$ are constants, hence $\partial_{x_l} \hat E_{i,m} = 0$ for $m \le k$. For $m \le k$, \begin{align*} [\hat E_i, \hat E_j]_m = \sum_{l=1}^n \hat E_{i,l}\, \partial_{x_l} \hat E_{j,m} - \hat E_{j,l}\, \partial_{x_l} \hat E_{i,m} = 0 - 0 = 0. \end{align*} Therefore \begin{align*} [\hat E_i, \hat E_j] = \sum_{l = k+1}^n c_{ij}^{l}(x)\, \partial_{x_l} \end{align*} for some smooth functions $c_{ij}^l \in C^\infty(V')$; in particular, $[\hat E_i, \hat E_j]$ has no component along $\partial_{x_1}, \ldots, \partial_{x_k}$. [/step]

[step:Use involutivity to force $[\hat E_i, \hat E_j] = 0$] By involutivity of $E$ and the fact that $\hat E_1, \ldots, \hat E_k$ form a local basis of $E|_{V'}$, the bracket $[\hat E_i, \hat E_j] \in \Gamma(E|_{V'})$ expands as \begin{align*} [\hat E_i, \hat E_j] = \sum_{s=1}^k \mu_{ij}^s(x)\, \hat E_s = \sum_{s=1}^k \mu_{ij}^s \partial_{x_s} + \sum_{s=1}^k \sum_{l=k+1}^n \mu_{ij}^s\, a_{sl}\, \partial_{x_l} \end{align*} for smooth coefficients $\mu_{ij}^s$. Comparing with the previous step, the components along $\partial_{x_s}$ for $s \le k$ must match: $\mu_{ij}^s = 0$ for every $s \in \{1, \ldots, k\}$. Substituting back, \begin{align*} [\hat E_i, \hat E_j] = 0 \qquad \text{on } V'. \end{align*} The frame $\hat E_1, \ldots, \hat E_k$ on $V'$ satisfies the hypothesis of the chart-building step. [/step]

[step:Assemble the Frobenius chart] Apply the chart-building step (with $U := V'$ and the frame $\hat E_1, \ldots, \hat E_k$ constructed in the last two steps) to obtain an open neighbourhood $W \subset V'$ of $p$ and a smooth chart $\psi: W \to (-\varepsilon, \varepsilon)^n$ with coordinates $(x_1, \ldots, x_n)$ (we reuse the symbol $x$; these are the new Frobenius coordinates) such that \begin{align*} E_q = \operatorname{Span}\langle \partial_{x_1}|_q, \ldots, \partial_{x_k}|_q\rangle \qquad \text{for every } q \in W. \end{align*} The slices $\{(x_1, \ldots, x_n) \in W : x_{k+1} = c_{k+1}, \ldots, x_n = c_n\}$ for each fixed $(c_{k+1}, \ldots, c_n)$ are smooth $k$-dimensional submanifolds of $W$ with tangent space at each point equal to $E_q$, so they are integral manifolds of $E$. Therefore $E$ is integrable, completing the proof. [/step]