How to construct the RHS of Wick's theorem

Question

I've been trying this for a few days now, but I just can't seem to find an efficient (speed-wise) functional way to code the following task.

Context:

I would like to construct an expressions that is the result of sorting a product of non-commuting objects (they can be expressions, but we'll use only Symbols for clarity):

input = str[j,i,h,g,f,e,d,c,b,a]

If these objects were commuting, then Sort will do the job. But I would like to write a (non-commutative sort) NCSort that will also include the appropriate commutators comm[x,y] for every required pairwise interchange. Here is the code that recursively performs this task:

NCSort[expr_str] := ReplaceRepeated[expr, 
   str[left___, a_, b_, right___] /; ! OrderedQ[{a, b}] :> 
   str[left, b, a, right] + comm[a, b] str[left, right]]

It takes 1 sec on my computer to give the result for NCSort[str[j,i,h,g,f,e,d,c,b,a]]. Furthermore, I would like the answer expanded out. Calling Expand on the output takes another 4 seconds.

Problem:

To speed things up, it is useful to know that the algebra problem can be turned into an expression-construction problem. The final result can be written in closed from (Wick's theorem),

$$ \begin{eqnarray}\text{str}(a,b,c,d,\ldots) &=& \text{sorted}(a,b,c,d,\ldots)\\ && + \hspace{-5mm}\sum_{\substack{\text{all req'd}\\\text{single comm's}}} \hspace{-5mm} \text{comm}(a,b)\text{sorted}(c,d,\ldots) + \ldots \\ && + \hspace{-5mm}\sum_{\substack{\text{all req'd}\\\text{double comm's}}} \hspace{-5mm} \text{comm}(a,b)\text{comm}(c,d)\text{sorted}(\ldots) + \ldots \\ &&\ldots \end{eqnarray}$$

I need help in directly constructing the entire expression in the RHS. Once again for clarification, this is not an algebra problem. This is an expression-construction problem. The end result should be identical to the output of NCSort[str[j,i,h,g,f,e,d,c,b,a]]//Expand given above.

---

What I have so far

Instead of using Sum I'm thinking it would be faster to write this functionally. I have pieces that would potentially be useful, but I can't quite put it all together, and that is where I need help.

This will select all possible subsets containing elements of the input string that would be involved in comm:

elementsInvolvedInCommutator[input_str, n_?EvenQ] := Subsets[List @@ input, {n}]

If n is 4 or greater, then they need to be partitioned into pairs (inspired by this answer). This involves recursion, so the spirit of the solution is lost here. I need help:

allPairwisePartitions[l : {_, _}] := {{l}}

allPairwisePartitions[list_] := 
 Join @@ Function[x, 
   Apply[{x, ##} &, 
     allPairwisePartitions[
       DeleteCases[list, Alternatives @@ x]], {1}]] /@ 
         Subsets[list, {2}, Length[list] - 1]

If I put the above two together, I can get all the pairwise commutators involved in a given sum. For example, the following gives a list of all potential double commutators (third line in the equation above):

Join @@ (allPairwisePartitions /@ elementsInvolvedInCommutator[input, 4])

I'm having trouble proceeding from here. I am unable to remove the commutators that are not needed, i.e. those that are already in order, and I need help writing allPairwisePartitions non-recursively.

Answer 1

This isn't an exact answer as for how to construct NCSort. However, it is the core of my approach for generating Wick contractions.

I ripped out all the bells and whistles.

It helps to have some background: I do lattice QCD, and so I mostly think about quark contractions. So, if you see "quark" you can parse that as "fermion". My code is designed to Wick contract a big NonCommutativeMultiply of quark operators, which are denoted quarkFlavor[spaceTimePoint][spinIndex,colorIndex] where quarkFlavor, spaceTimePoint, spinIndex and colorIndex can be whatever you want, as long as they make sense physically. In the lengthy version of this code, spin indices talk to gamma matrices, color indices talk to epsilon tensors, etc. But I've ripped out all of that in this answer because it's too much.

What really needs to be known is how to handle bar (ie. $\bar{q}$):

Clear[bar];
(* bar^2 is the identity: *)
bar[bar[f_]] := f;
(* Complex conjugate numbers: *)
bar[c_] := Conjugate[c] /; NumericQ[c];
(* Distribute over addition: *)
bar[Plus[a_, addends__]] := bar[a] + bar[Plus[addends]];
(* Distribute over multiplication: *)
bar[Times[f_, factors__]] := bar[Times[factors]] bar[f];
(* Be careful to reverse order of noncommutative multiplication: *)
bar[Times[s_., NonCommutativeMultiply[NCM___]]] := bar[s] NonCommutativeMultiply @@ bar /@ Reverse@{NCM}
(* quark flavor just bars the flavor: *)
bar[f_[x_][s_, c_]] := bar[f][x][s, c] /; Not[f === γ]

Then, I have some little helper functions:

appearanceOf[lst_][arg_] := FromDigits[
     1 - Boole@Through[(FreeQ /@ (Reverse@lst))@#], 2] &@arg;

Clear[flavor, spacetime, spin, color];
flavor[S[f_][x_][s_, c_]] := f;
flavor[f_[x_][s_, c_]] := f;
spacetime[f_[x_][s_, c_]] := x;
spin[f_[x_][s_, c_]] := s;
color[f_[x_][s_, c_]] := c;

And some utilities that help make writing big contractions easier:

(* 
 * Things for which ScalarQ are true are just numbers, in the sense that they are not anticommuting
 * once you've written out all the indices.
 *)
Clear[ScalarQ];
ScalarQ[c_?NumericQ] :=True; (* anything that Mma thinks is a number *)
ScalarQ[ϵ[idx___]] := True; (* epsilon tensors *)
ScalarQ[δ[idx___]] := True; (* kronecker deltas *)
ScalarQ[γ[_][__]] := True; (* gamma matrices *)
ScalarQ[S[__][__]] := True; (* propagators *)
ScalarQ[q[__][__]] := False;(* quarks *)

Unprotect[NonCommutativeMultiply];
Clear[NonCommutativeMultiply];
(* Flatness causes some pattern matching infinite loops.  Get rid of it: *)
ClearAttributes[NonCommutativeMultiply, Flat];
(* Pauli exclusion principle: *)
NonCommutativeMultiply[h___, a_, m___, a_, t___] := 0 
(* Distribute over addition: *)
NonCommutativeMultiply[h___, Plus[a_, addends__], t___] := NonCommutativeMultiply[h, a, t] + NonCommutativeMultiply[h, Plus[addends], t];
(* Scalars commute freely: *)
NonCommutativeMultiply[h___, Times[c_?ScalarQ, factors__], t___] :=  c NonCommutativeMultiply[h, Times[factors], t];
(* Scalars can be pulled out: *)
NonCommutativeMultiply[h___, c_?ScalarQ, t___] := c NonCommutativeMultiply[h, t]
(* One-way flatness is OK: *)
NonCommutativeMultiply[h___, NonCommutativeMultiply[NCM___], t___] := NonCommutativeMultiply[h, NCM, t];
(* The empty product is 1: *)
NonCommutativeMultiply[] := 1

And finally, the main show:

Clear[wickContract];
wickContract[Plus[a_, addends__]] := ExpandAll[wickContract[a] + wickContract[Plus[addends]]];
wickContract[Times[c_?ScalarQ, factors__]] := ExpandAll[c wickContract[Times[factors]]];
wickContract[] := 1;

wickContract[NonCommutativeMultiply[NCM___]] :=With[{
  quarks = {NCM},
  allFlavors = Union[flavor /@ {NCM} /. bar[f_] :> f]
  },
  (* Mandate that quarks appear in Overscript[q, _]q pairs: *)
  If[Not[SameQ[
    (allFlavors /. CountsBy[Select[FreeQ[bar]][quarks], flavor]),
    (bar /@ allFlavors) /. CountsBy[Select[Not@*FreeQ[bar]][quarks], flavor]]], 
    Print["Flavor mismatch."]; Return[0]];
  With[{
    sorted = Join[
      SortBy[Select[FreeQ[bar]][quarks], appearanceOf[allFlavors]],
      Reverse@SortBy[Select[Not@*FreeQ[bar]][quarks], appearanceOf[allFlavors]]
      ]},
  With[{
    qs = SplitBy[Select[FreeQ[bar]][sorted], flavor],
    qbars = Reverse@SplitBy[Select[Not@*FreeQ[bar]][sorted], flavor]
    },
    (Signature[List @@ quarks] Signature[sorted]) Times @@ MapThread[
      Function[{fs, fbars}, 
    Total[Map[
      Function[c, Signature[c] Times @@ MapThread[
        Function[{q, qbar}, S[flavor[q]][ spacetime[q] ← spacetime[qbar]][
          spin /@ {q, qbar}, color /@ {q, qbar}]], {fs, fbars[[c]]}]], 
      Permutations[Range[Length[fs]]]]]], {qs, qbars}]
]]]

Then, a proton at xf from a proton at xi looks like

wickContract[(
  u[xf][μ, a[xf]] ** u[xf][ρ[xf], b[xf]] ** d[xf][σ[xf], c[xf]])
  ** (bar[d][xi][σ[xi], c[xi]] ** bar[u][xi][ρ[xi], b[xi]] ** bar[u][xi][ν, a[xi]])
  ]

yields

S[d][xf ← xi][{σ[xf], σ[xi]}, {c[c[xi]}] (
  S[u][xf ← xi][{μ, ρ[xi]}, {a[xf], b[xi]}] 
  S[u][xf ← xi][{ρ[xf], ν}, {b[xf], a[xi]}] 
 -S[u][xf ← xi][{μ, ν}, {a[xf], a[xi]}] 
  S[u][xf ← xi][{ρ[xf], ρ[xi]}, {b[xf], b[xi]}]
  )

Where S[flavor][sink ← source][{mu,nu},{a,b}] represents a quark propagator of flavor flavor from source to sink with spin indices {mu,nu} at the sink and source and color indices {a,b} at the sink and source.

I have tested [the full version of] this code on a single proton correlator, a two-proton correlator, a proton-neutron correlator, proton-meson correlators, and a few others with correct results.

The rest of my code makes it easy to put in all the appropriate gamma matrices, epsilon tensors, etc etc etc., but that's definitely beyond the scope of this question. I would consider writing up a short note and putting the code on github with a manual on the arXiv if there's enough interest.