Section 3.1: In the paper we use that 

   G^ijkl_0 xi_i xi_j xi_k xi_l =   G^ijkl xi_i xi_j xi_k xi_l
  
for the Tamm-Rubilar tensor density G^ijkl and the corresponding non-symmetrised version.

In this notebook we check this result.

p=Permutations[{x1,x2,x3,x4}]

Length[p]

G[i_,j_,k_,l_]:=Module[
{tmp,perm,x1,x2,x3,x4},

p=Permutations[{x1,x2,x3,x4}];
p=p/.{x1->i,x2->j,x3->k,x4->l};

1/4!Sum[f[p[[s]][[1]],p[[s]][[2]],p[[s]][[3]],p[[s]][[4]]],{s,1,Length[p]}]

]

vars={xi1,xi2,xi3,xi4};

symmetricSum=Sum[G[a1,a2,a3,a4] vars[[a1]]vars[[a2]]vars[[a3]]vars[[a4]],{a1,1,4},{a2,1,4},{a3,1,4},{a4,1,4}];

nonSymmetricSum=Sum[f[a1,a2,a3,a4] vars[[a1]]vars[[a2]]vars[[a3]]vars[[a4]],{a1,1,4},{a2,1,4},{a3,1,4},{a4,1,4}];

Simplify[symmetricSum-nonSymmetricSum]