I have recently starting using Mathematica and have recently come to what seems as an impasse in my understanding of the language. If this is too "tell me how to do it" I would certainly understand.

I have written a program that applies commutator logic to a sequence of creation and annihilation operators to produce a string of Kronecker deltas. I then need to apply the Kronecker delta restrictions to my operators. This is purely symbolic algebra, no real numbers need be produced at the end of the day.

For example,

(* KD[a,b] is a Kronecker delta with indices a and b. *)
(* KD Attributes: Orderless,NumericFunction *)
(* A[{x},{y}] is an operator with creation indices x and annihilation indices y*)

A[{x},{y}] (KD[x,a1]*KD[y,a2]-KD[x,b1]*KD[y,b2])


What I have so far:



This, however, is not quite what I need and converting from KD's to the rules list seems overly complicated. The main issues with this is that I cannot generate operator expressions that have plus within them. The structure of KD is not easily changed, but the input of the operator expressions is quite flexible.

Another take on this, I am not sure if it will help, but its my current direction.

SetAttributes[A, NumericFunciton]
A[{x},{y}] (KD[x,a1]*KD[y,a2]-KD[x,b1]*KD[y,b2])

A[{x},{y}]*KD[x,a1]*KD[y,a2]- A[{x},{y}]*KD[x,b1]*KD[y,b2]

This should be a fairly simple replacement pattern, but I cannot get any further.

  • 1
    $\begingroup$ (1) If your KD is Orderless then how do you enforce that, say, A[{x},{a2}] KD[x,a1]*KD[y,a2] will evaluate to A[{a1},{a2}] rather than A[{x},{y}] or A[{x},{a2}] or A[{a1},{y}]? (2) What is your torules[] definition? $\endgroup$ Feb 6, 2013 at 22:20
  • $\begingroup$ Have you thought about using one of the built-in deltas? Say, replace KD[x,a] with DiracDelta[x-a] and then just Integrate over {x,-Infinity,Infinity}. $\endgroup$
    – Xerxes
    Feb 6, 2013 at 22:23
  • $\begingroup$ @DanielLichtblau (1)The deltas when applied to an operator in this way is similar to saying if x replace with a1 or if a1 replace with x in the example KD[x,a1]. The formalism is designed so that you will not have these issues. (2)It does not exist, I was trying to supply at least something that I thought of. $\endgroup$
    – Ophion
    Feb 6, 2013 at 22:33
  • 1
    $\begingroup$ @Xerxes My understanding points towards these functions being used to compute numerical expressions within Mathematica. I am trying to use Mathematica to generate symbolic expressions that will be used by another program. $\endgroup$
    – Ophion
    Feb 6, 2013 at 22:35
  • $\begingroup$ Integrate is symbolic. Perhaps you are thinking of NIntegrate? Try evaluating this to see what I mean: Integrate[ A[x, y] DiracDelta[x - a] DiracDelta[y - b], {x, -Infinity, Infinity}, {y, -Infinity, Infinity}, Assumptions -> {a, b} [Element] Reals] $\endgroup$
    – Xerxes
    Feb 6, 2013 at 22:44

2 Answers 2


Here are several rules that can help to implement the Einstein convention:

deltaSimplifyRules =
   (a_ + b_)*x_ :> (a*x + b*x),
   KD[i_Symbol , j_]*
      x_ /; (MemberQ[Attributes[i], Temporary] && 
       MemberQ[x, i, Infinity]) :> (x /. i :> j),
   KD[i_ , j_Symbol]*
      x_ /; (MemberQ[Attributes[j], Temporary] && 
       MemberQ[x, j, Infinity]) :> (x /. j :> i),
   KD[i_Symbol , i_] :>  3 /; MemberQ[Attributes[i], Temporary],
   KD[i_Symbol , j_]^2 /; MemberQ[Attributes[i], Temporary] :> KD[j, j],
   KD[x__] :> Sort[KD[x]]

simplifyKD[x_] := x //. deltaSimplifyRules

now you can use simplifyKD to simplify your expression:

 Module[{x, y}, 
  A[{x}, {y}] (KD[x, a1]*KD[y, a2] - KD[x, b1]*KD[y, b2])  ]]
A[{a1}, {a2}] - A[{b1}, {b2}]

