Plot Figure 1

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

Fresnel surface

Ax =- DiagonalMatrix[{1,2,3}];
Bx = DiagonalMatrix[{1,1,1}];
Cx = 0 IdentityMatrix[3];
Dx = 0 IdentityMatrix[3];

kappa = emABCDToKappa[Ax,Bx,Cx,Dx];

xi={xi0,xi1,xi2,xi3};
fresnel = Simplify[emKappaToFresnel[kappa, xi]];

FullSimplify[fresnel]

Simplify[Solve[(fresnel/.xi0->1)==0,xi3]]

sol[xi1_,xi2_,sigma_]:=(1/(Sqrt[6]))(√(9-4 (xi1^2)-5 (xi2^2)+sigma (Sqrt[4 (xi1^4)+4 (xi1^2) (-3+(xi2^2))+(((3+(xi2^2)))^2)])))

Check that solutions parameterise the Fresnel surface

Simplify[fresnel/.{xi0->1,xi3->sol[xi1,xi2,+1]}]
Simplify[fresnel/.{xi0->1,xi3->sol[xi1,xi2,-1]}]

Define PlotRange

pRange={{-1,Sqrt[3]+1},{-1,Sqrt[3]+1},{-1,Sqrt[3]+1}};
pPoints=80;
width=0.003;

Compute curves in xz-plane

Solve[sol[x,0,1]==0,x]
Solve[sol[x,0,-1]==0,x]

P1=ParametricPlot3D[{xi1,0,sol[xi1,0,1]},
{xi1,0,Sqrt[3]},Axes->False,Boxed->False,PlotRange->pRange,PlotPoints->pPoints,PlotStyle->{Blue,Thickness[width]}];
P2=ParametricPlot3D[{xi1,0,sol[xi1,0,-1]},
{xi1,0,Sqrt[2]},Axes->False,Boxed->False,PlotPoints->pPoints,PlotStyle->{Blue,Thickness[width]}];

Compute curves in yz-plane

Solve[sol[0,y,1]==0,y]
Solve[sol[0,y,-1]==0,y]

P3=ParametricPlot3D[{0,y,sol[0,y,1]},
{y,0,Sqrt[30]},Axes->False,Boxed->True,PlotPoints->pPoints,PlotStyle->{Blue,Thickness[width]}];
P4=ParametricPlot3D[{0,y,sol[0,y,-1]},
{y,0,10},Axes->False,Boxed->False,PlotPoints->pPoints,PlotStyle->{Blue,Thickness[width]}];

Find curves in xy-plane

fresnel/.{xi3->0,xi0->1}

Solve[%==0,xi2]

P5=ParametricPlot3D[{xi1,((Sqrt[2-(xi1^2)])/(Sqrt[2])),0},
{xi1,0,Sqrt[2]},Axes->False,Boxed->True,PlotPoints->pPoints,PlotStyle->{Blue,Thickness[width]}];
P6=ParametricPlot3D[{xi1,(Sqrt[3-(xi1^2)]),0},
{xi1,0,Sqrt[3]},Axes->False,Boxed->False,PlotPoints->pPoints,PlotStyle->{Blue,Thickness[width]}];

Solve singular point

(* Find singular points on Fresnel surface*)
FrNabla = { 
D[fresnel, xi[[1]]]==0, 
D[fresnel, xi[[2]]]==0, 
D[fresnel, xi[[3]]]==0, 
D[fresnel, xi[[4]]]==0
}/.xi0->1;

FullSimplify[FrNabla] //MatrixForm

(* there are more singular points in the complex domain, but let us only solve real singular points *)
Solve[FrNabla,{xi1,xi2,xi3},Reals]

(* The fresnel equation is mirror symmetric across all xi_i planes. Let us therefore assume that xi_i >0 for all i. *)
sPoint={(Sqrt[(3/2)]),0,(1/(Sqrt[2]))};

Check that the solutions intersect at the singular point

sol[Sqrt[3/2],0,1]
sol[Sqrt[3/2],0,-1]

singPoint =Graphics3D[ {AbsolutePointSize[5.5],Point[{sPoint}]}];

Coordinate axes

eps=0.5;
back=+0.35;
w=0.04;

xArr = Graphics3D[{
Line[{{0,0,0},{Sqrt[3]+eps,0,0}}],
Line[{{Sqrt[3]+back,+w,0},{Sqrt[3]+eps,0,0}}],
Line[{{Sqrt[3]+back,-w,0},{Sqrt[3]+eps,0,0}}]
}];

yArr = Graphics3D[{
Line[{{0,0,0},{0,Sqrt[3]+eps,0}}],
Line[{{w,Sqrt[3]+back,0},{0,Sqrt[3]+eps,0}}],
Line[{{-w,Sqrt[3]+back,0},{0,Sqrt[3]+eps,0}}]
}];

zArr = Graphics3D[{
Line[{{0,0,0},{0,0,Sqrt[2]+eps}}],
Line[{{0,w,Sqrt[2]+back},{0,0,Sqrt[2]+eps}}],
Line[{{0,-w,Sqrt[2]+back},{0,0,Sqrt[2]+eps}}]
}];

Show all curves on surface

plot=Show[{P1,P2,P3,P4,P5,P6,singPoint,xArr,yArr,zArr,singPoint}]

The image in the paper was further manipulated with image manipulation softwares.