I am trying to do some products involving some objects made out of Kronecker deltas. For example, taking an object like $x_{abcd}=\delta_{ab}\delta_{cd}$, where all the indices run from $1$ to $N$, I would like to compute a product like
$$x_{abmn}x_{mncd},$$ where repeated indices are being summed over. Is there a way in Mathematica to fully compute such a product? I mean, it is quite easy to just get the a symbolic product like $$x_{abmn}x_{mncd}=\delta_{ab}\delta_{mn}\delta_{mn}\delta_{cd}.$$ But I would like to further carry out the sums over $m$ and $n$ to obtain $$x_{abmn}x_{mncd}=\delta_{ab}\delta_{cd}N.$$ Of course this is a very simple and quick example, but eventually I would calculate products involving hundreds of terms so any tips would be appreciated!

Edit: Here is an example code of just the symbolic product (What I have manage to do). There is not much to it, it is just a product, and it produces the right result.


The rest I do not know how to do. What I want the code to do is the following $$x_{abmn}x_{mncd}=\delta_{ab}\delta_{mn}\delta_{mn}\delta_{cd}=\sum_{n=1}^N\sum_{m=1}^N\delta_{ab}\delta_{mn}\delta_{mn}\delta_{cd}=\sum_{m=1}^N\delta_{ab}\delta_{mm}\delta_{cd}=\delta_{ab}\delta_{cd}N$$

  • 2
    $\begingroup$ Please include an example in MMA syntax, together with the expected output. $\endgroup$
    – MarcoB
    Jul 11, 2023 at 21:04

3 Answers 3


Here's one way, using pattern matching. (I used nsize instead of N, since N is a reserved symbol). First define x, as you did:

x[a_, b_, c_, d_] := KroneckerDelta[a, b] KroneckerDelta[c, d];

