[proofplan] Assume $Y := f^{-1}(q)$ is non-empty and fix $p \in Y$. Working in charts around $p$ and $q$, we rearrange coordinates so that $(\partial f_i / \partial x_j)_{1 \le i, j \le k}$ is an invertible $k \times k$ minor of the Jacobian at $p$, then define an auxiliary map $F$ whose first $k$ components are those of $f$ and whose last $n - k$ components are the remaining coordinates of $M$. The [Inverse Function Theorem](/theorems/???) produces a chart on $M$ in which $Y$ is the zero locus of the first $k$ coordinates. Projection onto the last $n - k$ coordinates yields a submanifold chart, so $Y$ is locally an $(n-k)$-dimensional graph and the inclusion is smooth. Hausdorff and second-countable properties of $Y$ are inherited from $M$. [/proofplan]

[step:Reduce to the case where $Y$ is non-empty and fix an arbitrary point of $Y$] If $Y = \varnothing$ there is nothing to prove, so assume $Y \neq \varnothing$ and fix an arbitrary $p \in Y$. The conclusion that $Y$ is an embedded $(n-k)$-submanifold is local: it suffices to exhibit, for every $p \in Y$, an open neighbourhood $W$ of $p$ in $M$ and a chart $(W, \psi)$ on $M$ adapted to $Y$, meaning that $\psi(W \cap Y) = \psi(W) \cap (\{0\}^k \times \mathbb{R}^{n-k})$. Subspace topology inheritance: since $M$ is Hausdorff and second countable and $Y \subseteq M$ carries the subspace topology, $Y$ is automatically Hausdorff and second countable. [/step]

[step:Pass to local coordinates around $p$ and $q$ that send $p$ and $q$ to the origin] Choose a smooth chart $(U, \varphi)$ on $M$ with $p \in U$ and $\varphi(p) = 0 \in \mathbb{R}^n$, and a smooth chart $(V, \theta)$ on $N$ with $q \in V$ and $\theta(q) = 0 \in \mathbb{R}^k$. By shrinking $U$ we may assume $f(U) \subseteq V$. Set \begin{align*} \tilde f: \varphi(U) \subseteq \mathbb{R}^n &\to \mathbb{R}^k \\ x &\mapsto \theta \circ f \circ \varphi^{-1}(x). \end{align*} Then $\tilde f$ is smooth, $\tilde f(0) = 0$, and the regularity hypothesis that $df_p: T_p M \to T_q N$ is surjective translates, under the chart isomorphisms $d\varphi_p: T_p M \xrightarrow{\sim} \mathbb{R}^n$ and $d\theta_q: T_q N \xrightarrow{\sim} \mathbb{R}^k$, to the statement that $J\tilde f_0 \in \mathbb{R}^{k \times n}$ has rank $k$. Write $x = (x_1, \dots, x_n)$ for the standard coordinates on $\mathbb{R}^n$ and $\tilde f = (\tilde f_1, \dots, \tilde f_k)$ for the components. From now on we suppress tildes and write $f$ for $\tilde f$; we also identify $U$ with $\varphi(U) \subseteq \mathbb{R}^n$ and work on $\mathbb{R}^n$. [/step]

[step:Permute coordinates so the first $k \times k$ block of the Jacobian at $0$ is invertible] The Jacobian matrix $Jf_0 \in \mathbb{R}^{k \times n}$ has rank $k$, so it has a non-zero $k \times k$ minor. Let $\{j_1 < j_2 < \dots < j_k\} \subseteq \{1, \dots, n\}$ be a set of column indices for which the corresponding $k \times k$ submatrix of $Jf_0$ is invertible. Let $\sigma: \{1, \dots, n\} \to \{1, \dots, n\}$ be any permutation with $\sigma(i) = j_i$ for $1 \le i \le k$, and consider the linear diffeomorphism \begin{align*} P_\sigma: \mathbb{R}^n &\to \mathbb{R}^n \\ (x_1, \dots, x_n) &\mapsto (x_{\sigma(1)}, \dots, x_{\sigma(n)}). \end{align*} Replacing $\varphi$ by $P_\sigma \circ \varphi$ (a smooth chart because $P_\sigma$ is a linear isomorphism) we may assume, without loss of generality, that the $k \times k$ minor \begin{align*} \left( \frac{\partial f_i}{\partial x_j}(0) \right)_{1 \le i, j \le k} \end{align*} is invertible. [/step]

