Testing Wick's theorem
Part of SNEG package, http://nrgljubljana.ijs.si/sneg
(C) Rok Zitko, rok.zitko@ijs.si, 2006, 2007

<< sneg`sneg`;

Note: a number of functions, such as number[], depend on the fact that the last index is the spin index. Therefore the channel index must come before the spin index!

snegfermionoperators[a, c, d];
snegrealconstants[x]

Off[General::spell1];

(* 
      It is important not to have interfering indeces !!! *) 
snegfreeindexes[k, k1, k2, k3, k4, 
      kp, σ, τ, δp]

    snegfreeindexes[α, α1, α2, α3, α4, αp]

    snegfreeindexes[β]

Note: O=:O:+<O>

c is SEA ordered, a and d are EMPTY ordered

ordering[c] = SEA;
ordering[c] === SEA
ordering[d] === EMPTY

Unequal operators: always sort cannonically, creation operators to the left, etc.

nc[c[AN, k, 1, UP], d[CR, k, 1, UP]]

Equal operators with EMPTY ordering

nc[a[AN, 1], a[CR, 1]]

nc[c[AN, k, 1, UP], c[CR, k, 1, DO]]

nc[c[AN, k, 1, UP], c[CR, k, 1, DO], c[AN, k1, 1, UP], 
      c[CR, k1, 1, UP]]

Empty band: creation operators shifted to the left

nc[d[AN, k, 1, sp1], d[CR, k, 1, sp2]]

Let's test improved definition of vev[], which knows more about vacuum state

vev[a[CR, UP]~nc~a[AN, UP]]
vev[a[AN, 1]~nc~a[CR, 1]]

number[c[k, α]]
vev[%]

vev[
      c[CR, k, α, σ]~nc~
        c[AN, k, α, σ]]
vev[
      c[AN, k, α, σ]~nc~
        c[CR, k, α, σ]]

vev[
      c[CR, k1, α, σ]~nc~
        c[AN, k2, α, σ]]

This one can be shown to be equal to 0: different number of spins!

nc[c[AN, k, 1, UP], c[CR, k, 1, DO], c[AN, k1, 1, UP], 
      c[CR, k1, 1, UP]]
vev[%]

nc[c[AN, k, 1, UP], c[CR, k, 1, UP], c[AN, k1, 1, UP], 
      c[CR, k1, 1, UP]]
vevwick[%]
SimplifyKD[%]

nc[c[CR, k1, 1, UP], c[CR, k2, 1, UP], c[AN, k1, 1, UP], 
      c[AN, k2, 1, UP]]
vevwick[%]
SimplifyKD[%]

vevwick[
      nc[c[CR, k1, α, σ], c[AN, k2, α, σ], 
        c[CR, k3, α, σ], 
        c[AN, k4, α, σ]]]
SimplifyKD[%]

normalorder[
      nc[c[AN, k, 1, UP], c[CR, k, 1, UP], c[AN, k1, 1, UP], 
        c[CR, k1, 1, UP]]]
vev[%]

Tests of normal ordering

nc[c[CR, k1, α, σ], 
      c[AN, k2, α, σ]]
normalorder[
      nc[c[CR, k1, α, σ], 
        c[AN, k2, α, σ]]]

Fermion fields in continuous space

psi[α_, σ_] = 
      sum[Exp[I k x] c[AN, k, α, σ], {k}]

NumericQ[x]
NumericQ[k]
NumericQ[Exp[I k x]]
Conjugate[k] ^= k
conj[Exp[I k x]]

conj[psi[α, σ]]

x1 = 
      nc[c[CR, k1, α, σ], c[CR, k2, α, σ], 
        c[AN, k3, α, σ], 
        c[AN, k4, α, σ]]

contractedpairs[1, x4]
contractedpairs[2, x4]
contractedpairs[3, x4]

contraction[c[CR, k, α, σ], 
      c[AN, kp, αp, σp]]
contraction[c[AN, k, α, σ], 
        c[CR, kp, αp, σp]] // Simplify

contract[{a[AN, σ], 
        a[CR, σ]}, {{1, 2}}]
contract[{c[CR, k, α, σ], 
        c[AN, k, α, σ]}, {{1, 2}}]
contract[{c[CR, k1, α, σ], c[AN, k2, α, σ], 
        c[CR, k3, α, σ], 
        c[AN, k4, α, σ]}, {{1, 2}}]
contract[{c[CR, k1, α, σ], c[AN, k2, α, σ], 
        c[CR, k3, α, σ], 
        c[AN, k4, α, σ]}, {{1, 2}, {3, 
          4}}]
contract[{c[CR, k1, α, σ], 
        c[AN, k2, α, σ], c[CR, k3, α, σ], 
        c[AN, k4, α, σ]}, {{1, 4}, {3, 
          2}}]

      
contract[{c[CR, k1, α, σ], c[CR, k2, α, σ], 
        c[AN, k3, α, σ], 
        c[AN, k4, α, σ]}, {{1, 3}, {2, 
          4}}]
contract[{c[CR, k1, α, σ], c[CR, k2, α, σ], 
        c[AN, k3, α, σ], 
        c[AN, k4, α, σ]}, {{1, 4}, {2, 3}}]

x1
wickorder[0, x1]
wickorder[1, x1]
wickorder[2, x1] // Simplify

nc[c[CR, k1, 1, UP], c[CR, k2, 1, UP], c[AN, k1, 1, UP], 
      c[AN, k2, 1, UP]]
wickorder[2, %] // Simplify

wick[
      nc[c[AN, k, 1, UP], c[CR, k, 1, UP], c[AN, k1, 1, UP], 
        c[CR, k1, 1, UP]]]
SimplifyKD[%]

wick[
      nc[c[CR, k1, α, σ], c[AN, k2, α, σ], 
        c[CR, k3, α, σ], 
        c[AN, k4, α, σ]]]

wick[
        nc[c[CR, k1, α, σ], c[AN, k2, α, σ], 
          c[CR, k3, α, τ], c[AN, k4, α, τ]]] // 
      Expand

y1 = 
      normalorder[
        nc[c[CR, k1, α, σ], 
          c[AN, k2, α, σ]]]
y2 = 
      normalorder[
        nc[c[CR, k3, β, τ], 
          c[AN, k4, β, τ]]]

y3 = nc[y1, y2]
y4 = vevwick[y3]
SimplifyKD[%]

y5 = Simplify[y3 - y4]

Sum[
      y5, {α, 1, 2}, {β, 1, 2}, {σ, 0, 1}, {τ, 0, 
        1}]

Simplify[%]