Lagrangian of Phonon-photon

A quite interesting but also hard problem are Polaritons. As far as I have understand the concept it's about phonons coupling to light. The Lagrangian function should therefore have a term for the electrostatic potential, two for the phonos (harmonic oscillator) and one for the coupling.

The Lagrange function has so basically the form

\begin{align*} L_{t}=\frac{\epsilon_{0}}{2} \int \mathrm d^{3} r(\nabla \phi)^{2}+\frac{m}{2} \mathbf{u}_{1}^{2}-\frac{f}{2} \mathbf{u}_{1}^{2} .-\rho \mathbf{u}_{1} \cdot \nabla \phi \end{align*}

where $$\mathbf u_l$$ is the relative coordinate, $$\phi$$ the potential, $$\rho$$ the charge density and $$m,\epsilon_{0}, f$$ are some constants. I want to derive an equation of motion for $$\mathbf{u}_l$$, but I struggle a lot. My attempt was \begin{align*} \frac{\mathrm{d}}{\mathrm{d} t} \frac{\partial L_{t}}{\partial \dot{\mathbf{u}}_{l}}-\frac{\partial L_{t}}{\partial \mathbf{u}_{l}}=\frac{\epsilon_{0}}{2} \int \mathrm{d}^{3} r\left[m \mathbf{ü}_{l}-f \mathbf{u}_{l}-\rho \nabla \phi\right] \end{align*} but this seems to be wrong and not solvable. It would be good to first somehow cancel out the electrostatic potential with \begin{align*} \nabla^{2} \phi=\rho / \epsilon_{0} \end{align*} It would be great if someone can explain to me how to find the equation of motion.