How to find an algebraic representation

Question

I want to find matrix representation of some algebra (Clifford algebra Cl(3,1) in this case). Here is an example, which I would like to get and then extend to higher dimensional matrices.

Suppose we look for 4x4 matrix representation and want to find 4 matrices (matrix dimensions are not directly related to number of matrices)

em1 = Array[a1, {4, 4}]; em2 = Array[a2, {4, 4}]; 
em3 =  Array[a3, {4, 4}]; em4 = Array[a4, {4, 4}];

I know that the algebra representation can be found using SolveAlways as follows

SolveAlways[
{Det[c[1]*em1 + c[2]*em2 + c[3]*em3 + 
             c[4]*em4] == (c[1]^2 + c[2]^2 + c[3]^2 - c[4]^2)^2, 
          Det[em1] == 1, Tr[em1] == 0, Det[em2] == 1, Tr[em2] == 0, 
          Det[em3] == 1, Tr[em3] == 0, Det[em4] == 1, Tr[em4] == 0}, 
 {c[1],c[2], c[3], c[4]}]

Unfortunatelly, this runs forever (>2 days on my laptop). Below is the test that at least one of representation, which satisfy the the above conditions is

em1test = {{0, 0, 1, 0},
      {0, 0, 0, 1},
      {1, 0, 0, 0},
      {0, 1, 0, 0}}; 
em2test = {{0, 0, -I, 0},
      {0, 0, 0, -I},
      {I, 0, 0, 0},
      {0, I, 0, 0}}; 
em3test = {{0, -I, 0, 0},
      {I, 0, 0, 0},
      {0, 0, 0, I},
      {0, 0, -I, 0}}; 
em4test = {{I, 0, 0, 0},
      {0, -I, 0, 0},
      {0, 0, -I, 0},
      {0, 0, 0, I}};

Indeed,

testRules = 
 Join[Thread[Flatten[em1] -> Flatten[em1test]], 
  Thread[Flatten[em2] -> Flatten[em2test]], 
  Thread[Flatten[em3] -> Flatten[em3test]], 
  Thread[Flatten[em4] -> Flatten[em4test]]]

and

Det[c[1]*em1 + c[2]*em2 + c[3]*em3 + c[4]*em4] /. testRules // Simplify
(\* (c[1]^2 + c[2]^2 + c[3]^2 - c[4]^2)^2 \*)

Unfortunatelly separating real and imaginary parts seems do not reduce complexity of the problem (same number of unknowns).

From documentation I know that SolveAlways can be rewriten in terms of, say Solve and Eliminate. So, I hope, that it can be rewriten in terms of GroebnerBasis itself. How can this be done? How far I can hope to go for looking for representations of higher algebras? Can I hope say, 8x8 dimensional representation for 7 matrices to be reachable within this approach? Any other ideas?

Edit 1

Answering to the comment I can reformulate the problem in terms of hermitian Transpose[Conjugate[m]]=m and antihermitian Transpose[Conjugate[m]]=-m matrices as follows:

(reProto={{p[1,1],p[1,2],p[1,3],p[1,4]},{p[1,2],p[2,2],p[2,3],p[2,4]},{p[1,3],p[2,3],p[3,3],p[3,4]},{p[1,4],p[2,4],p[3,4],p[4,4]}})//MatrixForm

(imProto={{0,i[1,2],i[1,3],i[1,4]},{-i[1,2],0,i[2,3],i[2,4]},{-i[1,3],-i[2,3],0,i[3,4]},{-i[1,4],-i[2,4],-i[3,4],0}})//MatrixForm

and then

em1=(reProto+I*imProto)/.{p->ar1,i->ai1};
em2=(reProto+I*imProto)/.{p->ar2,i->ai2};
em3=(reProto+I*imProto)/.{p->ar3,i->ai3};
em4=(I*reProto+imProto)/.{p->ar4,i->ai4};

Unfortunatelly, counting variables we see that number of unknowns did't change, i.e. 64+4, except that now solutions should be real numbers. Unfortunatelly I don't know how to use this restriction in SolveAlways.

Edit2

Answering the comment about restructions I reformulated the problem to include all known constraints

 SolveAlways[Flatten[{Det[c[1]*em1+c[2]*em2+c[3]*em3+c[4]*em4]==(c[1]^2+c[2]^2+c[3]^2-c[4]^2)^2,
Tr[em1]==0,Tr[em2]==0,Tr[em3]==0,Tr[em4]==0,
    Thread[Flatten[em1.em2+em2.em1]==Flatten[Table[0,{4},{4}]]],
    Thread[Flatten[em1.em3+em3.em1]==Flatten[Table[0,{4},{4}]]],
    Thread[Flatten[em2.em3+em3.em2]==Flatten[Table[0,{4},{4}]]],
    Thread[Flatten[em1.em4+em4.em1]==Flatten[Table[0,{4},{4}]]],
    Thread[Flatten[em2.em4+em4.em2]==Flatten[Table[0,{4},{4}]]],
    Thread[Flatten[em3.em4+em4.em3]==Flatten[Table[0,{4},{4}]]],
    Thread[Flatten[em1.em1+em1.em1]==Flatten[2IdentityMatrix[4]]],
    Thread[Flatten[em2.em2+em2.em2]==Flatten[2IdentityMatrix[4]]],
    Thread[Flatten[em3.em3+em3.em3]==Flatten[2IdentityMatrix[4]]],
    Thread[Flatten[em4.em4+em4.em4]==Flatten[-2IdentityMatrix[4]]],
    Det[em1]==1,Det[em2]==1,Det[em3]==1,Det[em4]==1}],
        {c[1],c[2],c[3],c[4]}]