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

Suppose G^ijkl is a symmetric tensor and

    { xi in R^4 : gamma(xi)  = 0 }  supset { xi in R^4  : g(xi,xi) = 0 }
     
for gamma(xi) = G^ijkl xi_i xi_j xi_k xi_l,
     
      g = diag(-1,1,1,1).
    
Claim:  There exists a quadratic form h such that 

     gamma = g(xi,xi) h(xi,xi)      for all xi in R^4

First we define g = diag(-1,1,1,1) and gamma

g=DiagonalMatrix[{-1,1,1,1}];

symm[i_,j_,k_,l_]:=Module[{tmp},
tmp=Sort[{i,j,k,l}];
ToExpression[
"s"
<>ToString[tmp[[1]]]
<>ToString[tmp[[2]]]
<>ToString[tmp[[3]]]
<>ToString[tmp[[4]]]
]
]

symm[1,3,2,2]
symm[4,3,2,1]

(* define gamma as above *)
gamma[v_]:=Sum[
symm[i,j,k,l] v[[i]]v[[j]] v[[k]] v[[l]],
{i,1,4},
{j,1,4},
{k,1,4},
{l,1,4}
]

(* define quadratic form g *)
G[v_]:=v.g.v

gamma[{1,1,0,0}]
G[{1,1,0,0}]

Find null vectors for G.

Note: We do not need to trust the below code. We will later
check that the vectors are indeed null vectors for G.

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

LL=Length[coordValues];

nullVectors={};

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

frPoint={
coordValues[[i0]],
coordValues[[i1]],
coordValues[[i2]],
coordValues[[i3]]};

(*frPoint=vars/.frSub;*)

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

]
]
]
]

nullVectors

We found 124 null vectors for g. We can think of these as a discretisation of the null cone.

Length[nullVectors]

Verify that all vectors are indeed null vectors

Union[Table[
G[nullVectors[[i]]],
{i,1,Length[nullVectors]}
]]

Since these vectors are null for g, they are null for gamma

eqs=Table[
gamma[nullVectors[[i]]],
{i,1,Length[nullVectors]}
];

show[simp[eqs]]

We now have 62 linear equations constraining the coefficients in gamma.

sol=Solve[toEqs[eqs],Variables[eqs]]

Length[sol]

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

gammaConstrained=FullSimplify[gamma[vars]//.sol[[1]]]

HH=-(({{s1111, 2 s1112, 2 s1113, 2 s1114}, {2 s1112, s1111+6 s1122, 6 s1123, 6 s1124}, {2 s1113, 6 s1123, s1111+6 s1133, 6 s1134}, {2 s1114, 6 s1124, 6 s1134, s1111+6 s1144}}));

f1=G[vars]
f2= Simplify[vars.HH.vars]

Simplify[gammaConstrained-f1 f2]

For a similar argument for two Lorentz metrics, see: Richard A. Toupin, Elasticity and electro-magnetics, in: Non-Linear Continuum Theories, C.I.M.E. Conference, Bressanone, Italy 1965. C. Truesdell and G. Grioli coordinators. Pp.203-342.