SetDirectory["~/writing/WIP/KappaLib/"];
<<kappaLib.m
<<helper.m

Define the medium

(* We assume kappa = C1 ast_g  + C2 Id where

g = (1,-1,-1,-1) or g=diag(-1,1,1,1)

However, since the Hodge operator is conformally
  invariant, we assume that g=(-1,1,1,1).
*)
g=DiagonalMatrix[{-1,1,1,1}];
kappa=C1 emHodge[g] + C2 emIdentityKappa[];

FullSimplify[emDet[kappa]]

Compute Fresnel polynomial

xi={x0,x1,x2,x3};

fresnel=FullSimplify[emKappaToFresnel[kappa,xi]]

Find non-zero points on the Fresnel surface of kappa.

Note: We do not need to trust the below code. We will verify the
output in the below.

coordValues={1,Sqrt[2],Sqrt[3],0};

LL=Length[coordValues];

xiList={};

For[i0=1,i0≤LL,i0++,
For[i1=1,i1≤LL,i1++,
For[i2=1,i2≤LL,i2++,
For[i3=1,i3≤LL,i3++,

frSub={
x0->coordValues[[i0]],
x1->coordValues[[i1]],
x2->coordValues[[i2]],
x3->coordValues[[i3]]};
frPoint=xi/.frSub;

(* if point is non-zero and belongs 
to the Fresnel surface add it to list. *)
If[Simplify[frPoint.frPoint]≠0,
If[Simplify[fresnel/.frSub]==0,
xiList=Append[xiList,frSub];
];
];

]
]
]
]

xiList

If 

     xi = xi_i dx^i,
     alpha = alpha_i dx^i,
     kappa = 1/8 kappa^ij_rs dx^r /\ dx^s otimes d/dx^i /\ dx^j     

then  xi /\ kappa (xi /\ alpha ) = 0 holds if and only if

     xi_i xi_a kappa^ab_lm epsilon^ilmk alpha_b = 0,     k=0,...,3

Compute the 4x4 matrix Lxi such that the above equation is 
equivalent with Lxi.alpha =0.

Lxi=Table[
Sum[
xi[[i]] xi[[a]] emReadNormal[kappa,a,b,l,m] Signature[{i,l,m,k}],
{a,1,4},{i,1,4},{l,1,4},{m,1,4}]
,
{k,1,4},{b,1,4}
];

(* matrix is symmetric *)
Union[Flatten[Lxi-Transpose[Lxi]]]

Lxi//MatrixForm

For each xi find 2 linearly independent alpha such that Lxi. alpha = 0 and alpha /\ xi != 0.

Note: We do not need to trust the below code. We will verify that lists xiList and alphaList 
have the sought properties in the next step.

alphaList={};

For[ii=1,ii≤Length[xiList],ii++,

(*** Get a point on the Fresnel surface    ***)
iSub=xiList[[ii]];
frPoint=Simplify[xi/.iSub];

(*** Find alpha such that g(frPoint,alpha)=0  ***)
aa={a0,a1,a2,a3};
evecs=Eigenvectors[(Lxi/.iSub)];
evals=Eigenvalues[(Lxi/.iSub)];

alphas = {};

For[jj=1, jj≤4,jj++,
If[evals[[jj]]==0,
ej=evecs[[jj]];

(* If eigenvector is not proportional to xi, add it to list of
possible alpha:s *)

If[Length[alphas]==0,
(* First alpha: if alpha is not proportional to xi, then it *)
propToXi=Length[Solve[ej==Const frPoint,Const]];
If[propToXi==0,
alphas=Append[alphas,ej];
];

,
(* Subsequent alphas: add unless new alpha is in the
span of old alphas and xi.
*)

consts=Table[ToExpression["Const"<>ToString[co]],{co,1,Length[alphas]+1}];

spanVectors=Append[alphas,frPoint];
inSpan=Length[Solve[ej==consts.spanVectors,consts]];
If[inSpan==0,
alphas=Append[alphas,evecs[[jj]]];
];
];
];
];

(* collect *)
alphaList =Append[alphaList, alphas];

];

Points on Fresnel surface:

xiList

For each point, a list of three alphas:

alphaList//MatrixForm

Check that xiList and alphaList satisfy the sought properties:

vecNorm[v_]:=Simplify[v.v]

frVal={};
alphasForXi={};
C123Solutions={};
isInKernel={};