Module is used to create a temporary variable. Any temporary variable is summed over by these rules.

  • $\begingroup$ This works really well with some tweaking! One thing that I am confused about is if I set the KD string to a variable and then: simplifyKD[Module[{x, y}, A[{x}, {y}]*(KDstring)] This does not work. Its easy to work around, but interesting that its a problem. $\endgroup$
    – Ophion
    Feb 7, 2013 at 14:23

Assuming you want to normal-order products of creation and annihilation operators, followed by tidying up the Kronecker deltas and the ordering of indices, you could do something like this (using an obvious notation):

Clear[opprod, del];

Define the basic commutation rule for moving annihilation operators to the right.

opprod[u___, an[i_], cr[j_], v___] := opprod[u, cr[j], an[i], v] + del[i, j] opprod[u, v];

Tidy up a corner case.

opprod[] = id;

Define rules for simplifying and tidying Kronecker deltas - there are other possibilities that could be used here.

del[i_, i_] = 1;
del[i_, j_] /; ! OrderedQ[{i, j}] := del[j, i];
del /: del[i_, j_]^_ := del[i, j];
del /: del[i_, j_] del[j_, k_] := del[i, j] del[i, k];
del /: del[i_, k_] del[j_, k_] := del[i, j] del[i, k];

Define a rule for standardising the indices between Kronecker deltas and operator products.

opprod /: del[i_, j_] opprod[u___, (x : (an | cr))[j_], v___] :=
  del[i, j] opprod[u, x[i], v];

Define a rule for tidying up indices in operator products - this puts the indices of commuting operators in a canonical order.

opprod[u___, (x : (an | cr))[i_], (x : (an | cr))[j_], v___] /; !OrderedQ[{i, j}] :=
  opprod[u, x[j], x[i], v];

Now try out these definitions on a non-trivial operator product.

opprod[an[k], an[j], cr[i], an[i], cr[j], cr[k]] // Expand

(* id del[i, j] + id del[i, k] + 2 id del[i, j] del[i, k] + 
opprod[cr[i], an[i]] + 2 del[i, j] opprod[cr[i], an[i]] + 
2 del[i, k] opprod[cr[i], an[i]] + 
6 del[i, j] del[i, k] opprod[cr[i], an[i]] + 
del[j, k] opprod[cr[i], an[i]] + del[i, k] opprod[cr[j], an[j]] + 
del[i, j] opprod[cr[k], an[k]] + 
opprod[cr[i], cr[j], an[i], an[j]] + 
2 del[i, k] opprod[cr[i], cr[j], an[i], an[j]] + 
2 del[j, k] opprod[cr[i], cr[j], an[i], an[j]] + 
opprod[cr[i], cr[k], an[i], an[k]] + 
2 del[i, j] opprod[cr[i], cr[k], an[i], an[k]] + 
opprod[cr[i], cr[j], cr[k], an[i], an[j], an[k]] *)

I have quoted this result, not because it is necessarily interesting, but so you can verify that the above definitions are working as I had intended. I have cross-checked that this result is the same as what I get using another home-grown method that I have "up my sleeve". If it turns out to be incorrect then please tell me!

It is also very useful to use the Notation package to format the output in a human-readable way. For instance, you could try this (it looks much less messy when you enter it using the Notation package's palette):

<< Notation`

  RowBox[{"cr", "[", "i_", "]"}]]]
  RowBox[{"an", "[", "i_", "]"}]]]
InfixNotation[ParsedBoxWrapper["\[SmallCircle]"], opprod]
  RowBox[{"i_", ",", 
  RowBox[{"del", "[", RowBox[{"i_", ",", "j_"}], "]"}]]]
  • $\begingroup$ This is pretty close to what I came up with. The /: and Notation is a neat trick and will definitely clean up the code some. I am curious though why del /: del[i_, j_]^_ := del[i, j]; is necessary? $\endgroup$
    – Ophion
    Feb 7, 2013 at 14:31
  • $\begingroup$ When you use the commutator rule two or more times in succession it generates products of Kronecker deltas, so terms such as del[i,j]^2 (etc) can arise. The del /: del[i_, j_]^_ := del[i, j] definition simplifies these. $\endgroup$ Feb 7, 2013 at 17:17
  • $\begingroup$ Ah I was reading that as del /: del[i_, j_] ^= del[i, j]. One additional question, is there any benefit to the above over something like del[i_, j_]^_ := del[i, j] without the del /:? $\endgroup$
    – Ophion
    Feb 7, 2013 at 17:37
  • $\begingroup$ The effect of del /: is to make the definition associate itself with del, which is what you would want. If you omit del /: then it defaults to trying to associate itself with the Head of del[i_, j_]^_ which is Power, which is not what you you would want. Try a_^b_ := f[a,b] to see what happens. $\endgroup$ Feb 8, 2013 at 0:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.