# How to impose null divergence in the solution of NDEigensystem

I am trying to use NDEigensystem to find eigenvalues and eigenmodes of the following Laplacian with boundary conditions:

$$\vec{\nabla}^2 \vec{A}(x,z) = vals \; \vec{A},$$

$$A_{x}(x,z=\pm \frac{d}{2}) = 0,$$

$$A_{y}(x,z=\pm \frac{d}{2}) = A_{y}(x = \pm r,z) = 0,$$

$$A_{z}(x = \pm r,z) = 0,$$

$$\frac{\partial A_{x}(x,z = \pm d/2)}{\partial x} = 0,$$

$$\frac{\partial A_{z}(x = \pm r,z)}{\partial z} = 0,$$

where $$\vec{A} = \{A_x,A_y,A_z\}$$ and $$vals$$ is an eigenvalue found by NDEigensystem. The domain is a rectangle centred at $$(x,z) = (0,0)$$ and with basis $$2r$$ and height $$d$$. Importantly, to get these equations I used the fact that

$$\nabla \cdot \vec{A} = 0,$$

and as such I would like the solutions to respect this constraint. However, for (some) of the solutions found by NDEigensystem, $$\nabla \cdot \vec{A} \neq 0$$ (as you can see in the last plot if you ran the code below). Luckily, these solutions come in groups with the same eigenvalue, and a superposition of them actually can give null divergence. I am confident that I am using NeumannValue in the wrong way. I also tried to impose the null divergence directly in the operator by setting $$\partial_x A_x = -\partial_z A_z$$, but in that case NDEigensystem was giving a completely crazy solution. In particular, the boundary conditions involving the derivatives were not satisfied anymore (this further suggests that there is a problem with NeumannValue - or better with me trying to use NumannValue).

My question is the following: How can I impose null divergence of the solution? How does NeumannValue works, exactly, in case of coupled differential equations? I have tried to read documentation about NeumannValue, but unfortunately I am a lame physicists and have really not understood if and how it can be used to impose null divergence of a field in my case or in general.

As a small comment, I am aware that my problem can be solved analytically in a relatively simple manner. The reason for which I want to have a numerical solution is that I would like to increase the complexity once I am sure that I understood how NDEigensystem and NeumannValue work.

I have a "bonus question" here:

• There is a way to change the normalization of the eigenfunctions? I have read that NDEigensystem spits out eigenfunctions $$\vec{A}_i$$ s.t.: $$\int \vec{A}_i^* \cdot \vec{A}_j dxdz = \delta_{i,j}$$. I would like to change that to $$\int \epsilon_r(x,z) \vec{A}_i^* \cdot \vec{A}_j dxdz = \delta_{i,j}$$, where $$\epsilon_r(x,z)$$ is a certain function. I can do that by hand afterwards, but it seems that there is a method option in NDEigensystem that does that (called VectorNormalization), but I did not understand how it works...

In the following, the block of the code that I am currently using:

(*Parameters:*)
d = 3;
r = 1;
(*Domain:*)
\[CapitalOmega] = Rectangle[{-r, -(d/2)}, {r, d/2}];
Graphics[{Opacity[.3], Blue, \[CapitalOmega]}, Axes -> True,
AxesLabel -> {"x", "z"}]
(*Differential operator:*)
EqSt[x_, y_, z_] :=
Laplacian[{Ax[x, z], Ay[x, z], Az[x, z]}, {x, y, z}, "Cartesian"]
(*Dirichlet and Neumann boundary conditions:*)
BndCnd = {DirichletCondition[Ax[x, z] == 0, z^2 == (d/2)^2],
DirichletCondition[Ay[x, z] == 0, Or[x^2 == r^2, z^2 == (d/2)^2]],
DirichletCondition[Az[x, z] == 0, x^2 == r^2]};
NeumBnd = {NeumannValue[0, x^2 == r^2], 0,
NeumannValue[0, z^2 == (d/2)^2]};
(*Getting solutions:*)
{valsZ3, funsZ3} =
NDEigensystem[
Flatten[{EqSt[x, y, z] - NeumBnd, BndCnd}], {Ax, Ay,
Az}, {x, z} \[Element] \[CapitalOmega], 16];
(*Eigenvalues:*)
valsZ3
(*Checling the solutions:*)
ii = 3;
{Plot3D[funsZ3[[ii, 1]][x, z], {x, -r, r}, {z, -(d/2), d/2},
PlotRange -> All, AxesLabel -> {"x", "z"}, PlotLabel -> "Ax[x,z]"],
Plot3D[funsZ3[[ii, 2]][x, z], {x, -r, r}, {z, -(d/2), d/2},
PlotRange -> All, AxesLabel -> {"x", "z"}, PlotLabel -> "Ay[x,z]"],
Plot3D[funsZ3[[ii, 3]][x, z], {x, -r, r}, {z, -(d/2), d/2},
PlotRange -> All, AxesLabel -> {"x", "z"}, PlotLabel -> "Az[x,z]"],
Plot3D[Evaluate[
Div[{funsZ3[[ii, 1]][x, z], funsZ3[[ii, 2]][x, z],
funsZ3[[ii, 3]][x, z]}, {x, y, z}]], {x, -r, r}, {z, -(d/2), d/
2}, PlotRange -> All, AxesLabel -> {"x", "z"},
PlotLabel ->
"\!$$\*SubscriptBox[\(\[Del]$$, $${x, y, \ z}$$]\).{Ax[x,z],Ay[x,z],Ax[x,z]}",
ColorFunction -> Function[{x, y, z}, Hue[.65 (1 - z)]]]}