Einstein summation convention for symbolic vector calculus

Question

I am trying to do some vector calculus in Mathematica in index notation form because it gives a clear result that can be compared to pen and paper calculations. Since there is no built in Einstein summation convention this makes life rather difficult. I think I'm doing things the right way but if there is a better way please let me know!

I wish to evaluate vector derivatives of the tensor denoted $G$. Using the first of the linked questions below I make the definitions

R /: D[R, R[a_], NonConstants -> {R}] := R[a]/R
R /: D[R[a_], R[b_], NonConstants -> {R}] := KroneckerDelta[a, b]
G = KroneckerDelta[a, b]/R + TensorProduct[R[a], R[b]]/R^3

And evaluate

D[G, R[c], NonConstants -> {R}] // FullSimplify

So far so good, but suppose I wish to evaluate the Laplacian.

D[G, R[c], R[d], NonConstants -> {R}] //. d -> c // TensorExpand

Now we get terms like R[c]^2 that should automatically evaluate to R^2. This is not a problem, just add the lines

R /: R[a_] R[a_] := R^2
R /: R[a_] KroneckerDelta[a_, b_] := R[b]

And this works. However there are also terms like $\delta_{a,c}\delta_{b,c}$ that I want to automatically contract to $\delta_{a,b}$. These is no way to set a tag for KroneckerDelta since it is protected. Similarly if I want to take curls using the Levi Civita tensor this is going to get very ugly.

Tl;dr I have defined symbolic derivatives of tensors according to the first linked question since this allows me to get a result that I can visually compare to analytics. However I would also like to be able to automatically sum over all repeated indices. Is there a good way to do this? Thanks in advance!

Cartesian tensor gradient

symbolic summation involving kronecker delta

Implementing Einstein's summation convention for a particular case

EDIT: adding Levi-Civita functionality

I whipped up the ability to use the Levi-Civita tensor symbolically and have it work correctly (as far as I have tested). First its essential that $\delta_{a,a} = d = 3$ instead of 1 then I add total antisymmetry by modifying this question on totally antisymmetric functions.

Unprotect[KroneckerDelta]; KroneckerDelta /: 
 KroneckerDelta[a_, a_] := 3;
leviC /: MakeBoxes[leviC[a_, b_, c_], fmt_] := 
 MakeBoxes[Subscript[\[CurlyEpsilon], a, b, c], fmt]
leviC[a__] := Signature[{a}] (leviC @@ Sort@{a}) /; ! OrderedQ[{a}];
leviC[a__] := 0 /; ! DuplicateFreeQ[{a}];

Finally the Levi-Civita replacement rules are defined by

leviCrules := {leviC[indi_, indj_, indk_] leviC[indl_, indm_, 
     indn_] :> 
   KroneckerDelta[indi, 
      indl] (KroneckerDelta[indj, indm] KroneckerDelta[indk, indn] - 
       KroneckerDelta[indj, indn] KroneckerDelta[indk, indm]) - 
    KroneckerDelta[indi, 
      indm] (KroneckerDelta[indj, indl] KroneckerDelta[indk, indn] - 
       KroneckerDelta[indj, indn] KroneckerDelta[indk, indl]) + 
    KroneckerDelta[indi, 
      indn] (KroneckerDelta[indj, indl] KroneckerDelta[indk, indm] - 
       KroneckerDelta[indj, indm] KroneckerDelta[indk, indl]), 
  leviC[indi_indj _, indk_]^2 -> 6, 
  leviC[a___, b_, c_, d___] A_[Q___, b_, q___] A_[P___, c_, p___] :> 0}

Where the replacement rules here were taken from the code for VEST.

As far as I can tell this gives all the correct behaviour required for symbolic vector calculus! :)

Edit 2: dubious dummy index replacement

I add the rule

   A_[Q___, c__, q___] B_[P___, c__, p___] :> 
     A[Q, dummy, q] B[P, dummy, p]

so that Mathematica can recognize repeated indices.However this is dangerous for the obvious reason that $v_a v_a w_b w_b = v_a w_a v_b w_b$ in this scheme but I have been unable to do rule replacement such that each time the pattern is matched the dummy index is replaced by dummyi+1.

Answer 1

First let me remark that you can unprotect symbols by doing Unprotect[symbol], and then you can define whatever you want to it. This is however not always advised.

The way I usually deal with index contractions is to define a list of replacement rules

rules = {δ[X_, X_] :> d, 
 δ[X_, Y_]^2 :> d, δ[X_, Y_] δ[Y_, Z_] :> δ[X,Z],
 δ[X_, Y_] A_[a___, Y_, b___] :> A[a, X, b]};

which can then be applied to expr as expr//.rules. Actually often times if you have ugly expressions Mathematica won't recognize patterns. To avoid that one can define a function that not only replaces repeatedly but expands every time. For example

