Contracting indices in lists

index[_[x__]] := x (* <-- Extracts arguments from an expression. *)

muteIndexSum[list_] := Module[
  {indices, repeatedIndices, result},
  (*ALL indices: *)
  indices = Table[index[gamma], {gamma, list}];
  (*REPEATED indices: *)
  repeatedIndices = {};
  Do[If[Count[indices, i] == 2 && ContainsNone[repeatedIndices, {i}], 
    AppendTo[repeatedIndices, i]], {i, indices}];
  (*Dot product of gamma matrices: *)
  result = Apply[Dot, list];
  Do[result = Sum[result /. i -> k, {k, 0, 3}], {i, 
    repeatedIndices}];
  (* If there are no repeated indices, 
  this will just return the dot product of what you passed. *)

  Return[result];
  ]

Since you have already used γ, I am using γm as a notation (using γ would mess things up):

muteIndexSum[{γm[μ], γm[ν], γm[μ]}]

Output:

γm[0].γm[ν].γm[0] + γm[1].γm[ν].γm[1] + γm[2].γm[ν].γm[2] + γm[3].γm[ν].γm[3]

If you contract all indices,

muteIndexSum[{γm[μ], γm[ν], γm[μ], γm[ν]}]

The output will be:

γm[0].γm[0].γm[0].γm[0] + γm[0].γm[1].γm[0].γm[1] + γm[0].γm[2].γm[0].γm[2] + γm[0].γm[3].γm[0].γm[3] + γm[1].γm[0].γm[1].γm[0] + γm[1].γm[1].γm[1].γm[1] + γm[1].γm[2].γm[1].γm[2] + γm[1].γm[3].γm[1].γm[3] + γm[2].γm[0].γm[2].γm[0] + γm[2].γm[1].γm[2].γm[1] + γm[2].γm[2].γm[2].γm[2] + γm[2].γm[3].γm[2].γm[3] + γm[3].γm[0].γm[3].γm[0] + γm[3].γm[1].γm[3].γm[1] + γm[3].γm[2].γm[3].γm[2] + γm[3].γm[3].γm[3].γm[3]

If you want to use your pre-defined γ matrices, just use the rule to substitute:

muteIndexSum[{γm[μ], γm[ν], γm[μ], γm[ν]}]/.{γm -> γ}

Output:

{{4, 0, 0, 0}, {0, 4, 0, 0}, {0, 0, 4, 0}, {0, 0, 0, 4}}

What the shape of $ \gamma_\mu \gamma^\mu $ do you expect? I am not sure whether the code below works to your satisfaction or not:

Total @ MapThread[Dot, {γ[Range[0, 3]], γ[Range[0, 3]]}]

---

Update

OK, let me make it in a more formal and natural way. First let some points made clear:

• Contraction of a pair of two same Greek indexes is accompanied with Minkowski metric signatured $ \mathrm{diag} g_{\mu\nu} = \{+, -, -, -\} $ (g = DiagonalMatrix[SparseArray @ {1, -1, -1, -1}];);
• Contraction of a pair of two same Latin indexes is the same as that of a matrix dot product, or say the metric is just the identity matrix;
• Einstein convention is employed.

Then the quantity of interest in fact is

$$ \gamma^\mu \gamma_\mu = (\gamma^\mu)^{mn} (\gamma_\mu)_{nl} = (\gamma^{^{\substack{\color{red}{1} \\ \mu}}})^{^{\substack{2\; \color{blue}{3} \\ mn}}} g_{_{\substack{\mu\nu \\ \color{red}{4}\,\color{cyan}{5}}}} (\gamma^{^{\substack{\color{cyan}{6} \\ \nu}}})_{_{\substack{nl \\ \color{blue}{7}8}}} $$

So it should work by using TensorProduct with TensorContract (pay attention to which two indexes are paired to contract):

γμ = γ[Range[0, 3]];
TensorContract[TensorProduct[γμ, g, γμ], {{1, 4}, {5, 6}, {3, 7}}]

In the same spirit $$ \sigma^{\mu\nu}\sigma_{\nu\rho} = (\sigma^{\mu\nu})^{mn}\ g_{\nu\alpha}\ (\sigma^{\alpha\beta})_{nl}\ g_{\beta\rho} \\ \sigma^{\mu\nu}\sigma_{\mu\nu} = (\sigma^{\mu\nu})^{mn}\ g_{\mu\alpha}\ (\sigma^{\alpha\beta})_{nl}\ g_{\beta\nu} $$ should correspond to, respectively,

Clear[σ]
σ[μ_, ν_] := I (γ[μ].γ[ν] - γ[ν].γ[μ]) / 4
σμν = Outer[σ, Range[0, 3], Range[0, 3]] // SparseArray;
TensorContract[TensorProduct[σμν, g, σμν, g], {{2, 5}, {4, 9}, {6, 7}, {8, 11}}]
TensorContract[TensorProduct[σμν, g, σμν, g], {{1, 5}, {2, 12}, {4, 9}, {6, 7}, {8, 11}}]

Other contractions can be done similarly.