Note I am not sure if this post is of relevance for this platform, but I already asked the question in Physics Stack Exchange and in Mathematics Stack Exchange without success.
Setup Let's suppose a homogeneous dielectric with a (spatially) local dielectric response function $\underline{\underline{\epsilon}} (\omega)$ (in general a tensor), such that we have the linear response relation $$\mathbf{D}(\mathbf{x}, \omega) = \underline{\underline{\epsilon}} \left(\omega \right) \mathbf{E}(\mathbf{x}, \omega) \, ,$$ for the displacement field $\mathbf{D}$ in the dielectric. We can now write down a Lagrangian in Fourier space, describing the EM-field coupling to the dielectric body $$\mathcal{L}=\frac{1}{2}\left[\mathbf{E}^{*}\left(x,\omega\right) \cdot (\underline{\underline{\epsilon}} \left(\omega \right)-1) \mathbf{E}\left(x,\omega\right)+|\mathbf{E}|^{2}\left(x,\omega\right)-|\mathbf{B}|^{2}\left(x,\omega\right) \right] \, .$$ If we choose a gauge \begin{align} \mathbf{E} &= \frac{i \omega}{c} \mathbf{A} \\ \mathbf{B} &= \nabla \times \mathbf{A} \, , \end{align} such that we can write the Lagrangian (suppressing arguments) in terms of the vector potential $\mathbf{A}$ as $$\mathcal{L} =\frac{1}{2}\left[\frac{\omega^2}{c^2} \mathbf{A}^{*} \cdot (\underline{\underline{\epsilon}}- \mathbb{1}) \mathbf{A}+ \frac{\omega^2}{c^2} |\mathbf{A}|^{2}-|\nabla \times \mathbf{A}|^{2}\right] \, . $$ And consequently we have the physical action $$ S[\mathbf{A}] = \int d \mathbf{x} \int \frac{d \omega}{2 \pi} \; \mathcal{L} \left(\mathbf{A}\right) \, . $$
Goal My goal is to derive the EM-wave equation for the electric field in the dielectric media.
Idea So my ansatz is the following: If we use Hamilton's principle, we want the first variation of the action to be zero \begin{align} 0 = \delta S[\mathbf{A}] &= \left.\frac{\mathrm{d}}{\mathrm{d} \varepsilon} S[\mathbf{A} + \varepsilon \mathbf{h}] \right|_{\varepsilon=0} \\ &= \left.\frac{\mathrm{d}}{\mathrm{d} \varepsilon} \int d \mathbf{x} \int \frac{d \omega}{2 \pi} \; \mathcal{L} (\mathbf{A} + \varepsilon \mathbf{h}) \right|_{\varepsilon=0} \\ &= \int d \mathbf{x} \int \frac{d \omega}{2 \pi} \; \frac{1}{2} \Bigg( \frac{\omega^2}{c^2} \mathbf{A}^* \cdot ({\underline{\underline{\epsilon}}}-\mathbb{1}) \mathbf{h} + \frac{\omega^2}{c^2} \mathbf{h}^* \cdot ({\underline{\underline{\epsilon}}}- \mathbb{1}) \mathbf{A} + \frac{\omega^2}{c^2} \mathbf{A}^* \cdot \mathbf{h} + \frac{\omega^2}{c^2} \mathbf{h}^* \cdot \mathbf{A} \\ &\quad \quad \quad \quad \quad \quad \quad \quad- (\nabla \times \mathbf{A}^* ) \cdot ( \nabla \times \mathbf{h}) - (\nabla \times \mathbf{h}^* ) \cdot ( \nabla \times \mathbf{A}) \Bigg) \\ &= \int d \mathbf{x} \int \frac{d \omega}{2 \pi} \; \frac{1}{2} \Bigg( \frac{\omega^2}{c^2} \left[ ({\underline{\underline{\epsilon}}}^{\dagger}-\mathbb{1}) \mathbf{A} \right]^* \cdot \mathbf{h} + \frac{\omega^2}{c^2} \left[({\underline{\underline{\epsilon}}}- \mathbb{1}) \mathbf{A} \right] \cdot \mathbf{h}^* + \frac{\omega^2}{c^2} \mathbf{A}^* \cdot \mathbf{h} + \frac{\omega^2}{c^2} \mathbf{A} \cdot \mathbf{h}^* \\ &\quad \quad \quad \quad \quad \quad \quad \quad- (\nabla \times \nabla \times \mathbf{A}^* ) \cdot \mathbf{h} - (\nabla \times \nabla \times \mathbf{A} ) \cdot \mathbf{h}^* \Bigg) \\ &= \int d \mathbf{x} \int \frac{d \omega}{2 \pi} \; \frac{1}{2} \Bigg( \underbrace{\left[ \frac{\omega^2}{c^2} \left[ ({\underline{\underline{\epsilon}}}^{\dagger}-\mathbb{1}) \mathbf{A} \right]^* + \frac{\omega^2}{c^2} \mathbf{A}^* - \nabla \times \nabla \times \mathbf{A}^* \right]}_{\stackrel{!}{=} 0} \cdot \mathbf{h} \\ & \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad + \underbrace{\left[ \frac{\omega^2}{c^2} \left[({\underline{\underline{\epsilon}}}- \mathbb{1}) \mathbf{A} \right] + \frac{\omega^2}{c^2} \mathbf{A} - \nabla \times \nabla \times \mathbf{A} \right]}_{\stackrel{!}{=} 0} \cdot \mathbf{h}^* \Bigg) \, , \end{align} for all $\mathbf{h}(\mathbf{x}, \omega)$. And consequently we get the equations \begin{align} \frac{\omega^2}{c^2} \left[ ({\underline{\underline{\epsilon}}}^{\dagger}-\mathbb{1}) \mathbf{A} \right]^* + \frac{\omega^2}{c^2} \mathbf{A}^* - \nabla \times \nabla \times \mathbf{A}^* &= 0 \\ \frac{\omega^2}{c^2} \left[({\underline{\underline{\epsilon}}}- \mathbb{1}) \mathbf{A} \right] + \frac{\omega^2}{c^2} \mathbf{A} - \nabla \times \nabla \times \mathbf{A} &= 0 \, . \end{align} If we suppose a lossy dielectric body, such that $\underline{\underline{\epsilon}}^{\dagger} \neq \underline{\underline{\epsilon}}$, the equations are in contradiction.
Time-domain An analogues derivation in the time-domain (which I can post here on request), yields the wave equation (in Fourier space) $$\frac{\omega^2}{2 c^2} (\underline{\underline{\epsilon}}-1) \mathbf{A} + \frac{\omega^2}{2 c^2} (\underline{\underline{\epsilon}}^{\dagger} -1) \mathbf{A} + \frac{\omega^2}{c^2} \mathbf{A} - \left( \nabla \times \nabla \times \mathbf{A} \right) = 0 \, .$$ This result is also not resembling the expected result, for $\underline{\underline{\epsilon}}^{\dagger} \neq \underline{\underline{\epsilon}}$.
Question What went wrong in the calculation?