quantum mechanics - Variational proof of the Hellmann-Feynman theorem

I use the following notation and definition for the (first) variation of some functional $E[\psi]$: \begin{equation} \delta E[\psi][\delta\psi] := \lim_{\varepsilon \rightarrow 0} \frac{E[\psi + \varepsilon \delta\psi] - E[\psi]}{\varepsilon}. \end{equation} For a Hamiltonian that depends on some parameter, $H(\lambda)$, which fulfills the eigenequation $H(\lambda) |\psi(\lambda)\rangle = E(\lambda) |\psi(\lambda)\rangle$, one can quite easily show that \begin{equation} \frac{dE(\lambda)}{d\lambda} = \langle \psi(\lambda) | \frac{dH(\lambda)}{d\lambda} | \psi(\lambda) \rangle \end{equation} which is called the Hellmann-Feynman theorem. Now in quantum chemistry, one often seeks an approximation of the true eigenfunctions by making the functional for the energy expectation value stationary, \begin{equation} E[\psi] = \frac{\langle \psi | H | \psi \rangle}{\langle \psi | \psi \rangle} \end{equation} \begin{equation} \delta E[\psi][\delta \psi] = \frac{1}{\langle \psi | \psi \rangle}\left[ \langle \delta\psi | H - E[\psi] | \psi\rangle + \langle \psi | H - E[\psi] | \delta\psi\rangle \right] = 0 \end{equation} for all variations $\delta\psi$ in some subspace (called the variational space) of the full Hilbert space. In this case, one often reads that the Hellmann-Feynman theorem is still valid for the approximate wave functions and energies.

Unfortunately, I have problems finding a clear proof of this statement in books. Wikipedia tries a proof (section "Alternate proof"), but it assumes that the functional derivative \begin{equation} \frac{\delta E[\psi]}{\delta\psi} \end{equation} exists, which is not the case because of the antilinearity in bra vectors, which leads to \begin{equation} \delta E[\psi][\alpha \cdot \delta \psi] \neq \alpha \cdot \delta E[\psi][ \delta \psi] \end{equation} for complex $\alpha$. Can anyone provide me with a clear justification of why the Hellmann-Feynman theorem works for variational wavefunctions?

Answer

Comment to the question (v2): One should not require the functional derivative $\frac{\delta E[\psi]}{\delta\psi}$ to be complex differentiable/holomorhic. That's an impossible requirement for a real functional $$E[\psi]~:=~ \frac{\langle \psi | H | \psi \rangle}{\langle \psi | \psi \rangle}~\in~ \mathbb{R}$$(apart from trivial cases), and not necessary. If we rewrite the wave function $\psi$ as a real column vector $$\begin{pmatrix} {\rm Re}\psi \cr {\rm Im}\psi\end{pmatrix}~\in~ \mathbb{R}^2,$$ then the functional derivative $\frac{\delta E[\psi]}{\delta\psi}$ should be interpreted as the corresponding row vector $$\left(\frac{\delta E[\psi]}{\delta{\rm Re}\psi}, \frac{\delta E[\psi]}{\delta{\rm Im}\psi} \right)~\in~ \mathbb{R}^2.$$