SetDirectory["~/KappaLib/"];
<<kappaLib-1.1.m
<<helper.m

Metaclass VII:

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

kappa=emMatrixToKappa[(({{a1, 0, 0, a4, 0, 0}, {0, a2, 0, 0, a5, 0}, {0, 0, a3, 0, 0, a6}, {a4, 0, 0, a1, 0, 0}, {0, a5, 0, 0, a2, 0}, {0, 0, a6, 0, 0, a3}}))];

fr=emKappaToFresnel[kappa,vars];

We may assume that a4, a5, a6 != 0

FullSimplify[fr/.{a4->0,x2->0,x3->0}]
FullSimplify[fr/.{a5->0,x1->0,x3->0}]
FullSimplify[fr/.{a6->0,x0->0,x3->0}]

Define Fresnel equation using constants D0, D1, D2, D3

frSym=(x0^4)+(x1^4)+(x2^4)+(x3^4)+D0 x0 x1 x2 x3-D3 ((x1^2) (x2^2)+(x0^2) (x3^2))-D2 ((x1^2) (x3^2)+(x0^2) (x2^2))-D1 ((x2^2) (x3^2)+(x0^2) (x1^2));

subSym={
D0->(1/(a4 a5 a6))2 ((a1^2) (a2-a3)+(a2^2) a3+a3 (a4-a5) (a4+a5)-a2 ((a3^2)+(a4^2)-(a6^2))+a1 (-(a2^2)+(a3^2)+(a5^2)-(a6^2))),
D1->(((((a2-a3))^2)-(a5^2)-(a6^2))/(a5 a6)) ,
D2->(((((a1-a3))^2)-(a4^2)-(a6^2))/(a4 a6)),
D3->(((((a1-a2))^2)-(a4^2)-(a5^2))/(a4 a5)) 
};

Simplify[fr-(a4 a5 a6) frSym/.subSym]

Implicit expression for D0

D0sq=4 (-4+(D1^2)+(D2^2)+(D3^2)+D1 D2 D3);

Simplify[D0^2-D0sq/.subSym]

We assume that the Fresnel polynomial factorises:

A=Table[ToExpression["A"<>ToString[Min[{i,j}]]<>ToString[Max[{i,j}]]],{i,0,3},{j,0,3}];
B=Table[ToExpression["B"<>ToString[Min[{i,j}]]<>ToString[Max[{i,j}]]],{i,0,3},{j,0,3}];
A//MatrixForm
B//MatrixForm
factorised=(vars.A.vars) (vars.B.vars);

cons=Union[Flatten[CoefficientList[frSym-factorised,vars]]];

cons=simp[cons];
show[cons]//MatrixForm

Equation (2): By renaming and scaling, we may assume that A00 = 1.

