SetDirectory["~/Factorisation/"];
<<kappaLib.m
<<helper.m

Metaclass IV:

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

kappa=emMatrixToKappa[(({{a1, 0, 0, -b1, 0, 0}, {0, a2, 0, 0, -b2, 0}, {0, 0, a3, 0, 0, a4}, {b1, 0, 0, a1, 0, 0}, {0, b2, 0, 0, a2, 0}, {0, 0, a4, 0, 0, a3}}))];

vars={x0,x1,x2,x3};
fr=emKappaToFresnel[kappa,vars];
fr=Simplify[fr];

FullSimplify[emDet[kappa]]

If a4=0, then the Fresnel surface has 2-dim subspace, so a4 != 0.

Simplify[fr/.{a4->0,x0->0,x3->0}]

Write Fresnel polynomial using 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 b1 b2))2 ((-a1+a2) ((a1-a3) (-a2+a3)+(a4^2))- (-a2+a3) (b1^2)- (a1-a3) (b2^2)),
D1->(((((a2-a3))^2)-(a4^2)+(b2^2))/(a4 b2)),
D2->(((((a1-a3))^2)-(a4^2)+(b1^2))/(a4 b1)),
D3->(((((a1-a2))^2)+(b1^2)+(b2^2))/(b1 b2))
};

Simplify[fr-(b1 b2 a4) frSym/.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=FullSimplify[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[tmp],{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]

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

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

show[simp[gb]]

The first equation implies that D0 = 0

gb=simp[gb/.D0->0];

show[gb]

By equation (2) we have D3 = +2/-2, but since D3>=2, it follows that D3 = 2.

gb=simp[gb/.D3->2];
show[gb]

D3 = 2 implies that 

    a1 = a2      &&    b1 = b2.

We have proven: If the Fresnel surface decomposes, then

 D3 = 2, 
 a1 = a2,   
 b1 = b2.

It follows that D1 = D2

D1-D2/.subSym/.a1->a2/.b1->b2

Check that the given Lorentz metrics decompose the Fresnel surface

AA= (({{1, 0, 0, 0}, {0, (1/2) (-D1+(Sqrt[4+(D1^2)])), 0, 0}, {0, 0, (1/2) (-D1+(Sqrt[4+(D1^2)])), 0}, {0, 0, 0, -1}}));
BB= (({{1, 0, 0, 0}, {0, (1/2) (-D1-(Sqrt[4+(D1^2)])), 0, 0}, {0, 0, (1/2) (-D1-(Sqrt[4+(D1^2)])), 0}, {0, 0, 0, -1}}));

verify=(vars.AA.vars) (vars.BB.vars);
Simplify[verify-frSym/.{D0->0,D3->2,D1->D2}]

printNotebook["Metaclass_IV.pdf"]