<<VectorAnalysis`
SetCoordinates[Cartesian[x,y,z]];

Define TE and TM solutions

ee[x_,y_,z_,m_,n_] := Sin[Pi m x/a] Sin[Pi n y/b]

hh[x_,y_,z_,m_,n_]:=Cos[Pi m x/a] Cos[Pi n y/b]

kc[m_,n_]:=Sqrt[(m Pi/a)^2+(n Pi/b)^2]

ETM[x_,y_,z_,m_,n_] := ({0,0,1} ee[x,y,z,m,n] +I β/kc[m,n]^2 Grad[ee[x,y,z,m,n]]) E^(I(β z-ω t))

ETE[x_,y_,z_,m_,n_]:= -I ω Sqrt[ϵ μ] Sqrt[μ/ϵ] /kc[m,n]^2 Cross[{0,0,1}, Grad[hh[x,y,z,m,n]]]  E^(I(β z-ω t))

(* check that both ETM and ETM are solutions to the Helmholtz equation *)FullSimplify[Curl[Curl[ETM[x,y,z,m,n]]]- (kc[m,n]^2+β^2) ETM[x,y,z,m,n]]
FullSimplify[Curl[Curl[ETE[x,y,z,m,n]]]-(kc[m,n]^2+β^2) ETE[x,y,z,m,n]]

Define Beltrami fields

(* lambda = +1/-1 *)
BelTM[x_,y_,z_,m_,n_,λ_]:= 1/2(ETM[x,y,z,m,n] + λ/ Sqrt[kc[m,n]^2+β^2] Curl[ETM[x,y,z,m,n]])
BelTE[x_,y_,z_,m_,n_,λ_]:=1/2(ETE[x,y,z,m,n]+λ/Sqrt[kc[m,n]^2+β^2] Curl[ETE[x,y,z,m,n]])

kgen[m_,n_]:=Sqrt[kc[m,n]^2+β^2]

FullSimplify[kgen[m,n]BelTE[x,y,z,m,n,+1]- Curl[BelTE[x,y,z,m,n,+1]]]
FullSimplify[kgen[m,n]BelTE[x,y,z,m,n,-1]+Curl[BelTE[x,y,z,m,n,-1]]]
FullSimplify[kgen[m,n]BelTM[x,y,z,m,n,+1]- Curl[BelTM[x,y,z,m,n,+1]]]
FullSimplify[kgen[m,n]BelTM[x,y,z,m,n,-1]+Curl[BelTM[x,y,z,m,n,-1]]]

TE_11 contact forms

(* lambda = +1 or -1 *)Re[BelTE[x,y,z,m,n,λ]].Re[Curl[BelTE[x,y,z,m,n,λ]]]/.{ϵ->1, μ-> 1, a-> 1, b-> 1, β-> 1, ω-> 1};
nonZero = FullSimplify[ComplexExpand[%]]

dA =nonZero /.{m-> 1, n-> 1}

(* rescale *)
dA = dA 32 (π^2) (Sqrt[1+2 (π^2)])

Aplus =dA/. λ-> 1
Aminus =dA/. λ-> -1;
FullSimplify[Aplus+Aminus==0] 
(* the last computation shows that TE+ is a contact form if and only if TE- is a contact form *)

polP[x_,y_,z_,t_]:=2 ((Cos[π x])^2) ((Sin[π y])^2) (((Cos[t-z])^2)+(1+2 (π^2)) ((Sin[t-z])^2))+((Cos[π y])^2) (2 ((Cos[t-z])^2) ((Sin[π x])^2)+(1+6 (π^2)+(-1+2 (π^2)) Cos[2 π x]) ((Sin[t-z])^2))

Table[Plot3D[polP[x,y,2Pi k/16,0],{x,0,1},{y,0,1},PlotRange->Automatic],{k,0,16}]

(* the contact condition fails in the middle of the wave guide for both forms. This can also be seen from the graphs above. We only need to check this for polP since polM = -polP.
*)
polP[1/2,1/2,z ,t]

(* condition fails in all corners when z = 0, z = Pi and z=2Pi and t = 0 *)
Plot[{polP[0,0,z ,0],polP[0,1,z ,0]+1,polP[1,0,z ,0]+2,
polP[1,1,z ,0]+3},{z,0,2 Pi}]

Extra material: a more rigorous proof of location for zeroes of alpha wedge dalpha.

(* see derivation below *)
Alt = 2 ((Cos[π y])^2) ((Sin[π x])^2)+2 (π^2) (3+Cos[2 π x]) ((Cos[π y])^2) ((Sin[t-z])^2)+2 ((Cos[π x])^2) ((Sin[π y])^2) (((Cos[t-z])^2)+(1+2 (π^2)) ((Sin[t-z])^2));

Simplify[Alt-dA]

Plot[{(3+Cos[2 π x]),(((Cos[x])^2)+(1+2 (π^2)) ((Sin[x])^2))},{x,0,2Pi}]

(* we can write Alt as: *)
Alt = 2 ((Cos[π y])^2) ((Sin[π x])^2)+posFunc[x] ((Cos[π y])^2) ((Sin[t-z])^2)+2 ((Cos[π x])^2) ((Sin[π y])^2) posFuncII[t-z]
(* where posFunc and posFuncII are strictly positive functions *)

Derivation of expression for dA:

dA =nonZero /.{m-> 1, n-> 1}

term1=2 ((Cos[π x])^2) ((Sin[π y])^2) (((Cos[t-z])^2)+(1+2 (π^2)) ((Sin[t-z])^2));
term2=((Cos[π y])^2) (2 ((Cos[t-z])^2) ((Sin[π x])^2)+(1+6 (π^2)+(-1+2 (π^2)) Cos[2 π x]) ((Sin[t-z])^2));

Simplify[dA == term1 + term2]

term2/.{Cos[2 π x]-> Cos[Pi x]^2- Sin[Pi x]^2}

Expand[%]

%/.{((Cos[t-z])^2) -> 1- Sin[t-z]^2}

Collect[%,Sin[t-z]^2]

Simplify[(((Cos[π y])^2)+6 (π^2) ((Cos[π y])^2)-((Cos[π x])^2) ((Cos[π y])^2)+2 (π^2) ((Cos[π x])^2) ((Cos[π y])^2)-((Cos[π y])^2) ((Sin[π x])^2)-2 (π^2) ((Cos[π y])^2) ((Sin[π x])^2))]

term2Alt = 2 ((Cos[π y])^2) ((Sin[π x])^2)+2 (π^2) (3+Cos[2 π x]) ((Cos[π y])^2) ((Sin[t-z])^2);

Simplify[term2 == term2Alt]

Simplify[dA == term1 + term2Alt]

term1 + term2Alt

TM_11 plus field contact forms

(* lambda = +1 or -1 *)Re[BelTM[x,y,z,m,n,λ]].Re[Curl[BelTM[x,y,z,m,n,λ]]]/.{ϵ->1, μ-> 1, a-> 1, b-> 1, β-> 1, ω-> 1};
FullSimplify[ComplexExpand[%]]

dA =% /.{m-> 1, n-> 1}

Aplus =dA/. λ-> 1
Aminus =dA/. λ-> -1;
FullSimplify[Aplus+Aminus] 
(* the last computation shows that TE+ is a contact form if and only if TE- is a contact form *)

(* expression for lam = +1*)
pol[x_,y_,z_,t_]:=-(1/(128 (π^2)))(Sqrt[1+2 (π^2)]) (2 (-2-4 (π^2)+Cos[2 π (x-y)]+Cos[2 π (x+y)]-4 (π^2) Cos[2 (t-z)])+2 Cos[2 π x] (2 (π^2)+2 (π^2) Cos[2 (t-z)])+2 Cos[2 π y] (2 (π^2)+2 (π^2) Cos[2 (t-z)]))

Table[Plot3D[pol[x,y,2Pi k/20,0],{x,0,1},{y,0,1}],{k,0,20}]

Plot[pol[1/2,1/2,z,0],{z,0,2Pi}]

(* condition fails in all corners when z = kPi, t = 0 *)
Plot[{pol[0,0,z ,0],pol[0,1,z ,0]+1,pol[1,0,z ,0]+2,
pol[1,1,z ,0]+3},{z,0,2 Pi}]