For[ii=1,ii≤Length[xiList],ii++,
xiSub=xiList[[ii]];
alphas=alphaList[[ii]];

(* Check that xi is on Fresnel surface *)
(* should be zero *)
frVal=Append[frVal,Simplify[fresnel/.xiSub]];

(* check that there are two alpha:s associated to each xi. *)
(* should be 2 *)
alphasForXi=Append[alphasForXi,Length[alphas]];

(* check that alpha:s and xi are linearly independent *)
(* only solution should be C1=C2=C3=0 *)
eqs=C1 alphas[[1]]+C2 alphas[[2]]+C3 (xi/.xiSub);
C123Solutions=Append[C123Solutions,Solve[toEqs[eqs],{C1,C2,C3}]];

(* check that xi /\ kappa(xi /\ alpha) = 0 *)
(* should output (0) *)
cond ={};
For[k=1,k≤2,k++,
LL=Lxi/.xiSub;
isInKernel=Append[isInKernel,vecNorm[LL.alphas[[k]]]];

];

];

Union[frVal] (* should be 0*)
Union[alphasForXi ](* should be 3 *)
Union[C123Solutions](* should be C1=C2=C3=0 *)
Union[isInKernel] (* should be 0:s *)

Length[xiList]
Dimensions[alphaList]

Define A and B bivectors

Abivector=(({{0, A12, A13, A14}, {-A12, 0, A23, A24}, {-A13, -A23, 0, A34}, {-A14, -A24, -A34, 0}}));
Bbivector=(({{0, B12, B13, B14}, {-B12, 0, B23, B24}, {-B13, -B23, 0, B34}, {-B14, -B24, -B34, 0}}));
Abivector+Transpose[Abivector]
Bbivector+Transpose[Bbivector]

(* For bivector 'bivector' and 1-forms 'co1' and 'co2' compute

     bivector (co1 /\ co2 )

*)
contract[bivector_,co1_,co2_]:=Module[{i,j},
Sum[
bivector[[i]][[j]](co1[[i]] co2[[j]]-co1[[j]]co2[[i]]),
{i,1,4},{j,1,4}
]
]

eqs={};

For[ii=1,ii≤Length[xiList],ii++,

(*** Get a point on the Fresnel surface    ***)
iSub=xiList[[ii]];
xiVec=Simplify[xi/.iSub];

(* Go through both alpha:s associated to this xi. *)
For [k=1, k≤2, k++,

alpha=alphaList[[ii]][[k]];

ee=contract[Abivector, xiVec,alpha]contract[Bbivector,xiVec,alpha];
ee=Simplify[ee];
eqs=Append[eqs,ee];
];
];

Analysis of constraints

Length[eqs]

eqs=simp[eqs];//Timing

show[eqs]

gb=GroebnerBasis[eqs,Variables[eqs]];//Timing
gb=simp[gb];//Timing

Length[gb]

(* Routine to extract equations that depend on a given variable  *)
emEqsWithVariable[eqs_,var_]:=Module[
{ii,res,eq},
res={};

For[ii=1,ii≤Length[eqs],ii++,

eq=eqs[[ii]];
If[Count[Variables[eq],var]>0,
res=Append[res,eq];
];

];

res
]

Solve the Gröbner basis equations. We know that the Gröbner basis equations have the same solution as the original equations (in the complex domain). See for example D.Cox, J.Little, D.O’Shea, “Ideals, Varieties, and Algorithms”.

show[Take[simp[emEqsWithVariable[gb,A34]],6]]

If A34 != 0, then B=0. Thus A34 = 0.

subs={A34->0}

tmp=simp[emEqsWithVariable[gb//.subs,B12]];
tmp=Take[tmp,7];
show[tmp]

If B12 != 0, we have A=0. Thus B12 = 0

subs = Append[subs,B12->0]

tmp=simp[emEqsWithVariable[gb//.subs,B13]];
tmp=Take[tmp,7];
show[tmp]

If B13 != 0, then A = 0. We can therefore assume that B13 = 0.

subs=Append[subs,B13->0]

tmp=simp[emEqsWithVariable[gb//.subs,B23]];
tmp=Take[tmp,9];
show[tmp]

If B23 != 0, then A = 0. Thus B23 = 0

subs=Append[subs,B23->0]

tmp=simp[emEqsWithVariable[gb//.subs,B14]];
show[tmp]

If B14 != 0, then A = 0. Thus B14 = 0.

subs=Append[subs,B14->0]

tmp=simp[emEqsWithVariable[gb//.subs,B34]];
show[tmp]

If B34 != 0, then A = 0. Thus B34 = 0.

subs=Append[subs,B34->0]

show[simp[gb//.subs]]

If B24 != 0, then A = 0. Thus B24 = 0

subs=Append[subs,B24->0]

Abivector/.subs//MatrixForm
Bbivector/.subs//MatrixForm

In conclusion: Either A = 0 or B = 0.

Extra: Doublecheck result with Mathematica’s internal Solve routine

sol=Solve[toEqs[eqs],Variables[eqs]];//Timing

matrixNorm[mat_]:=Simplify[1/2Tr[mat.Transpose[mat]]]

Table[
{
matrixNorm[Abivector]//.sol[[i]],
matrixNorm[Bbivector]//.sol[[i]]
}
,
{i,1,Length[sol]}
]//MatrixForm

We get the same result: Either A=0 or B=0.