Does this version of Clairaut-Schwarz theorem hold when mixed partial derivatives are of order greater than $2$? I asked this question on MSE here. One person gave an answer but then he deleted it because my version of Clairaut-Schwarz theorem is stronger than his. I meant my version only requires the continuity of one mixed partial derivative while his may requires the continuity of all mixed partial derivatives.
It seems that this question will receive no answer in MSE, so I have no choice but to post it on mathoverflow.net.

I usually encounter Clairaut-Schwarz theorem where the mixed partial derivatives are of order $2$, i.e.

$\textbf{Clairaut-Schwarz Theorem:}$ Let $X$ be open in $\mathbb R^n$, $f:X \to F$, and $i, j \in\{1,\ldots,n\}$. Suppose that $\partial_j \partial_i f$ is continuous at $a$ and that $\partial_j f$ exists in a neighborhood of $a$. Then $\partial_i \partial_j f (a)$ exists and $$\partial_i \partial_j f (a) = \partial_j \partial_i f (a)$$

I would like to ask if Clairaut-Schwarz theorem holds in case the mixed partial derivatives are of arbitrary order $m$, i.e.

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} \partial_{j_2} \cdots \partial_{j_m} f (a)= \partial_{j_{\sigma(1)}} \partial_{j_{\sigma(2)}}  \cdots \partial_{j_{\sigma(m)}} f(a)$$


