In Appendix A of Density Functional Theory An Advanced Course by Eberhard Engel, Reiner M. Dreizler functional derivative is introduced by studying the Taylor expansion of some functional $F$ around function $f$. $F$ is treated as a function of scalar $\epsilon$, expanded around $\epsilon = 0$:
\begin{equation} F[f + \epsilon\eta] = F[f] + \frac{dF[f + \epsilon\eta]}{d\epsilon} \Bigg \rvert_{\epsilon=0} \epsilon + \frac{1}{2}\frac{d^2F[f + \epsilon \eta]}{d\epsilon^2}\Bigg \rvert_{\epsilon=0} \epsilon^2 + \ldots \end{equation}
In the next step, the following equality is imposed, and if the integral on the right actually exists, the object $\frac{\delta F[f]}{\delta f(x_1)}$ is said to be the (first order) functional derivative of $F$.
\begin{equation} \frac{dF[f + \epsilon\eta]}{d\epsilon}\Bigg \rvert_{\epsilon=0} =: \int dx_1 \frac{\delta F[f]}{\delta f(x_1)} \eta(x_1) \end{equation}
One way to think about functional derivative, the text goes, is to understand it as an expansion of the concept of the total differential of a multivariate function to the case where the number of variables $x_i$ goes to infinity:
\begin{equation} df = \sum_i \frac{\partial{f}}{\partial{x_i}} dx_i \end{equation}
Following this reasoning, the second order functional derivative is defined analogous to how second order total differential is defined.
\begin{equation} \frac{d^2F[f + \epsilon \eta]}{d\epsilon^2}\Bigg \rvert_{\epsilon=0} =: \int dx_1 dx_2 \frac{\delta^2 F[f]}{\delta f(x_1) \delta f(x_2)} \eta(x_1) \eta(x_2) \end{equation}
Although this does intuitively make sense to me, I wonder if there is a more direct way of defining the second order functional derivative, possibly by somehow applying the definition of the first order one? Similarly as the second order total differential can be understood as a differential of the first order total differential:
\begin{equation} d^2f = d (df) = \sum_j \frac{\partial (df)}{\partial x_j} = \sum_{i,j} \frac{\partial^2 f}{\partial x_i \partial x_j} \end{equation}
My attempt:
Assuming that we can find some functional $G$ such that:
\begin{equation} G[f + \epsilon \eta] = \frac{d F[f + \epsilon \eta]}{d \epsilon} \end{equation}
Now, it is not clear why this should be the case in general (although for $\rvert_{\epsilon=0}$, this assumption seems to be fulfilled if the first order functional derivative exists), but assuming it is, we can write:
\begin{equation} \frac{d^2F[f + \epsilon \eta]}{d\epsilon^2}\Bigg \rvert_{\epsilon=0} = \frac{dG[f + \epsilon \eta]}{d\epsilon}\Bigg \rvert_{\epsilon=0} = \int dx_1 \frac{\delta G[f + \epsilon\eta]}{\delta f(x_1)} \eta(x_1)\Bigg \rvert_{\epsilon=0} \\= \int dx_1 \eta(x_1) \frac{\delta}{\delta f(x_1)} \int dx_2 \eta(x_2) \frac{\delta F[f]}{\delta f(x_2)} = \int dx_1 \int dx_2 \frac{\delta^2 F[f]}{\delta f(x_1)\delta f(x_2)} \eta(x_2) \eta(x_1) \end{equation}
Where in the third term, $\epsilon$ in re-introduced in but effectively only evaluated at $\epsilon = 0$. Although the resulting term does correspond to the definition of the second order functional derivative, the reasoning behind the steps does not seem very clear and formal.
Is there is a better way to think about this?