[step:Extend $f$ to a map $F: \mathbb{R}^n \to \mathbb{R}^n$ whose differential at $0$ is an isomorphism] Define \begin{align*} F: U \subseteq \mathbb{R}^n &\to \mathbb{R}^n \\ x = (x_1, \dots, x_n) &\mapsto (f_1(x), \dots, f_k(x), x_{k+1}, \dots, x_n). \end{align*} Then $F$ is smooth on $U$ and $F(0) = 0$. The Jacobian matrix of $F$ at $0$ is the block matrix \begin{align*} JF_0 = \begin{pmatrix} A & B \\ 0 & I_{n-k} \end{pmatrix}, \end{align*} where $A = \left( \partial f_i / \partial x_j (0) \right)_{1 \le i, j \le k} \in \mathbb{R}^{k \times k}$, $B = \left( \partial f_i / \partial x_j (0) \right)_{1 \le i \le k,\, k+1 \le j \le n} \in \mathbb{R}^{k \times (n-k)}$, and $I_{n-k}$ is the identity matrix of size $n-k$. By the block triangular structure, \begin{align*} \det(JF_0) = \det(A) \cdot \det(I_{n-k}) = \det(A) \neq 0, \end{align*} so $JF_0$ is invertible and consequently $DF_0: \mathbb{R}^n \to \mathbb{R}^n$ is a linear isomorphism. [/step]

[step:Apply the Inverse Function Theorem to obtain a local diffeomorphism] We apply the [Inverse Function Theorem](/theorems/???) to $F$ at $0$. The hypotheses are: $F$ is smooth on an open neighbourhood of $0$ (satisfied, since $F$ is a smooth map $U \to \mathbb{R}^n$ with $0 \in U$), and $DF_0$ is a linear isomorphism (verified in the previous step). The theorem provides an open neighbourhood $W \subseteq U$ of $0$ such that $F(W)$ is open in $\mathbb{R}^n$ and $F|_W: W \to F(W)$ is a smooth diffeomorphism. Consequently the components $(f_1, \dots, f_k, x_{k+1}, \dots, x_n)$ of $F$ form a smooth chart on $M$ centred at $p$, with chart domain $W$ (identified back with the corresponding open subset of $M$ via $\varphi$) and coordinate map $F: W \to F(W)$. [/step]

[step:Identify $Y \cap W$ as the zero locus of the first $k$ new coordinates] In the new coordinates $y = F(x) = (y_1, \dots, y_n)$, we have $y_i = f_i(x)$ for $1 \le i \le k$ and $y_i = x_i$ for $k+1 \le i \le n$. Therefore, for $x \in W$, \begin{align*} x \in Y \cap W \iff f(x) = 0 \iff f_1(x) = \dots = f_k(x) = 0 \iff y_1 = \dots = y_k = 0. \end{align*} Hence \begin{align*} F(Y \cap W) = F(W) \cap \left( \{0\}^k \times \mathbb{R}^{n-k} \right). \end{align*} This is precisely the condition that $(W, F)$ is a submanifold chart for $Y$ adapted to the embedding $Y \hookrightarrow M$. [/step]

[step:Conclude that $Y$ is an $(n-k)$-dimensional embedded submanifold with smooth inclusion] We have exhibited, for every $p \in Y$, a submanifold chart $(W, F)$ such that $F(Y \cap W) = F(W) \cap (\{0\}^k \times \mathbb{R}^{n-k})$. Composing $F|_{Y \cap W}$ with the projection \begin{align*} \pi: \{0\}^k \times \mathbb{R}^{n-k} &\to \mathbb{R}^{n-k} \\ (0, \dots, 0, y_{k+1}, \dots, y_n) &\mapsto (y_{k+1}, \dots, y_n) \end{align*} yields a chart $(W \cap Y, \pi \circ F|_{Y \cap W})$ on $Y$ taking values in an open subset of $\mathbb{R}^{n-k}$. Two such charts $(W_\alpha \cap Y, \pi \circ F_\alpha)$ and $(W_\beta \cap Y, \pi \circ F_\beta)$ have smooth transition: their composition factors through $F_\beta \circ F_\alpha^{-1}$, which is smooth because $F_\alpha$ and $F_\beta$ are smooth charts on $M$. Therefore these charts form a smooth atlas making $Y$ an $(n-k)$-dimensional smooth manifold. Finally, in the chart $(W, F)$ on $M$ and the chart $(W \cap Y, \pi \circ F|_{Y \cap W})$ on $Y$, the inclusion $\iota: Y \hookrightarrow M$ is represented by \begin{align*} (y_{k+1}, \dots, y_n) \mapsto (0, \dots, 0, y_{k+1}, \dots, y_n), \end{align*} which is a smooth linear map. Thus $\iota$ is smooth, completing the proof that $Y = f^{-1}(q)$ is an $(n-k)$-dimensional embedded submanifold of $M$. [/step]