It turns out to be relatively easy to define trigonometric addition rules for ReIm four-vectors. Let's start with a simple proof of the property that even makes it possible to have addition rules:

M1 = Dot[{Cos[x1], I c1 Sin[x1], I c2 Sin[x1], I c3 Sin[x1]}, 
   PauliMatrix[Range[0, 3]]];
M2 = Dot[{Cos[x2], I c1 Sin[x2], I c2 Sin[x2], I c3 Sin[x2]}, 
   PauliMatrix[Range[0, 3]]];
M3 = PolynomialMod[
   TrigReduce[M1 . M2], {c1, c2, c3} . {c1, c2, c3} - 1];
Expand@MapThread[
  Tr[M3 . PauliMatrix[#1] #2/2] &, {Range[0, 3], {1, -I, -I, -I}}]

Out[]= {Cos[x1 + x2], c1 Sin[x1 + x2], c2 Sin[x1 + x2], c3 Sin[x1 + x2]}

Now choose a system of coordinates where {1,0} falls on the unit circle, and use {1,0} as the fixed point for the addition rules. Derivation of the 4x4 matrices is programmatic:

v1 = Catenate[ReIm[M1 . {1, 0} ] /. {
     Im[x_] :> 0, Re[x_] :> x
     }];

v2 = Catenate[ReIm[M2 . {1, 0}] /. {
     Im[x_] :> 0, Re[x_] :> x
     }];

 v3 = Expand[M2 . M1 . {1, 0}];

v3b = Expand[Riffle[1/2 (v3 + ReplaceAll[v3, {I -> -I, -I -> I}]),
    -I/2 (v3 - ReplaceAll[v3, {I -> -I, -I -> I}])]];

trigAddRules = 
 Function[comp, Outer[Coefficient[comp, #1*#2] &,
    v1, v2, 1]] /@ v3b

Out[]= {
    {{1, 0, 0, 0}, {0, -1, 0, 0}, {0, 0, -1, 0}, {0, 0, 0, -1}}, 
    {{0, 1, 0, 0}, {1, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}}, 
    {{0, 0, 1, 0}, {0, 0, 0, 0}, {1, 0, 0, 0}, {0, 0, 0, 0}}, 
    {{0, 0, 0, 1}, {0, 0, 0, 0}, {0, 0, 0, 0}, {1, 0, 0, 0}}
  }

The proof is complete as follows:

v1 . # . v2 & /@ trigAddRules - v3b
Out[]= {0, 0, 0, 0}

Indeed, the first and third matrices are Pauli matrices when restricting to real numbers (odd indexed rows and columns). Maybe it's obvious that these matrices exist, but I've never seen such things computed before. Anyone else with a textbook reference?