Validity of functional derivative using the Dirac delta function In physics, it's customary to compute the functional derivative as $$\frac{\delta F[\rho(x)]}{\delta \rho(y)}=\lim_{\varepsilon\to 0}\frac{F[\rho(x)+\varepsilon\delta(x-y)]-F[\rho(x)]}{\varepsilon}.$$ 
The Wikipedia states that "This works in cases when $F[\rho(x)+\varepsilon f(x)]$ formally can be expanded as a series (or at least up to first order) in $\varepsilon$. ¿Where could I find a proof or a reference to this fact?
 A: In   many variational problems one is given an action functional  $f\mapsto S[f]$, described by an integral
$$ S[f]=\int_\Omega L\bigl(\;x,f(x),D f(x),\dotsc, D^k f(x)\;\bigr) dx $$
in  which 


*

*$\Omega$ is a region in some Euclidean space $\mathbb{R}^n$, $x\in\Omega$, 

*$f$ is a  $k$-times differentiable  function $f: \Omega\to\mathbb{R}^m$,  and 

*the Lagrangian $L$ is a function of the appropriate number of variables.  


For example,  in   classical mechanics $n=1$,  $\Omega\subset \mathbb{R}$  is an interval and the function 
$$f:\Omega\to \mathbb{R}^m, \;\;\Omega\ni t\mapsto f(t)\in\mathbb{R}^m$$
describes a path in  $\mathbb{R}^m$. The Lagrangia has the form $L: \mathbb{R}^m\times\mathbb{R}^m\to \mathbb{R}$, $\newcommand{\bR}{\mathbb{R}}$
$$L(y,v)=\frac{1}{2}|v|^2-U(y), \;\; (y,v)\in\bR^m\times\bR^m $$
where the potential $U$ is a function $U:\bR^m\to\bR$. Then
$$ L(f,\dot{f})= \frac{1}{2}|\dot{f}|^2-U(f), $$ 
where the dot indicates  the derivative with respect to the time parameter $t$ on $\Omega\subset\mathbb{R}$.
The functional (or variational) derivative of $S$ calculated at  $f_0:\Omega\to\mathbb{R}^m$ is  a gadget $\delta S[f_0]$ that  feeds on an infinitesimal deformation $\delta f$ of $f_0$ and returns a scalar
$$ \langle \delta S[f_0], \delta f\rangle =\lim_{h\to 0}\frac{1}{h} \bigl(S[f_0+h\delta f]-S[f_0]\;\bigr). \tag{1}$$
The deformation $\delta f$ is also a function $\Omega\to\mathbb{R}^m$. It is often desirable  to identify  $\delta S[f_0]$ with a function  $g:\Omega\to\mathbb{R}^m$   which, if it exists, is uniquely determined by the equality
$$ \langle \delta S[f_0], \delta f\rangle=\langle g(x), \delta f(x)\rangle =\int_\Omega \bigl( g(x), \delta f(x) \bigr) dx,\tag{2}  $$
where $(-,-)$ denotes the natural inner product on $\mathbb{R}^m$.  The value of $g$ at $x_0$ can be obtained from the  equality
$$ g(x_0)= \langle g(x), \delta(x-x_0)\rangle. \tag{3} $$
This means that the value of $g$ at $x_0$ is obtained by formally replacing $\delta f$ with $\delta(y-x_0)$ in (2).
Making the same formal replacement $\delta f(x)\to\delta(x-x_0)$ in (1) one obtains  the physicists' functional derivative in  your question.
How does one identify $\delta S[f_0]$ with a function?   In the example from classical mechanics  one has
$$ S[f_0](t)=-\frac{d}{dt}\frac{\partial L}{\partial v}(f_0(t), \dot{f_0}(t))+\frac{\partial L}{\partial y}(f_0,\dot{f}_0). $$
How can one see this?  Fix a function (path) $f_0=f_0(t)$. For simplicity set  $\alpha:=\delta f$ and assume $\Omega=[0,1]$. We have a Taylor approximation
$$ L(f_0+h \alpha, \dot{f}_0+h\dot{\alpha}) = L(f_0,\dot{f}_0)+ h\frac{\partial L}{\partial y}(f_0,\dot{f_0})\alpha+h\frac{\partial L}{\partial v}(f_0,\dot{f_0})\dot{\alpha} +O(h^2). $$
Hence
$$\frac{1}{h} \bigl(\; S[f_0+h\alpha]-S[f_0]\;\bigr)=\int_0^1 \frac{\partial L}{\partial y}(f_0,\dot{f_0})\alpha+\frac{\partial L}{\partial v}(f_0,\dot{f_0})\dot \alpha  dt +O(h). $$
Letting $h\to 0$ we deduce
$$\langle \delta S[f_0],\alpha\rangle = \int_0^1\frac{\partial L}{\partial y}(f_0,\dot{f_0})\alpha+\frac{\partial L}{\partial v}(f_0,\dot{f_0})\dot{\alpha}  dt. $$
If we further assume that $\alpha(0)=\alpha(1)$ then upon integrating by parts we deduce
$$ \langle \delta S[f_0],\alpha\rangle=\int_0^1\Bigl(\frac{\partial L}{\partial y}(f_0,\dot{f_0}) -\frac{d}{dt}\frac{\partial L}{\partial v}(f_0,\dot{f_0})\Bigr) \alpha dt. $$
A good place to look for more details is the book "Calculus of Variations" by Gelfand and Fomin, Dover 2000.