ReplaceExpand[expr_, rule_] := FixedPoint[Expand[#] /. rule &, expr];

---

Edit

I forgot a few things. First for this to work δ[X, Y] must be defined as an Orderless symbol

SetAttributes[δ,Orderless];

Second, the question also needed the Levi Civita tensor. Here is an implementation that I use (for four dimensions). The rules I define use OrderlessPatternSequence in order to deal with all possible ordering of the indices, but at the same time keeping track of the right signature. Here is the idea: if you define a rule like

A[x:OrderlessPatternSequence[a_,b_]] B[y:OrderlessPatternSequence[a_,b_]] :> f[x,y]

this will match both A[a,b]B[a,b] and A[a,b]B[b,a] but x and y will contain the right ordering. And if you want to know the signature of the permutation of a,b with respect to, say, x you can call Signature[PermutationList[FindPermutation[{x},{a,b}]]]. I'll write down examples. If four indices are contracted.

ε4[x:OrderlessPatternSequence[μ_,ν_,ρ_,λ_]]ε4[μ_,ν_,ρ_,λ_] :> -24*Signature[PermutationList[FindPermutation[{x},{μ,ν,ρ,λ}]]]

If three

ε4[x:OrderlessPatternSequence[μ_,ν_,ρ_,λ_]]ε4[y:OrderlessPatternSequence[κ_,ν_,ρ_,λ_]] :> -6*δ[μ,κ]*Signature[PermutationList[FindPermutation[{x},{μ,ν,ρ,λ}]]]*Signature[PermutationList[FindPermutation[{y},{κ,ν,ρ,λ}]]]

If two

ε4[x:OrderlessPatternSequence[μ_,ν_,ρ_,λ_]] ε4[y:OrderlessPatternSequence[κ_,τ_,ρ_,λ_]] :> -2*Det[Table[δ[ii,jj],{ii,{μ,ν}},{jj,{κ,τ}}]]*Signature[PermutationList[FindPermutation[{x},{μ,ν,ρ,λ}]]]*Signature[PermutationList[FindPermutation[{y},{κ,τ,ρ,λ}]]]

Hopefully you get the picture, if less indices are contracted you need a similar expression for the determinant of δ.

---

In my opinion, neither this is the best way possible. I don't know if this is going to be useful for you, but there is a very convenient formalism that gets rid of indices. Let me explain how it goes: In order to describe a bunch of vectors $x_1^{a_1}\cdots x_n^{a_n}$ with free indices $a_1\ldots a_n$, introduce auxiliary polarizations $w_1^{a_1}\cdots w_n^{a_n}$ that contracts those indices. You will end up with expressions involving only scalar products $w_i\cdot w_j$, $x_i\cdot w_j$ and $x_i\cdot x_j$. It's easy to recover the initial expressions since $$ w_i\cdot w_j\to\delta^{a_ia_j}\,,\qquad x_i\cdot w_j\to x_i^{a_j}\,,\qquad x_i\cdot x_j\to x_i^b x_j^b\,. $$ This is also convenient in that you don't have to come up with a new name of indices when they are summed over. If you have several $x$'s that are equal and you want to enforce the fact that their indices are symmetrized, you can contract all of them with the same $w$. So how to do it concretely?

Define the constants

You can define a function ConstQ that tells Mathematica what is a constant and what is a vector, in order to define linearity over them. These definitions should do the trick.

ConstQ[x_?NumberQ] := True;
ConstQ[x_Plus] := And @@ (ConstQ[#] & /@ x);
ConstQ[x_Times] := And @@ (ConstQ[#] & /@ x);

To which one can add any number of custom symbols with TagSet, e.g.

A /: ConstQ[A] := True;

will take A out of every linear operator. Similarly define what is a scalar, i.e. a fully contracted expression

ScalarQ[x_Power] := True;
ScalarQ[x_CenterDot] := True;
ScalarQ[x_?ConstQ] := True;

As you see I used the function CenterDot. This is the function that we are going to use for defining our scalar product. You can also come up with any name, it's just that CenterDot has a nicer built in appearance.

Define the scalar product

We need to define its linearity properties and the fact that it's symmetric:

SetAttributes[CenterDot, Orderless];
CenterDot[a_,b_Plus] := Plus @@ (CenterDot[a,#] & /@ b);
CenterDot[a_,c_*d_] /; ScalarQ[c] := c*CenterDot[a,d];
CenterDot[0,x_] := 0;

CenterDot has no built in meaning other than the fact that it looks like $a\cdot b$ and it can be entered as esc.esc.

Define the derivatives

I will not be using D as a derivative. Rather I will define an operator d[w,x,expr] that performs $$ w^a\frac{\partial}{\partial x^a}(\mathrm{expr})\,, $$ in the spirit of leaving no free indices around. This set of rules defines linearity

d[w_, x_, expr_Plus] := Plus @@ (d[w, x, #] & /@ List @@ expr);
d[w_, x_, y_?ConstQ] := 0;

See? ConstQ is useful now. This set of rules is the Leibniz rule for product and powers

d[w_, x_, expr_Times] :=   Plus @@ Table[MapAt[d[w, x, #] &, expr,i],{i, 1, Length[List @@ expr]}];
d[w_, x_, expr_^power_] := power expr^(power - 1) d[w, x, expr];

and this final set of rules actually defines the action on scalar products

d[w_, x_, CenterDot[x_,x_]] := 2 CenterDot[x,w];
d[w_, x_, CenterDot[x_,y_]] /; y =!= x := CenterDot[y,w];
d[w_, x_, CenterDot[y_,z_]] /; y =!= x && z =!= x := 0;

Define the Laplacian

This introduces a new thing: removing auxiliary polarizations. Let me explain: in order to do the Laplacian we can differentiate with two different new polarizations $$ w^b\frac{\partial}{\partial x^b}w^a\frac{\partial}{\partial x^a}(\mathrm{expr})\,, $$ and then replace $w^a w^b \to \delta^{ab}$. This function does exactly that

dSq[x_, expr_] := Expand[d[v, x, d[w, x, expr]]] //. 
 {CenterDot[v,w] -> dimensions, CenterDot[e_,v] CenterDot[f_,w] :> CenterDot[e,f]};

where dimensions is the dimension of the space.

As a bonus you can very easily check that $|x-y|^{(2-D)/2}$ is harmonic in $D$ dimensions.

Simplify[dSq[x, (CenterDot[x,x])^((2 - dimensions)/2)]]

0