Mixed second functional derivative not symmetrical with respect to order of differentiation

Question

I have a functional $E$, which is a functional of two different functions and their gradients: $$ E[\psi_1,\psi_2] = \int d^3\mathbf{x} ~ \varepsilon\{\mathbf{x}\}$$ where I'm using $\varepsilon\{\mathbf{x}\}$ as short-hand for $$ \varepsilon\{\mathbf{x}\} = \varepsilon\left(\mathbf{x};\psi_1(\mathbf{x}),\nabla\psi_1(\mathbf{x}),\psi_2(\mathbf{x}),\nabla\psi_2(\mathbf{x})\right) .$$ (This is coming from a density functional theory problem where the $\psi$ are orbitals, so there are actually many of them, but just two illustrates the problem.)

I would like to calculate the mixed second functional derivative of $E$ $$ \frac{\delta^2 E}{\delta \psi_1 (\mathbf{r}) \delta \psi_2 (\mathbf{r}')}. $$ As I understand it, the two functional derivatives here should commute, so the order of taking these derivatives shouldn't matter - however, by following what seems like reasonably simple steps, I find that this is not the case (see below for my work so far). I cannot find anything on mixed functional derivatives like this after a lot of searching, but potentially I have missed something. Can anyone explain what I am doing wrong?

Derivation:

First, I take the functional derivative of $E$ wrt $\psi_1(\mathbf{r})$. The result of this kind of operation is fairly easy to find (e.g. http://mbpt-domiprod.wikidot.com/calculation-of-gga-kernel), and we get: $$ \frac{\delta E}{\delta \psi_1(\mathbf{r})} = \frac{\partial \varepsilon\{\mathbf{r}\}}{\partial \psi_1 (\mathbf{r})} - \nabla.\frac{\partial\varepsilon\{\mathbf{r}\}}{\partial\nabla\psi_1(\mathbf{r})} = f(\mathbf{r}). $$ Now, we need to take the functional derivative of $f(\mathbf{r})$ wrt $\psi_2(\mathbf{r}')$. We write $f(\mathbf{r})$ as a functional: $$ f(\mathbf{r}) \rightarrow F_\mathbf{r}[f] = \int d^3\mathbf{x} ~ f(\mathbf{x}) \delta(\mathbf{r}-\mathbf{x}). $$ Using the same result as before, with $\varepsilon\{\mathbf{x}\}$ replaced by $f(\mathbf{x})\delta(\mathbf{r}-\mathbf{x})$, we can write the second derivative as \begin{align} \frac{\delta^2 E}{\delta \psi_1 (\mathbf{r}) \delta \psi_2 (\mathbf{r}')} &= \frac{\partial}{\partial \psi_2 (\mathbf{r}')} \left[f(\mathbf{r}')\delta(\mathbf{r}-\mathbf{r}')\right] - \nabla'.\frac{\partial}{\partial\nabla'\psi_2(\mathbf{r}')} \left[f(\mathbf{r}')\delta(\mathbf{r}-\mathbf{r}')\right] \\ &= \frac{\partial}{\partial \psi_2 (\mathbf{r}')} \left[\frac{\partial \varepsilon\{\mathbf{r'}\}}{\partial \psi_1 (\mathbf{r'})}\delta(\mathbf{r}-\mathbf{r}') - \nabla'.\frac{\partial\varepsilon\{\mathbf{r}'\}}{\partial\nabla'\psi_1(\mathbf{r}')}\delta(\mathbf{r}-\mathbf{r}')\right] \\ & ~ ~ ~ ~ - \nabla'.\frac{\partial}{\partial\nabla'\psi_2(\mathbf{r}')} \left[\frac{\partial \varepsilon\{\mathbf{r'}\}}{\partial \psi_1 (\mathbf{r'})}\delta(\mathbf{r}-\mathbf{r}') - \nabla'.\frac{\partial\varepsilon\{\mathbf{r}'\}}{\partial\nabla'\psi_1(\mathbf{r}')}\delta(\mathbf{r}-\mathbf{r}')\right] \\ &= \frac{\partial^2 \varepsilon\{\mathbf{r'}\}}{\partial \psi_1 (\mathbf{r'}) \partial \psi_2 (\mathbf{r'})}\delta(\mathbf{r}-\mathbf{r}') - \nabla'.\frac{\partial^2 \varepsilon\{\mathbf{r}'\}}{\partial\nabla'\psi_1(\mathbf{r}')\partial\psi_2(\mathbf{r}')}\delta(\mathbf{r}-\mathbf{r}') \\ & ~ ~ ~ ~ - \nabla'.\left[ \frac{\partial^2 \varepsilon\{\mathbf{r'}\}}{\partial \psi_1 (\mathbf{r'})\partial \nabla'\psi_2(\mathbf{r}')}\delta(\mathbf{r}-\mathbf{r}') \right] + \nabla_2'.\left[ \nabla_1'.\frac{\partial^2 \varepsilon\{\mathbf{r}'\}}{\partial\nabla_1'\psi_1(\mathbf{r}')\partial\nabla_2'\psi_2(\mathbf{r}')}\delta(\mathbf{r}-\mathbf{r}')\right]. \end{align} I've used subscripts to make it clear which $\nabla$ is dotted with which in the 4th term. For the final step, I've allowed $\nabla$ and the derivative wrt $\psi_2$ or $\nabla\psi_2$ to commute in the 2nd and 4th terms, which I believe is valid, but this is the step I'm least sure on. However, I don't think this affects the main issue here, which is to do with the $\delta$-functions.

This expression is clearly not symmetrical in the order of derivation, because of the $\delta$-function inside the derivatives - if I did the whole thing again but differentiating wrt $\psi_2$ first, the $\delta$-function would move inside the square brackets in the 2nd term, and move outside them in the 3rd term, whilst $\nabla_1$ and $\nabla_2$ would be flipped in the 4th term.

Answer 1

The problem is that OP incorrectly applies the Euler-Lagrange formula for the functional derivative of a local functional to a non-local functional. Instead it is better to use $$ \frac{\delta \phi^{\alpha}(x)}{\delta \phi^{\beta}(y)}~=~\delta^{\alpha}_{\beta}~\delta^n(x-y)$$ and appropriate generalizations of elementary rules in calculus, such as, e.g., the Leibniz rule, cf. e.g. my Phys.SE answer here. Then it is straightforward to check that the 2nd functional derivative is symmetric for a sufficiently smooth functional.

Example:

If $$S[\phi,\psi]~=~\int\!d^nx~{\cal L}(\phi(x),\psi(x),\partial\phi(x),\partial\psi(x),x)$$ is a local functional of first order, then the 1st functional derivative is the Euler-Lagrange (EL) formula $$ \frac{\delta S}{\delta\phi(x)}~=~\frac{\partial{\cal L}(x)}{\partial\phi(x)}-\frac{d}{dx^{\mu}}\frac{\partial{\cal L}(x)}{\partial\phi_{,\mu}(x)}, $$ and then the 2nd functional derivative becomes symmetric: $$\begin{align}&\frac{\delta^2 S}{\delta\psi(y)\delta\phi(x)}\cr ~=~&\frac{\partial^2{\cal L}(x)}{\partial\psi(x)\partial\phi(x)}\delta^n(x\!-\!y)+\frac{\partial^2{\cal L}(x)}{\partial\psi_{,\nu}(x)\partial\phi(x)}\frac{d}{dx^{\nu}}\delta^n(x\!-\!y)\cr ~-~&\frac{d}{dx^{\mu}}\left[\frac{\partial^2{\cal L}(x)}{\partial\psi(x)\partial\phi_{,\mu}(x)}\delta^n(x\!-\!y)+\frac{\partial^2{\cal L}(x)}{\partial\psi_{,\nu}(x)\partial\phi_{,\mu}(x)}\frac{d}{dx^{\nu}}\delta^n(x\!-\!y)\right]\cr ~=~&\frac{\partial^2{\cal L}(x)}{\partial\psi(x)\partial\phi(x)}\delta^n(x\!-\!y)\cr ~-~&\frac{d}{dy^{\nu}}\left[\frac{\partial^2{\cal L}(x)}{\partial\psi_{,\nu}(x)\partial\phi(x)}\delta^n(x\!-\!y)\right]-\frac{d}{dx^{\mu}}\left[\frac{\partial^2{\cal L}(x)}{\partial\psi(x)\partial\phi_{,\mu}(x)}\delta^n(x\!-\!y)\right]\cr ~+~&\frac{d}{dx^{\mu}}\frac{d}{dy^{\nu}}\left[\frac{\partial^2{\cal L}(x)}{\partial\psi_{,\nu}(x)\partial\phi_{,\mu}(x)}\delta^n(x\!-\!y)\right]\cr ~=~& (\phi(x)\leftrightarrow\psi(y)). \end{align}$$ In the last equality we used the identity $$ f(x)\delta^n(x\!-\!y)~=~f(y)\delta^n(x\!-\!y)$$ for the Dirac delta distribution, and that partial derivative commutes for a sufficiently smooth Lagrangian density ${\cal L}$.