After reading the answer here to understand how to apply difference operator, I've figured out how to generalize my proof of Clairaut-Schwarz theorem here to arbitrary order of mixed partial derivatives.
Could you please verify if my proof looks fine or contains logical gaps/errors? Thank you so much for your help!
$\textbf{Generalized Clairaut-Schwarz Theorem:}$ Let $X$ be open in $\mathbb R^n$, $f:X \to F$, and $m \in \mathbb N$. Suppose $j_1, j_2, \ldots, j_m \in\{1,\ldots,n\}$ and $\sigma$ is a permutation of $\{1, \ldots, m\}$. If $\partial_{j_1} \partial_{j_2} \cdots \partial_{j_m} f$ is continuous at $a$ and $\partial_{j_{\sigma(2)}} \cdots \partial_{j_{\sigma(m)}} f$ exists in a neighborhood of $a$, then $$\partial_{j_1} \cdots \partial_{j_m} f (a)= \partial_{j_{\sigma(1)}} \cdots \partial_{j_{\sigma(m)}} f(a)$$
In my proof, I utilize two below lemmas:
Let $\{e_1,\ldots, e_n\}$ be the standard basis of $\mathbb R^n$. For $h \in \mathbb R$ and $j \in \{1,\ldots,n\}$, we define a map $\Delta_j^h f$ by $$\Delta_j^h f: X \to F, \quad x \mapsto f(x+he_j)-f(x)$$
$\textbf{Lemma 1:}$ $$\partial_{j_1} \cdots \partial_{j_m} f (a) = \lim_{h_1 \to 0} \left ( \lim_{h_2 \to 0} \left( \cdots \left ( \lim_{h_m \to 0} \left( \frac{ \Delta_{j_1}^{h_1} \cdots\Delta_{j_m}^{h_m} f (a)}{h_1 \cdots h_m} \right ) \right ) \cdots \right ) \right)$$
$\textbf{Lemma 2:}$ The finite difference operator is commutative, i.e. $$ \Delta_{j_1}^{h_1} \cdots\Delta_{j_m}^{h_m} f (a) = \Delta_{j_{\sigma(1)}}^{h_{\sigma(1)}} \cdots\Delta_{j_{\sigma(m)}}^{h_{\sigma(m)}} f (a)$$
$\textbf{My attempt:}$
By Mean Value Theorem, we have $$\begin{aligned} & \quad \frac{\Delta_{j_{\sigma(1)}}^{h_{\sigma(1)}} \cdots\Delta_{j_{\sigma(m)}}^{h_{\sigma(m)}} f (a)}{h_{\sigma(1)} \cdots h_{\sigma(m)}}\\ =& \quad \frac{\Delta_{j_1}^{h_1} \cdots\Delta_{j_{m}}^{h_{m}} f (a)}{h_1 \cdots h_m} \quad \text{by} \,\, \textbf{Lemma 2} \\ =& \quad \frac{\partial_{j_{1}} \cdots \partial_{j_{m}} f (a + t_1 e_{j_1} + \cdots + t_{m} e_{j_{m}}) h_1 \cdots h_{m}}{h_1 \cdots h_m} \quad \text{by} \,\, \textbf{MVT} \\ =& \quad\partial_{j_{1}} \cdots \partial_{j_{m}} f (a + t_1 e_{j_1} + \cdots + t_{m} e_{j_{m}}) \end{aligned}$$ in which
$$\begin{aligned} \min\{0,h_1\} < t_1 < \max\{0,h_1\} \\ \vdots\quad\quad\quad\quad\quad\quad \,\,\, \\ \min\{0,h_m\} < t_1 < \max\{0,h_m\} \end{aligned}$$
Hence
$$\begin{aligned} & \quad \left \|\frac{\Delta_{j_{\sigma(1)}}^{h_{\sigma(1)}} \cdots\Delta_{j_{\sigma(m)}}^{h_{\sigma(m)}} f (a)}{h_{\sigma(1)} \cdots h_{\sigma(m)}} - \partial_{j_{1}} \cdots \partial_{j_{m}} f (a) \right \|\\ =& \quad \left \| \partial_{j_{1}} \cdots \partial_{j_{m}} f (a + t_1 e_{j_1} + \cdots + t_{m} e_{j_{m}}) - \partial_{j_{1}} \cdots \partial_{j_{m}} f (a) \right \|\end{aligned}$$
Let $t = |t_1| + \cdots+|t_{m}|$ and $h= |h_1| + \cdots+|h_{m}|$. It follows from the continuity of $\partial_{j_1} \partial_{j_2} \cdots \partial_{j_m} f$ at $a$ that for all $\delta > 0$ there is $\epsilon > 0$ such that $$\left \| \partial_{j_{1}} \cdots \partial_{j_{m}} f (a + t_1 e_{j_1} + \cdots + t_{m} e_{j_{m}}) - \partial_{j_{1}} \cdots \partial_{j_{m}} f (a) \right \| <\ \delta$$ for all $t < \epsilon$. As such, for all $h <\ \epsilon$, we have $$ \left \|\frac{\Delta_{j_{\sigma(1)}}^{h_{\sigma(1)}} \cdots\Delta_{j_{\sigma(m)}}^{h_{\sigma(m)}} f (a)}{h_{\sigma(1)} \cdots h_{\sigma(m)}} - \partial_{j_{1}} \cdots \partial_{j_{m}} f (a) \right \| < \delta$$
Take the limit $h_{\sigma(m)} \to 0$, we have $$\lim_{h_{\sigma(m)} \to 0} \left \|\frac{\Delta_{j_{\sigma(1)}}^{h_{\sigma(1)}} \cdots\Delta_{j_{\sigma(m)}}^{h_{\sigma(m)}} f (a)}{h_{\sigma(1)} \cdots h_{\sigma(m)}} - \partial_{j_{1}} \cdots \partial_{j_{m}} f (a) \right \| \le \delta$$
and consequently $$ \left \| \lim_{h_{\sigma(m)} \to 0} \left (\frac{\Delta_{j_{\sigma(1)}}^{h_{\sigma(1)}} \cdots\Delta_{j_{\sigma(m)}}^{h_{\sigma(m)}} f (a)}{h_{\sigma(1)} \cdots h_{\sigma(m)}} \right ) - \partial_{j_{1}} \cdots \partial_{j_{m}} f (a) \right \| \le \delta$$
and consequently $$\left \| \frac{\Delta_{j_{\sigma(1)}}^{h_{\sigma(1)}} \cdots\Delta_{j_{\sigma(m-1)}}^{h_{\sigma(m-1)}} \partial_{j_{m}} f (a)}{h_{\sigma(1)} \cdots h_{\sigma(m-1)}} - \partial_{j_{1}} \cdots \partial_{j_{m}} f (a) \right \| \le \delta \quad \text{by} \,\, \textbf{Lemma 1}$$
Iterating this process of taking limit, we get $$\left \| \lim_{h_{\sigma(1)} \to 0} \frac{\Delta_{j_{\sigma(1)}}^{h_{\sigma(1)}} \left ( \partial_{j_{\sigma(2)}} \cdots \partial_{j_{\sigma(m)}} f \right) (a)}{h_{\sigma(1)}} - \partial_{j_{1}} \cdots \partial_{j_{m}} f (a) \right \| \le \delta$$ or equivalently $$\left \| \lim_{h_{\sigma(1)} \to 0} \frac{ \left ( \partial_{j_{\sigma(2)}} \cdots \partial_{j_{\sigma(m)}} f \right) (a + h_{\sigma(1)} e_{\sigma(1)}) - \left (\partial_{j_{\sigma(2)}} \cdots \partial_{j_{\sigma(m)}} f \right)(a)}{h_{\sigma(1)}} - \partial_{j_{1}} \cdots \partial_{j_{m}} f (a) \right \| \le \delta$$
For all $\delta >0$, there is $\epsilon >0$ such that for all $|h_{\sigma(1)}| < \epsilon$, the last inequality holds. It follows that $$\partial_{j_{\sigma(1)}}\left (\partial_{j_{\sigma(2)}} \cdots \partial_{j_{\sigma(m)}} f \right)(a) = \partial_{j_{1}} \cdots \partial_{j_{m}} f (a)$$ and consequently $$\partial_{j_{\sigma(1)}} \partial_{j_{\sigma(2)}} \cdots \partial_{j_{\sigma(m)}} f (a) = \partial_{j_{1}} \cdots \partial_{j_{m}} f (a)$$
This completes the proof.