Thank you so much for your help!
 A: Assume $f:X\to F$ has the partial derivatives $\partial_m\partial_{m-1}\dots\partial_1f(x)$ at any point $x$ of the open rectangle $X:=\prod_{1\le i\le n}]b_i,c_i[$.
For $u\in \mathbb{R}^n$, let $\delta_u $ denote the finite difference $\delta_uf(x):=f(x+u)-f(x)$, defined for $x$ and $x+u\in X$, and let $(e_i)_i$ be the standard basis of $\mathbb{R}^n$.  
Iterating the Mean Value Theorem, we have, say with $|t_j|<r$  and $[x_j-r,x_j+r]\subset\; ]b_j,c_j[$ for $j=1\dots m$
$$\big\|\delta_{t_me_m}\dots \delta_{t_1e_1}f(x)\big\|_F\le |t_1\dots t_m| \sup_{|x_j-y_j| \le r\atop 1\le j\le m}\big\|\partial_m \dots\partial_1f(y) \big\|_F.$$ 
If $p\in F$ and we apply the latter to the function $g(x):=f(x)-px_1\dots x_m $ we also get
$$\big\|\delta_{t_me_m}\dots \delta_{t_1e_1}f(x)-p\,t_1\dots t_m\big\|_F\le |t_1\dots t_m| \sup_{|x_j-y_j| \le r\atop 1\le j\le m}\big\|\partial_m \dots\partial_1f(y) -p\big\|_F.$$
In particular if we take $p:=\partial_m \dots\partial_1f(a)$, and we assume that $\partial_m \dots\partial_1f$ is continuous at $x=a$, we have 
$$\big\|\frac{ \delta_{t_me_m}\dots \delta_{t_1e_1}f(a)} {t_1\dots t_m} -\partial_m \dots\partial_1f(a)\big\|_F=o(1) $$ for $\|t\|\to0$. Assume now there exists the $(m-1)$-order partial derivative  $\partial_{\sigma_{m-1}}\dots\partial_{\sigma_1}f(x)$ for all $x\in X$ and for a given permutation $\sigma$ of $\{1,\dots,m\}$. Since all $\delta_ {t_j e_j}$ commute, taking the limit  for $ t_{\sigma_{j}}\to0$ from $j=1$ to $j=m-1$, we arrive at
$$\big\|\frac{\partial_{\sigma_{m-1}}\dots\partial_{\sigma_1}f(a+t_{\sigma_m})-\partial_{\sigma_{m-1}}\dots\partial_{\sigma_1}f(a)} {t_{\sigma_m}} -\partial_m \dots\partial_1f(a)\big\|_F=o(1) $$ for $t_{\sigma_m}\to0,$ 
whence the thesis: there exists $\partial_{\sigma_{m}}\dots\partial_{\sigma_1}f$ at $x=a$ and coincides with $\partial_m \dots\partial_1f(a)$.
(Note: To simplify the notation I assumed $(j_1,\dots,j_m)$ are  $(1,\dots,m)$; nothing changes in general even if some of the $j_k$ coincide, because any mixed derivative can be written $\partial_{j_m}\dots\partial_{j_1}f(a)=\frac{\partial }{\partial s_m}\dots\frac{\partial }{\partial s_1}f(a+s_1e_{j_1}+\dots+s_m e_{j_m})|_{s_1=0\dots s_m=0}$).  
A: Suppose $\partial_{j_1} \ldots \partial_{j_m} f$ exists in a neighbourhood of $a$ and it is continuous at $a$. Let $\Delta^h_k f(x) = f(x + h e_k) - f(x)$ be the difference operator, with $e_k$ the $k$-th vector in the standard basis. By Fubini's theorem (and induction),
$$ \Delta_{j_1}^{h_1} \ldots \Delta_{j_m}^{h_m} f(a) = \idotsint\limits_{[0,h_1]\times\ldots\times[0,h_m]} \partial_{j_1} \ldots \partial_{j_m} f(a + y) dy_{j_1} \ldots dy_{j_m} ,$$
where we understand that $y_k = 0$ for $k \notin \{j_1, \ldots, j_m\}$. It follows that
$$ \lim_{(h_1, \ldots, h_m) \to 0} \frac{\Delta_{j_1}^{h_1} \ldots \Delta_{j_m}^{h_m} f(a)}{h_1 \ldots h_m} = \partial_{j_1} \ldots \partial_{j_m} f(a) .$$
The difference operators commute, and therefore
$$ \lim_{(h_{\sigma(1)}, \ldots, h_{\sigma(m)}) \to 0} \frac{\Delta_{j_{\sigma(1)}}^{h_{\sigma(1)}} \ldots \Delta_{j_{\sigma(m)}}^{h_{\sigma(m)}} f(a)}{h_{\sigma(1)} \ldots h_{\sigma(m)}} = \partial_{j_1} \ldots \partial_{j_m} f(a) .$$
In particular,
$$ \lim_{h_{\sigma(1)} \to 0} \lim_{(h_{\sigma(2)}, \ldots, h_{\sigma(m)}) \to 0} \frac{\Delta_{j_{\sigma(1)}}^{h_{\sigma(1)}} \ldots \Delta_{j_{\sigma(m)}}^{h_{\sigma(m)}} f(a)}{h_{\sigma(1)} \ldots h_{\sigma(m)}} = \partial_{j_1} \ldots \partial_{j_m} f(a) ,$$
that is,
$$ \lim_{h_{\sigma(1)} \to 0} \frac{1}{h_{\sigma(1)}} \Delta_{j_{\sigma(1)}}^{h_{\sigma(1)}} \biggl(\lim_{(h_{\sigma(2)}, \ldots, h_{\sigma(m)}) \to 0} \frac{\Delta_{j_{\sigma(2)}}^{h_{\sigma(2)}} \ldots \Delta_{j_{\sigma(m)}}^{h_{\sigma(m)}} f(\cdot)}{h_{\sigma(2)} \ldots h_{\sigma(m)}}\biggr)(a) = \partial_{j_1} \ldots \partial_{j_m} f(a) .$$
By assumption, the parenthesized expression is equal to $\partial_{j_{\sigma(2)}} \ldots \partial_{j_{\sigma(m)}} f(\cdot)$ in a neighbourhood of $a$, and so finally
$$ \lim_{h_{\sigma(1)} \to 0} \frac{1}{h_{\sigma(1)}} \Delta_{j_{\sigma(1)}}^{h_{\sigma(1)}} \bigl(\partial_{j_{\sigma(2)}} \ldots \partial_{j_{\sigma(m)}} f\bigr)(a) = \partial_{j_1} \ldots \partial_{j_m} f(a) ,$$
that is,
$$ \partial_{j_{\sigma(1)}} \partial_{j_{\sigma(2)}} \ldots \partial_{j_{\sigma(m)}} f(a) = \partial_{j_1} \ldots \partial_{j_m} f(a) ,$$
as desired.