Next define "contract using Einstein Summation" - i.e. sum on repeated indices. There are three rules, one for two different deltas with at least one repeated index, a second for powers, since we get delta^2 terms, and one to terminate the recursion:

   Times[x___, KroneckerDelta[a_, b_], y___, KroneckerDelta[c_, b_], 
    z___]] := nsize contractES[Times[KroneckerDelta[a, c], x, y, z]];

  Times[x___, KroneckerDelta[a_, b_]^n_?(# > 1 &), y___]] := 
 nsize^2 contractES[Times[x, y, KroneckerDelta[a, b]^(n - 2)]];

contractES[x___] := x;

KroneckerDelta has the Orderless Attribute, so these patterns suffice even if the indices are in a different order in the input.

Now runs tests, first for x:

x[a, b, m, n] x[m, n, c, d]

Out[7]= KroneckerDelta[a, b] KroneckerDelta[c, d] KroneckerDelta[m, 

And now for your example:

contractES[x[a, b, m, n] x[m, n, c, d]]

Out[8]= nsize^2 KroneckerDelta[a, b] KroneckerDelta[c, d]

I get nsize^2 because we get one factor of nsize for each index, and we contract over both m and n. So my nsize^2 is basically what you're calling "N".

Finally do a couple more cases, to show that the order of the indices doesn't matter:

contractES /@ {x[a, b, a, b] x[c, d, a, b], 
  x[m, n, a, b] x[c, d, m, n]}

Out[9]= {nsize^2 KroneckerDelta[a, b] KroneckerDelta[c, d], 
 nsize^2 KroneckerDelta[a, b] KroneckerDelta[c, d]}

And one last one with more x's:

contractES[x[a, b, m, n] x[m, n, c, d] x[c, d, g, h]]

Out[10]= nsize^4 KroneckerDelta[a, b] KroneckerDelta[g, h]

EDIT: Request for additional functionality.

As requested, here's code to (crudely) cover the case where you want differently dimensioned deltas. The basic idea is to build objects with constructors for both a group of products, all containing objects of the same size, and then call another constructor for each symbolic index, using TagSet to store the dimension common to this group. I also added code to allow for mixed integer/symbolic indices, and validity checking, of sorts.

You don't have to create all the "groups" at once. You do have to make all the objects in a given group the same size, since I pass in a single nMax parameter for each call. Assign the output from the constructor to either a symbol or an indexed symbol, as in the examples.

One important limitation: make sure your index names are unique. I didn't add any protection in the index constructor code, so if you repeat a symbol name, you'll overwrite the size[index] value stored, which will affect objects you create earlier. I suppose it would be possible to check to see if an index has already undergone processing, throw an error message and abort, but I'll leave that to you if you want it.

Another less critical limitation: I can't do validity checking on the indices when you either supply a symbolic size or a symbolic index. I..e if you define nMax = 4 and use 'a' for your index, or choose nMax = m and use '8' as an index, I have no way of knowing if it's in range or not, so that's for you to ensure.

Predicates to test for valid indices:

kdIndex::invindx = "Index `1` is invalid for size `2`";

validIndexQ[nMax_Integer][i_Integer] := IntervalMemberQ[Interval[{1,nMax}],i];

validIndexQ[_][_Symbol] = True;

validIndexQ[n_][i_] := (Message[kdIndex::invindx,i,n];False);

Index and group constructors:

makeKDIndex[nMax_][a_Symbol] /; validIndexQ[nMax][a] := size[a]^=nMax;

makeKroneckerIndexGroup[nMax_][texp:Times[expr__]] :=
        indices = Union@Flatten@Cases[expr, KroneckerDelta[a_,b_]:>{a,b},Infinity];

Transformation rules for contraction. I added two new rules to cover the case that we collapse down to a single delta raised to a power with no other terms, and hence no multiplication. The other "Times" based rules don't match for that case.

contractES[Times[x___,KroneckerDelta[a_,b_Symbol],y___,KroneckerDelta[c_,b_Symbol],z___]] := 
    size[b] contractES[Times[KroneckerDelta[a,c],x,y,z]];

contractES[Times[x___,KroneckerDelta[a_Integer,b_Symbol]^n_?(#>1&),y___]] := 
    size[b] contractES[Times[x,y,KroneckerDelta[a,b]^(n-2)]];

contractES[Times[x___,KroneckerDelta[a_Symbol,b_Symbol]^n_?(#>1&),y___]] := 
    size[a]size[b] contractES[Times[x,y,KroneckerDelta[a,b]^(n-2)]];

contractES[KroneckerDelta[a_Integer,b_Symbol]^n_?(#>1&)] :=
    size[b] contractES[KroneckerDelta[a,b]^(n-2)];

contractES[KroneckerDelta[a_Symbol,b_Symbol]^n_?(#>1&)] := 
    size[a] size[b] contractES[KroneckerDelta[a,b]^(n-2)];

contractES[x___] := x;

And some tests:

Construct some groups. Assign to some symbol or indexed symbol.

x[1] = makeKroneckerIndexGroup[4][
  KroneckerDelta[a, b] KroneckerDelta[c, b]]
Out[13]= KroneckerDelta[a, b] KroneckerDelta[b, c]

x[2] = makeKroneckerIndexGroup[5][KroneckerDelta[3,d]KroneckerDelta[e,3]KroneckerDelta[3,e]]
Out[14]= KroneckerDelta[3, d] KroneckerDelta[3, e]^2

x[3] = makeKroneckerIndexGroup[6][
  KroneckerDelta[q, r] KroneckerDelta[r, q]]
Out[16]= KroneckerDelta[q, r]^2

x[4] = makeKroneckerIndexGroup[mMax][
  KroneckerDelta[s, t] KroneckerDelta[u, s]]
Out[17]= KroneckerDelta[s, t] KroneckerDelta[s, u]

Check the sizes:

size /@ {a, b, c, d, e, q, r, s, t, u}
Out[23]= {4, 4, 4, 5, 5, 6, 6, mMax, mMax, mMax}

Some contractions:

contractES /@ {x[1] x[2], x[2] x[3], x[3], x[4], x[4] x[1] x[4]}
Out[24]= {20 KroneckerDelta[3, d] KroneckerDelta[a, c], 
 180 KroneckerDelta[3, d], 36, mMax KroneckerDelta[t, u], 
 4 mMax^4 KroneckerDelta[a, c]}

You can mix numbers (and other symbols) in with the deltas, as in the earlier answer.

contractES[x[1] 4.2 x[2] 5.03 x[3]]
Out[20]= 15210.7 KroneckerDelta[3, d] KroneckerDelta[a, c]

contractES[x[3] x[3] x[1] x[2] x[1] x[2]]
Out[21]= 41472000
  • $\begingroup$ Thank you, this has worked quite well so far. Would there be a way to implement this such that some of the indices have different ranges? For my end goal the different indices denote pairs like $a=(a_1,a_2)$ where $a_1=1,\ldots,N_1$ and $a_2=1,\ldots,N_2$. I tried adapting your sugestion by defining two different contractES for each, but haven't gotten it to work. $\endgroup$
    – dorrel
    Jul 13, 2023 at 12:25

Define a function and $N$ (recall that N is reserved/already defined)

nmax = 10;

Then we can write

\!\(\*UnderoverscriptBox[\(\[Sum]\), \(n = 1\), \(nmax\)]
  \(\*UnderoverscriptBox[\(\[Sum]\), \(m = 1\), \(nmax\)]
  \((x[a,b,m,n] x[m,n,c,d])\)\)\)

I copy-pasted that after typing escsumesc, but an equivalent input is

Sum[x[a,b,m,n] x[m,n,c,d], {m, nmax}, {n, nmax}]
  • $\begingroup$ Thank you for the suggestion. I tried it and iy works when, as you did I set nmax to a specific value. In my case I am inyerested in leaving it general, and in that case I haven't managed to get it, how would I need to change it? Also I assume this method would still work if $m$ and $n$ have different ranges nmax and mmax? $\endgroup$
    – dorrel
    Jul 13, 2023 at 15:02

Its standard job of replacement regardless of any algebra

  { Sum[x___  * KroneckerDelta[n_,p_]  KroneckerDelta[n_,q_] ,
d1___, {n_ , ___},  d2___ ]  :>  Sum[x  KroneckerDelta[p,q] ,d1,d2] , ...} 

plus same rule for all permutations of n,p,q in the arguments of the deltas


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.