I'd like to ask for a clarification regarding the Taylor expansion for the potential energy in the case of a 3D harmonic crystal. Let $\vec{R}\in\mathbb{R}^3$ be a Bravais lattice vector, pointing to the equilibrium position of an atom in the crystal and let $\vec{u}(\vec{R},t)$ be the displacement from equilibrium so that we may write: $$\vec{r} = \vec{R} + \vec{u}(\vec{R},t)$$ as the time-dependent atomic position.
Further, suppose that the interaction potential can be written as sum of pair interactions (where $\Phi:\mathbb{R}^3\rightarrow\mathbb{R}$ is the pair potential) so that the total potential energy is given by: $$U = \frac{1}{2}\sum_{\vec{r},\vec{r}'} \Phi \left(\vec{r}-\vec{r}'\right)$$ If we rewrite the potential as: $$U = \frac{1}{2}\sum_{\vec{R},\vec{R}'} \Phi \left(\vec{R}-\vec{R}' + \vec{u}(\vec{R})-\vec{u}(\vec{R}')\right)$$ the Taylor expansion in power series of $\epsilon(\vec{R},\vec{R}') = \vec{u}(\vec{R})-\vec{u}(\vec{R}')$ up to second order is: $$U = \frac{1}{2}\sum_{\vec{R},\vec{R}'} \left(\Phi(\vec{R}-\vec{R}')\right) + \frac{1}{2}\sum_{\vec{R},\vec{R}'} (\vec{\epsilon}\cdot\vec{\nabla})\Phi(\vec{R}-\vec{R}') + \frac{1}{4}\sum_{\vec{R},\vec{R}'} (\vec{\epsilon}\cdot\vec{\nabla})^2\Phi(\vec{R}-\vec{R}')$$ Since each term is evaluated at equilibrium ($\vec{\epsilon}=\vec{0}$), the first order term vanishes since it represents the total force acting on the atom in $\vec{R}$ due to all the others. My problem arises when I try to write the second order term in a more compact manner. I proceed as follows: $$ (\vec{\epsilon}\cdot\vec{\nabla})^2\Phi = \begin{pmatrix}\epsilon_1 & \epsilon_2 & \epsilon_3\end{pmatrix}\cdot\begin{pmatrix}\frac{\partial^2\Phi}{\partial{(r_1-r_1')}^2}&\frac{\partial^2\Phi}{\partial{(r_2-r_2')}\partial{(r_1-r_1')}}&\frac{\partial^2\Phi}{\partial{(r_3-r_3')}\partial{(r_1-r_1')}}\\\ \frac{\partial^2\Phi}{\partial{(r_1-r_1')}\partial{(r_2-r_2')}}&\frac{\partial^2\Phi}{\partial{(r_2-r_2')}^2}&\frac{\partial^2\Phi}{\partial{(r_3-r_3')}\partial{(r_2-r_2')}}\\\ \frac{\partial^2\Phi}{\partial{(r_1-r_1')}\partial{(r_3-r_3')}}&\frac{\partial^2\Phi}{\partial{(r_2-r_2')}\partial{(r_3-r_3')}}&\frac{\partial^2\Phi}{\partial{(r_3-r_3')}^2} \end{pmatrix}_{\vec{\epsilon}=\vec{0}} \cdot \begin{pmatrix}\epsilon_1\\\ \epsilon_2 \\\ \epsilon_3\end{pmatrix}$$ where $\epsilon_i = u_i(\vec{R}) - u_i(\vec{R}')$.
In a more compact form this yields: $$ \sum_{\mu,\nu} (u_{\mu}(\vec{R}) - u_{\mu}(\vec{R}'))\frac{\partial^2\Phi(\vec{R}-\vec{R}')}{\partial{(r_{\mu} - r_{\mu}')}\partial{(r_{\nu}-r_{\nu}')}}(u_{\nu}(\vec{R}) - u_{\nu}(\vec{R}')) $$
However, in the lectures this was instead written as: $$ \sum_{\mu,\nu} (u_{\mu}(\vec{R}) - u_{\mu}(\vec{R}'))\frac{\partial^2\Phi(\vec{R}-\vec{R}')}{\partial{u_{\mu}(\vec{R})}\partial{u_{\nu}(\vec{R}')}}(u_{\nu}(\vec{R}) - u_{\nu}(\vec{R}')) $$ and no further explanation was given.
I find myself unable to see why (or if) these two are the same. Do you have any ideas? Any help would be much appreciated.