Numerically solving 3D Maxwell equations with NDEigensystem

I am trying to get electric $ \vec{E}$ and magnetic $ \vec{B}$ fields in a cylindrical cavity with a dielectric, as in the following Figure. Cavity (pink) with a dielectric inside (blue).

Both the cavity (pink) and dielectric (blue, with dielectric constant \epsilon_r) are cylindrical and share the main axis. The cavity is assumed to be conducting, such as the field at its surface has to be null.

By resorting to the vector potential $ \vec{A}$ , and using the generalized Coulomb Gauge transformation such that:

$ \vec{E}(t,\vec{r})=-\partial_t \vec{A}(t,\vec{r}),$

$ \vec{B}(t,\vec{r})=\vec{\nabla} \times \vec{A}(t,\vec{r}),$

$ \vec{\nabla} \cdot \left[ \epsilon_r(\vec{r}) \vec{A}(t,\vec{r}) \right]= 0,$

I got the system of differential equations that $ \vec{A}$ has to satisfy: \begin{equation} \vec{\nabla}^2 \vec{A}(t,\vec{r}) – \frac{1 + \epsilon_r(\vec{r})}{c^2}\partial_t^2 \vec{A}(t,\vec{r})=0. \end{equation} Here, $ \vec{r} = (x,y,z)$ .

I would like to have the eigenfrequencies and (spatial) eigenfunctions of this last operator. As far as I have understood, NDEigensystem with DirichletCondition is what I have to use; with the first output being the frequencies of the modes and the second output the spatial envelopes of the eigenfunctions.

I tried that, and failed (quite miserably). Not only the eigenfrequencies are imaginary and the spatial envelopes of the eigenfunctions do have an imaginary part$ ^{*}$ . If I plug the outcome of NDEigensystem into the differential equations, I find that these equations are not even satisfied. I am sure that it is me failing somewhere, but after long time spent trying, I am starting being really frustrated. My code is the following:

{vals, funs} = NDEigensystem[{EqSt[t, x, y, z, z1, z2, e, c] == {0, 0, 0}, BndCnd}, {Ax[t, x, y, z], Ay[t, x, y, z], Az[t, x, y, z]}, t, {x, y, z} \[Element] Cylinder[{{0, 0, 0}, {0, 0, d}}, r], 16, Method -> {"PDEDiscretization" -> {"FiniteElement", "MeshOptions" -> {"MaxCellMeasure" -> 0.01}}}}]; 

where EqSt[t, x, y, z, z1, z2, e, c] is the system of differential equations

EqSt[t_, x_, y_, z_, z1_, z2_, er_, c_] := Laplacian[{Ax[t, x, y, z], Ay[t, x, y, z], Az[t, x, y, z]}, {x, y, z}] - (1 + fer[z, z1, z2, er])/c^2 {D[Ax[t, x, y, z], {t, 2}], D[Ay[t, x, y, z], {t, 2}], D[Az[t, x, y, z], {t, 2}]}; 

and BndCnd the boundary conditions

BndCnd = {DirichletCondition[Ax[t, x, y, z] == 0, True], DirichletCondition[Ay[t, x, y, z] == ,True], DirichletCondition[Az[t, x, y, z] == 0, True]}; 

Finally, the dielectric function fer[z, z1, z2, er] that I use to mimic the dielectric is the following:

fer[z_, z1_, z2_, e_] := e (HeavisideTheta[z - z1] - HeavisideTheta[z - z2]); 

which is a step function in the z-th coordinate (= the axis of both cylinders) having value e betweem z1 and z2 and zero elsewhere.

I have tried different methods and different values of "MaxCellMeasure" (there is little improvement, but the error which I obtain by plugging the solution into the differential equations is above one!!!). Do you have any idea what is wrong here?

I have few "bonus questions" here.

-First, I tried increasing the mesh resolution for getting better outcomes, but all the times the eigenfunctions look quite bad, not smooth at all. I guess that a 3D mesh is quite demanding, but I know that the most difficult part is going to be at the dielectric. How can I tell NDEigensystem to have a finer mesh there?

-Second, there is a way to change the normalization of the eigenfunctions? I have read that NDEigensystem spits out eigenfunctions $ \vec{\phi}_i$ s.t.: $ \int \vec{\phi}_i^* \cdot \vec{\phi}_j d \vec{r} = \delta_i,j$ . I would change that to $ \int \epsilon_r(\vec{r}) \vec{\phi}_i^* \cdot \vec{\phi}_j d \vec{r} = \delta_i,j$ , would this be possible?

-Finally, I guess that it might be easier for NDEigensolver to remove the time dependence from the differential equations. This can be easily done by assuming that $ \vec{A}(t,\vec{r}) \rightarrow \vec{A}(\vec{r})e^{-i \omega t})$ . In this case, one can rewrite the operator equation above as: \begin{equation} \vec{\nabla}^2 \vec{A}(\vec{r}) + \frac{\omega^{2}}{c^2}[1 + \epsilon_r(\vec{r})] \vec{A}(\vec{r})=0. \end{equation} However, I do not know how to tell NDEigensystem that $ \omega$ is not a parameter, but should be found with the system constraints… There is a way to remove the time dependence from the equation to be given in NDEigensystem?

$ ^*$ I guess this is a stupid problem; if I change the sign of the Laplacian I get real eigenfrequencies and eigenfunctions… But I am quite confident that the sign in the above equation is correct. There is something that I am missing here?