Spectral problem for differential vector operator (calculation of EM field in a cavity)

Question

I know that mathematica has a DEigensystem and NDEigensystem which allow one to find eigenfunction and eigenvalue for scalar differential operator. But I want to find eigenfunction and eigenvalue for vector operator with additional relation between component of that wector field.

I need to solve Maxwell's equation in arbitrary shape cavity. $$ \Delta \bf{E}=-\omega^2 \bf{E}\quad \bf{\nabla}\cdot \bf{E}=0$$ With boundary conditions $\bf{E}_{\tau}=0$ in a boundary of a cavity, where $\bf{E}_{\tau}$ is a tangent vector of a cavity boundary.

I don't now how to correctly take into account the equation $\bf{\nabla}\cdot \bf{E}=0$. And how two express the boundary condition correctly?

I can express one component of vector field $\bf E$ through another two, but in this case the boundary condition take a more complex form. This boundary conditions can not be expressed via DirichletCondition and NeumannValue.

Edit:

I can express one component of the vector field through another two. For example $$E_z=-\int dz(\partial_x E_x+\partial_y E_y)+g(x,y)$$ Function $g(x,y)$ should obey the same equation $\Delta g=-\omega^2 g$.

The boundary condition can be rewrited in the following form $\bf{E}_\tau=\bf{E}-(\bf{E}\cdot \bf{n})\,n=0$ where $\bf{n}$ is a unit vector perpendicular to the boundary.

It is possible to solve it for rectangle. But it does not make sense because in such case there is a analytical solution. I don't know how to realize this boundary condition in the mathematica.

I agree, that the additional conditions $\bf{\nabla}\cdot \bf{E}=0$ it is a not big problem here. But the boundary conditions it is a big problem, at least for me.

Edit2 For simplisity I try to solve the same 2 dimentional problem. If the boundary can be set in the following form $y=g(x)$, one can calcuete the tangent vector analyticaly $\bf{\tau}=\bf{e}_x+g'(x)\bf{e}_y$. I try the realize it in the Mathematica:

gg[x_] := 1 - x^2/2
region = ImplicitRegion[x > -1 && x < 1 && y < gg[x] && y > -gg[x], {x, y}];
NDEigensystem[{-Laplacian[{ux[x, y], uy[x, y]}, {x, y}], 
DirichletCondition[uy[x, y] == 0, x == -1], 
  DirichletCondition[uy[x, y] == 0, x == 1], 
  DirichletCondition[{ux[x, y], uy[x, y]}.{1, D[gg[x], x]} == 0, 
   y > 0 && x > -1 && x < 1], 
  DirichletCondition[{ux[x, y], uy[x, y]}.{1, -D[gg[x], x]} == 0, 
  y > 0 && x > -1 && x < 1]}, {ux[x, y], 
 uy[x, y]}, {x, y} \[Element] region, 6]

But it does not work.

Also for solution inside 3D sphere

region = ImplicitRegion[x^2 + y^2 + z^2 <= 1, {x, y, z}];
g1 = (x* Ex[x, y, z] + y *Ey[x, y, z] + z* Ez[x, y, z])/(x^2 + y^2 + z^2);
NDEigensystem[{-Laplacian[{Ex[x, y, z], Ey[x, y, z], Ez[x, y, z]}, {x,
  y, z}], DirichletCondition[{Ex[x, y, z] - x*g1 == 0, 
Ey[x, y, z] - y*g1 == 0, Ez[x, y, z] - z*g1 == 0}, True]}, {Ex[x,  y, z], Ey[x, y, z], Ez[x, y, z]}, {x, y, z} \[Element] region, 6]

It also does not work.