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?

A simple observation which proves that the Lagrangian does not represent the reality of physical phynomemes, we can divide the Lagrangian into two parts, the part which contains the free field (without the presence of charges) and the part which contains the interaction part between the EM field and the material:
A simple observation that proves that the Lagrangian does not represent the reality of physical phenomena, we can divide the Lagrangian into two parts, the part that contains the free field (without the presence of charges) and the part that contains the part of the interaction between the EM field and matter :
$ \mathcal{L}=\left[\textbf{E*}(\textbf{r},\omega).(\underline{\underline{\epsilon}}-1)\textbf{E}(\textbf{r},\omega)+|\textbf{E}|^{2}(\textbf{r},\omega)-|\textbf{B}|^{2}(\textbf{r},\omega)\right]$
$\mathcal{L}=\mathcal{L}_{field}+\mathcal{L}_{interaction}$
where :
$\begin{cases}\mathcal{L}_{f}=|\textbf{E}|^{2}(\textbf{r},\omega)-|\textbf{B}|^{2}(\textbf{r},\omega)\\ \mathcal{L}_{i}=\textbf{E*}(\textbf{r},\omega).(\underline{\underline{\epsilon}}-1).\textbf{E}(\textbf{r},\omega)\end{cases}$
it is sufficient that $\underline{\underline{\epsilon}}=1$ so that it is not interaction (  $\mathcal{L}_{i}=0)$
This means that vacuum is the only medium that does not absorb or disperses EM waves for all values of the pulse, which is physically false, there are transparent media such as glass that allow light to pass for a certain pulse band.
I suggest to replace the term $(\underline{\underline{\epsilon}}-1)$ by $(\underline{\underline{\epsilon}}- \underline{\underline{\epsilon}}^{\dagger})$ in the interaction Lagrangian to at least contain the case of transparent media, i.e. in the case where  $(\underline{\underline{\epsilon}}= \underline{\underline{\epsilon}}^{\dagger})$ , we obtain a symmetrical tensor ( diagonalizable).
I do not assume that the suggested Lagrangian is the right mold in which to pour the physical optics, but it has the merit of containing more information.
