Edit: this question is related to other already asked questions, like Symmetry of second functional derivatives , but a clear and definitive answer has never been given, and I am here giving a straightforward, explicit example rather than asking for a general answer/proof as is done in the above link.
I am puzzled because I always thought that functional derivatives were commuting, but I am considering an example where it seems me that it is not the case. Please help me. The action is $$ S\left[g_{\mu\nu},\phi\right]=\int\mathrm{d}^4x\sqrt{-g}F\left(X\right)R\tag{1} $$ for a scalar field $\phi$ and a metric $g_{\mu\nu}$. $R$ is the Ricci scalar while $$X=\nabla_\mu\phi\nabla^\mu\phi.\tag{2}$$ I know that the field equations are of order higher than two but this is not the point. The point is rather that the action depends on $\phi$ only through its derivatives. So both field equations, $$ \frac{\delta S}{\delta g^{\mu\nu}},\quad \frac{\delta S}{\delta \phi}=\sqrt{-g}\Bigl\{-2R\Box\phi F_X-2\nabla_a\phi\nabla^aRF_X-4R\nabla^a\phi\nabla_a\nabla_b\phi\nabla^b\phi F_{XX}\Bigr\}\tag{3} $$ also depend on $\phi$ only through its derivatives. Subscript $X$ is for derivation w.r.t. $X$. I have not written down explicitly the metric field equations since it will not be useful. On the one hand, as regards the second functional derivative $$ \frac{\delta^2 S}{\delta\phi\left(x'\right)\delta g^{\mu\nu}\left(x\right)},\tag{4} $$ it is the functional derivative with respect to $\phi$ of something which depends on $\phi$ only through its derivatives. So it has terms, schematically, like $\nabla\delta\left(x,x'\right)$, $\nabla\nabla\delta\left(x,x'\right)$, but no terms where the Dirac distribution $\delta\left(x,x'\right)$ has no derivatives. On the other hand, $$ \frac{\delta^2 S}{\delta g^{\mu\nu}\left(x\right)\delta\phi\left(x'\right)}\tag{5} $$ contains such pure $\delta\left(x,x'\right)$ terms. Indeed, when the derivative $\frac{\delta}{\delta g^{\mu\nu}}$ hits the term in $F_{XX}$ of $\frac{\delta S}{\delta \phi}$ above, it gives a term including $F_{XXX}$ of the form $$ \sqrt{-g}\Bigl\{-4R\nabla^a\phi\nabla_a\nabla_b\phi\nabla^b\phi\nabla_\mu\phi\nabla_\nu\phi F_{XXX}\delta\left(x,x'\right)\Bigr\},\tag{6} $$ and it is obvious that this is the unique term proportional to $F_{XXX}$, so it cannot be canceled by other terms. This shows that $$ \delta\left(x,x'\right)\subset \frac{\delta^2 S}{\delta g^{\mu\nu}\left(x\right)\delta\phi\left(x'\right)},\quad \delta\left(x,x'\right)\not\subset\frac{\delta^2 S}{\delta\phi\left(x'\right)\delta g^{\mu\nu}\left(x\right)},\tag{7} $$ so the functional derivatives are not commuting. Is this correct?
