EM-wave equation in matter from Lagrangian 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?

Expanding on Carlo Beenakker's comment, one can't just expect to substitute in a complex dielectric function to properly describe absorption. Rather, the relevant question is how to structure a Lagrangean such as to generate the desired damping term in the equation of motion. For example, to generate linear damping in a harmonic oscillator, one can use
$$
L=e^{\gamma t} \left( \frac{m\dot{q}^2 }{2} -\frac{kq^2 }{2} \right) .
$$
Since the electromagnetic field is a collection of oscillators, one might be able to use a generalization of this.