sub={A00->1,B00->1};
cons=simp[cons//.sub];
show[cons]

tmp=Join[Take[cons,{5,7}],Take[cons,{17,19}],Take[cons,{23,25}]];
tmp//MatrixForm

Solve[toEqs[%],{B01,B02,B03,B11,B22,B33,B12,B13,B23}]

sub=Join[sub,%[[1]]]

cons=simp[cons//.sub];
show[cons]

Eliminate

Variables[cons]

elimVars=Variables[A]
condVars={D0,D1,D2,D3}

gb=simp[GroebnerBasis[cons,condVars,elimVars]];//Timing

show[gb]

Case 1: D0 = 0

show[simp[gb/.{D0->0}]]

It follows: If D0=0 then there exists an i in {1,2,3} such that

      D0=0,         D_i = -2 sigma,      D_i’ = sigma D_i’’

Case i=1:

dsub={D1->-sigma 2,D3->sigma D2,D0->0};

LPlus=(({{1, 0, 0, 0}, {0, sigma, 0, 0}, {0, 0, (1/2) (-D2+(Sqrt[-4+(D2^2)])), 0}, {0, 0, 0, (1/2) (-D2+(Sqrt[-4+(D2^2)])) sigma}}));LMinus=(({{1, 0, 0, 0}, {0, sigma, 0, 0}, {0, 0, (1/2) (-D2-(Sqrt[-4+(D2^2)])), 0}, {0, 0, 0, (1/2) (-D2-(Sqrt[-4+(D2^2)])) sigma}}));

dd=(vars.LPlus.vars)(vars.LMinus.vars);
simp[Flatten[CoefficientList[dd-frSym,vars]//.dsub]]

Case: i=2

dsub={D2->-sigma 2,D3->sigma D1,D0->0};

LPlus=(({{1, 0, 0, 0}, {0, (1/2) (-D1+(Sqrt[-4+(D1^2)])), 0, 0}, {0, 0, sigma, 0}, {0, 0, 0, (1/2) (-D1+(Sqrt[-4+(D1^2)])) sigma}}));LMinus=(({{1, 0, 0, 0}, {0, (1/2) (-D1-(Sqrt[-4+(D1^2)])), 0, 0}, {0, 0, sigma, 0}, {0, 0, 0, (1/2) (-D1-(Sqrt[-4+(D1^2)]))sigma}}));

dd=(vars.LPlus.vars)(vars.LMinus.vars);
simp[Flatten[CoefficientList[dd-frSym,vars]]//.dsub]

Case: i=3

dsub={D3->-sigma 2,D2->sigma D1,D0->0};

LPlus=(({{1, 0, 0, 0}, {0, (1/2) (-D1+(Sqrt[-4+(D1^2)])), 0, 0}, {0, 0, (1/2) (-D1+(Sqrt[-4+(D1^2)])) sigma, 0}, {0, 0, 0, sigma}}));LMinus=(({{1, 0, 0, 0}, {0, (1/2) (-D1-(Sqrt[-4+(D1^2)])), 0, 0}, {0, 0, (1/2) (-D1-(Sqrt[-4+(D1^2)]))sigma, 0}, {0, 0, 0, sigma}}));

dd=(vars.LPlus.vars)(vars.LMinus.vars);
simp[Flatten[CoefficientList[dd-frSym,vars]//.dsub]]

Case 2: D0 != 0

gb//show

show[simp[gb/.D1->s1 2/.D2-> s2 2]]

DB12={D1-> s1 2,D2-> s2 2,D3->(1/2) (-4 s1 s2+s3 D0)};
DB23={D2-> s1 2,D3-> s2 2,D1->(1/2) (-4 s1 s2+s3 D0)};
DB13={D1-> s1 2,D3-> s2 2,D2->(1/2) (-4 s1 s2+s3 D0)};

For[i=0,i≤1,i++,
For[j=0,j≤1,j++,
For[k=0,k≤1,k++,
loopSub={s1->(-1)^i,s2->(-1)^j,s3->(-1)^k};
t13=Simplify[D0sq/.DB13/.loopSub];
t23=Simplify[D0sq/.DB23/.loopSub];
t12=Simplify[D0sq/.DB12/.loopSub];
Print[{t13,t23,t12}]
]
]
]

Decomposition for (i,j)= (1,2)

LPlus=(({{1, 0, 0, ((Sqrt[D0])/(2 (Sqrt[2]) (Sqrt[s3])))}, {0, -s1, -(((Sqrt[D0]) (Sqrt[s3]))/(2 (Sqrt[2]))), 0}, {0, -(((Sqrt[D0]) (Sqrt[s3]))/(2 (Sqrt[2]))), -s2, 0}, {((Sqrt[D0])/(2 (Sqrt[2]) (Sqrt[s3]))), 0, 0, s1 s2}}));
LMinus=(({{1, 0, 0, -((Sqrt[D0])/(2 (Sqrt[2]) (Sqrt[s3])))}, {0, -s1, (((Sqrt[D0]) (Sqrt[s3]))/(2 (Sqrt[2]))), 0}, {0, (((Sqrt[D0]) (Sqrt[s3]))/(2 (Sqrt[2]))), -s2, 0}, {-((Sqrt[D0])/(2 (Sqrt[2]) (Sqrt[s3]))), 0, 0, s1 s2}}));
dd=(vars.LPlus.vars)(vars.LMinus.vars);
t=Union[Flatten[CoefficientList[Simplify[(dd-frSym)],vars]]];
FullSimplify[t//.DB12]

qDet = FullSimplify[Det[LPlus]]

FullSimplify[Det[LMinus]-qDet]

D13

LPlus=(({{1, 0, (((Sqrt[D0]) (Sqrt[s3]))/(2 (Sqrt[2]))), 0}, {0, -s1, 0, -((Sqrt[D0])/(2 (Sqrt[2]) (Sqrt[s3])))}, {(((Sqrt[D0]) (Sqrt[s3]))/(2 (Sqrt[2]))), 0, s1 s2, 0}, {0, -((Sqrt[D0])/(2 (Sqrt[2]) (Sqrt[s3]))), 0, -s2}}));
LMinus=(({{1, 0, -(((Sqrt[D0]) (Sqrt[s3]))/(2 (Sqrt[2]))), 0}, {0, -s1, 0, ((Sqrt[D0])/(2 (Sqrt[2]) (Sqrt[s3])))}, {-(((Sqrt[D0]) (Sqrt[s3]))/(2 (Sqrt[2]))), 0, s1 s2, 0}, {0, ((Sqrt[D0])/(2 (Sqrt[2]) (Sqrt[s3]))), 0, -s2}}));
dd=(vars.LPlus.vars)(vars.LMinus.vars);
t=Union[Flatten[CoefficientList[Simplify[(dd-frSym)],vars]]];
FullSimplify[t//.DB13]

FullSimplify[qDet-Det[LPlus]]
FullSimplify[qDet-Det[LMinus]]

D23

LPlus=(({{1, ((Sqrt[D0])/(2 (Sqrt[2]) (Sqrt[s3]))), 0, 0}, {((Sqrt[D0])/(2 (Sqrt[2]) (Sqrt[s3]))), s1 s2, 0, 0}, {0, 0, -s1, -(((Sqrt[D0]) (Sqrt[s3]))/(2 (Sqrt[2])))}, {0, 0, -(((Sqrt[D0]) (Sqrt[s3]))/(2 (Sqrt[2]))), -s2}}));
LMinus=(({{1, -((Sqrt[D0])/(2 (Sqrt[2]) (Sqrt[s3]))), 0, 0}, {-((Sqrt[D0])/(2 (Sqrt[2]) (Sqrt[s3]))), s1 s2, 0, 0}, {0, 0, -s1, (((Sqrt[D0]) (Sqrt[s3]))/(2 (Sqrt[2])))}, {0, 0, (((Sqrt[D0]) (Sqrt[s3]))/(2 (Sqrt[2]))), -s2}}));
dd=(vars.LPlus.vars)(vars.LMinus.vars);
t=Union[Flatten[CoefficientList[Simplify[(dd-frSym)],vars]]];
FullSimplify[t//.DB23]

FullSimplify[qDet-Det[LMinus]]
FullSimplify[qDet-Det[LPlus]]