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

coords={xi0,xi1,xi2,xi3};

Notebook to find algebraically decomposable medium tensors with beta^2-alpha gamma = 0.

A general feature of such medium tensors seems to be that their Fresnel surfaces always have a linear factor.

Define general kappa and put X variables to zero to simplify the computations. 

Change X for speed/generality

X=10;
kappa=emGeneralKappa["k"];
For[ii=1,ii≤X,ii++,
vKappa=Variables[kappa];
v=Random[Integer,{1,Length[vKappa]}];
kappa=kappa/.vKappa[[v]]->0;
];
emKappaToMatrix[kappa]//MatrixForm

Formulate constraints:

*   kappa is algebraically decomposable
*   beta^2 -alpha = 0
*   det(kappa)=1
*   trace(kappa)=0

Abivector=emMatrix["a",4,Structure->"AntiSymmetric"];
Bbivector=emMatrix["b",4,Structure->"AntiSymmetric"];

LHS =alpha emIdentityKappa[] + beta (kappa +emPoincare[kappa]) + emCompose[emPoincare[kappa],kappa]; 
RHS=emBiProduct[rho,Abivector,Bbivector]+emBiProduct[rho,Bbivector,Abivector];
eqs=Union[Simplify[Flatten[(LHS-RHS)]]];

eqs=Append[eqs,beta^2-alpha];
eqs=Append[eqs,emDet[kappa]-1];
eqs=Append[eqs,emTrace[kappa]];

Randomly assign some variables and solve the rest.

Note: The below can take a long time.

For[done=0,done==0,(* nop *),
subs={rho->1};

(* pick values for Y variables randomly from 0/1  *)
(* Change Y for speed/generality *)
For[ii=1,ii≤20,ii++,
vars=Variables[eqs//.subs];
v=Random[Integer,{1,Length[vars]}];
v=vars[[v]];
subs=Append[subs,v->Random[Integer,{0,1}]];
];
vars=Variables[eqs//.subs];

eqs0=eqs//.subs;

sol=TimeConstrained[Solve[toEqs[eqs0],Variables[eqs0]],2];
If[Length[sol]==0,
,
Print["Found "<>ToString[Length[sol]]<>" solutions."];
done=1;
];
]
subs

index=2;

(* Compute alpha,beta,A,B,kappa for solution *)
kSubs=Join[subs,sol[[index]]];
kSubs=Append[kSubs,kx45->1];
kSubs=Append[kSubs,kx64->1/2];

a0=alpha//.kSubs;
b0=beta//.kSubs;
A0=Abivector//.kSubs;
B0=Bbivector//.kSubs;
k0=Simplify[kappa//.kSubs];

(* output *)
Variables[k0]
a0//MatrixForm
b0//MatrixForm
A0//MatrixForm
B0//MatrixForm
emKappaToMatrix[k0]//MatrixForm

fr=FullSimplify[emKappaToFresnel[k0,coords]]

Does factor!

Write equations that should be satisfied for D

(* contract kappa with a bivector from the left *)
contract[biv_,kappa_]:=Table[
1/2 Sum[
biv[[i]][[j]]emReadNormal[kappa,a,b,i,j]
,
{i,1,4},{j,1,4}
]
,
{a,1,4},{b,1,4}
]

Solve D (when possible) such that 

   D (kappa + beta Id) = 1/2 trace (rho bar(D) otimes D)  A + B.

(* contract kappa with a bivector from the left *)
contract[biv_,kappa_]:=Table[
1/2 Sum[
biv[[i]][[j]]emReadNormal[kappa,a,b,i,j]
,
{i,1,4},{j,1,4}
]
,
{a,1,4},{b,1,4}
]

Dbivector=emMatrix["d",4,Structure->"AntiSymmetric"];
dLHS = contract[Dbivector,k0+b0 emIdentityKappa[]];
dRHS = 1/2 emTrace[emBiProduct[rho,Dbivector,Dbivector]]A0+B0;

deqs=Union[Flatten[Simplify[dLHS-dRHS]]];
dsol=Solve[toEqs[deqs],Variables[deqs]];

Dbivector//.dsol[[1]]//MatrixForm