Basic Methods in Mathematica
Ancillary file to 2nd Summer School on String theory and Holography, August 2021

Introduction

Why Mathematica?

Because Stephen Wolfram was a friend of Feynman, Gell-Mann, Anderson, Wilczek, ...

More suitable for symbolic / analytic computation in theoretical physics than other softwares, including
- integrate / series expand functions
- solving algebraic / differential equations
- guessing general formula from integer patterns

But there are troubles
- expensive
- math kernel is a blackbox, sometimes with bugs

How to get Mathematica?
- ask your University / Institute
- download trial version
- Raspberry Pi has it

Other options

Maple
popular for engineers; useful to double-check the numerical precision of Mathematica (GNU Multiple Precision Library)

SageMath, https://www.sagemath.org/index.html
free, but recommends paid plan for the cloud service (CoCalc)

Python
free, dynamic typing (or “Mathematica style”)
User-defined variables are initially all pointers, and Python determines the variable type during runtime.
for using GPUs, Deep learning, quantum computing, Google Colab, ...

References

Textbooks:
- Stephen Wolfram, “The Mathematica Book”
- Sal Mangano, “Mathematica Cookbook”

Courses, focusing on the application to HEP-TH
- Mathematica Summer School on Theoretical Physics, http://msstp.org/

Check Mathematica StackExchange

Notation, Structure and Tips

Notation

Use the proper notation

(* abbreviate built-in commands *)
FS=FullSimplify;
FF=FullForm;

(* Notation Package *)
<<Notation`

(* Human notation (LHS) = Machine notation (RHS) *)
{Notation[(NotationTemplateTag[(Subscript[δ, k_,l_])]) ⟺ (NotationTemplateTag[KroneckerDelta[k_,l_]])]};

KroneckerDelta[a,b]//FF
(KroneckerDeltaSeq[a,b])//FF

MMA interprets the two expressions in the same way

(* LHS is Subsuperscript, not Power *)
{Notation[(NotationTemplateTag[(Subsuperscript[(Superscript[X, a_]), i_, j_])]) ⟺ (NotationTemplateTag[X[a_,i_,j_]])]};

(Subsuperscript[(Superscript[Y, a_]), i_, j_])//FF

(A^2)//FF

Functional vs Procedural programming

Mathematica (Wolfram Language, or Lisp) is a functional programming, whereas C or Fortran is a procedural programming.
See https://community.wolfram.com/groups/-/m/t/1681993

RandomReal[1,10]

(* a small number of repetitions *)
ClearAll[repeat,sum]
repeat=10;
array=RandomReal[1,repeat];
sum=0;

FullForm[array]

AbsoluteTiming[For[i=1,i<=repeat,i++,sum+=array[[i]]];sum]
AbsoluteTiming[array/.List->Plus]
AbsoluteTiming[Total[array]]

The for loop is fast, while the list operation is slow

(* increase the number of repetitions *)
ClearAll[repeat,sum]
repeat=10^4;
array=RandomReal[1,repeat];
sum=0;

AbsoluteTiming[For[i=1,i<=repeat,i++,sum+=array[[i]]];sum]
AbsoluteTiming[array/.List->Plus]
AbsoluteTiming[Total[array]]

The for loop is slow, while the list operation is fast

Operators and Operations

Prefix, Postfix and Infix

(* Use "@" *)
Table[i,{i,3}]
{Total[%],Total@%}

(* Use "//" *)
Table[(Subscript[a, i]),{i,3}]
%//Total
%*Sum[(Subscript[b, i]),{i,2}]//Expand

(* use "~" *)
x~f~y
%/.f->Plus

(* Join text strings *)
a6=30;
"a"<>"6"//FF
ToExpression[%]

Collection of useful notations

Table[(Subscript[a, i]),{i,3}]
%//FF

(* "@@" rewrites the Head of the expressions *)
Apply[Times,Table[(Subscript[a, i]),{i,3}]]
Times@@Table[(Subscript[a, i]),{i,3}]

%//FF

(* "/@" distributes the operation *)
Map[Sqrt,Table[(Subscript[a, i]),{i,3}]]
Sqrt/@Table[(Subscript[a, i]),{i,3}]

(* Apply functions at higher levels of Table , or "@@@" *)
mat=Table[(Subscript[a, i])(Subscript[b, j])(Subscript[c, k]),{i,2},{j,3},{k,4}];
MatrixForm/@%

{Times@@mat,Times@@@mat,Apply[Times,mat,{2}]};
MatrixForm/@%

Multiply a’s, multiply b’s or multiply c’s, depending on the level “Times[]” is applied

Functions

f[x_]:=3+x;
f[y]

Function[3+#][x]
(3+#)&@x

Two variables

(3+#1+#2)&@@{x,y}

Abstract notation is useful when there are derivatives

eqϕ=-D[ϕ[x],{x,2}]-(m^2)ϕ[x]

eqϕ;
%/.ϕ[x_]:>Exp[I ω x]
%%/.ϕ:>(Exp[I ω #]&)//Factor

10 Tips for Writing Fast Mathematica Code

See https://blog.wolfram.com/2011/12/07/10-tips-for-writing-fast-mathematica-code/

1. Use floating-point numbers if you can, and use them early.  (★★★)
2. Learn about Compile…  (★)
3. Use built-in functions.   (★★)
4. Use Wolfram Workbench.  (★)
5. Remember values that you will need in the future.  (★★★★)

(* compare *)
fib[1]=1;
fib[2]=2;
fib[n_]:=fib[n]=fib[n-1]+fib[n-2];

fib2[1]=1;
fib2[2]=2;
fib2[n_]:=fib2[n-1]+fib2[n-2];

AbsoluteTiming[fib[25]]
AbsoluteTiming[fib2[25]]

Table[fib[n],{n,30}]

FindSequenceFunction[%]

6. Parallelize.  (★★★★)
6.5. Think about CUDALink and OpenCLLink.  (★)
7. Use Sow and Reap to accumulate large amounts of data (not AppendTo).  (★)
8. Use Block or With rather than Module.   (★★)
9. Go easy on pattern matching.  (★★)
10. Try doing things differently.  (★★★★)

LaunchKernels[]

Range[11,21]

Clear[array,array2];
array=Range[1000001,3000001];
array2=Range[3000001,5000001];
DistributeDefinitions[array,array2]

PrimeQ/@array//Tally//AbsoluteTiming
PrimeQ/@array2//Tally//AbsoluteTiming

ParallelMap[PrimeQ,array]//Tally//AbsoluteTiming
ParallelMap[PrimeQ,array2]//Tally//AbsoluteTiming

Parallelize is not efficient for simple tasks (like checking the prime number)

fl[n_]:=Length[FactorInteger[(10^n-1)/9]];
DistributeDefinitions[fl]

Map[fl,Range[55,65]]//AbsoluteTiming

ParallelMap[fl,Range[55,65]]//AbsoluteTiming

Parallelize is faster!

(Subscript[XXX, 1/2]) spin chains

Definitions

(* notation *)
FS=FullSimplify;
MF=MatrixForm;

(* commutation *)
comm[a_,b_]:=a.b-b.a;

(* 2x2 Pauli matrices *)
id[n_]:=DiagonalMatrix[ConstantArray[1,n]];
sp=Sqrt[2]{{0,1},{0,0}};
sm=Transpose[sp];
sx= {{0,1},{1,0}};
sy={{0,-I},{I,0}};
sz=DiagonalMatrix[{1,-1}];
cas=sp.sm+sm.sp+sz.sz;

Check notation

MF/@{sx,sy,sz,sp,sm}
FS[(#==0*IdentityMatrix[2])&/@{sx+I sy-Sqrt[2]sp,sx-I sy-Sqrt[2] sm}]

Check commutations

(#==0*IdentityMatrix[2])&/@(comm@@@{{sp,sm},{sz,sp},{sz,sm}}-{2sz,2sp,-2sm})

Construct an operator acting on the i-th and j-th site.  (total L sites)
(Subscript[L, ij]) acts trivially on the other sites.

(* substitute Pauli matrices at the sites i and j *)
generateTensors[i_,j_,LL_]:=Module[{iTab=Table[id[2],{n,LL}],spm,smp,szz},
spm=ReplacePart[iTab,{i->sp,j->sm}];
smp=ReplacePart[iTab,{i->sm,j->sp}];
szz=ReplacePart[iTab,{i->sz,j->sz}];
Array[(Subscript[a, #])&,{4}].(KroneckerProduct@@@{iTab,szz,spm,smp})
];

MF/@Table[id[2],{n,5}]
MF/@ReplacePart[%,3->sz]
Array[(Subscript[a, #])&,{4}]

Solving Yang-Baxter relations

Our ansatz for the R-matrix of XXX type


which should satisfy the Yang-Baxter relations

R[i_,j_]:=generateTensors[i,j,3]/.(Subscript[a, k_]):>(Subscript[b, k])((Subscript[u, i])-(Subscript[u, j]))+(Subscript[c, k])//FS;

R12=R[1,2];
R23=R[2,3];
R13=R[1,3];

generateTensors[1,2,2]//MF

R12//MF

R12.R13.R23-R23.R13.R12//FS//Flatten//Tally;
eqYB=%[[All,1]]

This equation should vanish for any u’s

Fold[CoefficientList,eqYB,{(Subscript[u, 1]),(Subscript[u, 2]),(Subscript[u, 3])}]//Flatten//FS//Tally;
eqYB2=%[[All,1]]

Array[u (Subscript[b, #])+(Subscript[c, #])&,{4}]/.Solve[eqYB2==0,Join[Array[(Subscript[b, #])&,{4}],Array[(Subscript[c, #])&,{4}]]];
solYB=%//DeleteDuplicates

Pick up the solution with (Subscript[a, 1])(Subscript[a, 2])(Subscript[a, 3])(Subscript[a, 4 ])!=0

Position[solYB,0];
solPos=Complement[Range[Length[solYB]],DeleteDuplicates[%[[All,1]]]];
Table[Table[((solYB[[a,k]])/(solYB[[a,4]])),{k,4}],{a,solPos}]//DeleteDuplicates//FS

Only two non-trivial solutions, and we can remove (b, c) by rescaling rapidity

(* Faddeev's solution, hep-th/9605187 *)
I generateTensors[1,2,2]/.{(Subscript[a, 1])->(1/2)-I u,(Subscript[a, 2])->(1/2),(Subscript[a, 3])->(1/2),(Subscript[a, 4])->(1/2)}//FS;
(* something else, perhaps related by change of basis *)
I generateTensors[1,2,2]/.{(Subscript[a, 1])->(1/2),(Subscript[a, 2])->(1/2)-I u,(Subscript[a, 3])->(1/2),(Subscript[a, 4])->(1/2)}//FS;
MF/@{%%,%}

Diagonalize Hamiltonian numerically

XXX spin chain Hamiltonian:

HXXX[L_]:=HXXX[L]=Sum[generateTensors[i,i+1,L],{i,L-1}]+generateTensors[L,1,L]/.(Subscript[a, 1])->1/.(Subscript[a, i_]):>-1;

HXXX[2]//MF
HXXX[2]//Eigenvalues//Tally

HXXX[3]//MF
HXXX[3]//Eigenvalues//Tally

HXXX[4]//MF
HXXX[4]//Eigenvalues
%//Tally

In general, the eigenvalues are highly irrational numbers

Table[Eigenvalues[HXXX[L]],{L,2,7}]//FS

Solution of the Bethe Ansatz

Bethe Ansatz Equations (k=1,2,...,M)


Periodic boundary conditions.
exp( i (Subscript[p, k])L) -> Sum[(Subscript[p, k]),{k,L}] = 0  (mod 2 π)

Regular :  all (Subscript[u, k]) are finite, distinct
Some u’s are infinite => su(2) descendants;  
XXX = su(2) of N=4 SYM,  in the full N=4 SYM the global symmetry is psu(2,2|4)

Energy of M-magnon states

eqBAE[L_,M_]:=Table[(((((u[k]+I/2)/(u[k]-I/2))))^L)+Product[((u[k]-u[j]+I)/(u[k]-u[j]-I)),{j,M}],{k,M}];
enBAE[M_]:=Sum[(8/(4((u[j])^2)+1)),{j,M}];

Take a closer look at the L=4 states, with M=1,2

Solve[eqBAE[4,1]==0,Table[u[i],{i,1}]]//FS
enBAE[1]/.%//FS

Consistent with the previous results modulo degeneracy. The regular solutions of BAE correspond to the su(2) HWS

Eigenvalues@HXXX[4]
Length@%

MM=2;
Solve[eqBAE[4,MM]==0,Table[u[i],{i,MM}]]//FS;
Table[u[i],{i,MM}]/.%
enBAE[MM]/.%%//FS

The solutions with irrational energies are singular (coincident roots) and unphysical

Clear[MM];
MM=3;
NSolve[eqBAE[4,MM]==0,Table[u[i],{i,MM}]]//Chop[#,(10^(-5))]&;
enBAE[MM]/.%//Chop[#,(10^(-5))]&

Not sensible solutions for M>L/2 (beyond the equators problem)

The classification of the physical solutions are discussed in 
Hao, Nepomechie, Sommese, 1308.4645
Marboe, Volin, 1608.06504
Complete solutions = (2^L) states (dim of irrep as (Subscript[S, L]))

Compute higher conserved charges

The next conserved charge is given by

which commute with the XXX Hamiltonian. Equivalently, (Subscript[Q, 3]) is generated by the boost operator


Further references:
Grabowski, Mathieu, hep-th/9412039
Beisert et al.  1308.1584,
Jiang et al. 1401.0384
de Leeuw et al. 1904.12005

BXXX[L_]:=BXXX[L]=Sum[i generateTensors[i,i+1,L],{i,L-1}]+L generateTensors[L,1,L]/.(Subscript[a, 1])->1/.(Subscript[a, i_]):>-1;

(* for periodic spin chain, we need to subtract the boundary term *)
Qbd[3,L_]:=Qbd[3,L]=comm[generateTensors[L,1,L],generateTensors[1,2,L]]/.(Subscript[a, 1])->1/.(Subscript[a, i_]):>-1;
QXXX[3,L_]:=QXXX[3,L]=-I(comm[BXXX[L],HXXX[L]]-L Qbd[3,L]);

(* explicit result *)
Qsr[3,L_]:=Qsr[3,L]=I(Sum[comm[generateTensors[i,i+1,L],generateTensors[i+1,i+2,L]],{i,L-2}]+comm[generateTensors[L-1,L,L],generateTensors[L,1,L]]+comm[generateTensors[L,1,L],generateTensors[1,2,L]])/.(Subscript[a, 1])->1/.(Subscript[a, i_]):>-1;

Check the equivalence of two results, and commutation with the XXX Hamiltonian

Table[QXXX[3,L]==Qsr[3,L],{L,2,5}]
Table[comm[Qsr[3,L],HXXX[L]]==0*IdentityMatrix[(2^L)],{L,2,5}]

Evaluate (Subscript[Q, 3]) in the basis where (Subscript[H, XXX]) is diagonal

Do[
{eval[ll],evec[ll]}=FS@Eigensystem[HXXX[ll]];
QsrDiag[3,ll]=Inverse[Transpose@evec[ll]].Qsr[3,ll].Transpose[evec[ll]]//FS;
,{ll,2,5}]

(* higher charge not diagonal *)
MF/@FS[{Inverse[Transpose@evec[4]].HXXX[4].Transpose[evec[4]],QsrDiag[3,4]}]

Resolve the ambiguity in the degenerate eigenvectors

eval[4]
Position[eval[4],4]//Flatten

{oval[4],ovec[4]}=Eigensystem[QsrDiag[3,4][[5;;11,5;;11]]];
rot[4]=ArrayFlatten[{{IdentityMatrix[4],0,0},{0,ovec[4],0},{0,0,IdentityMatrix[5]}}];
Length/@{evec[4],ovec[4],rot[4]}

(* now Q3 is diagonal at L=4 *)
MF/@FS[{Inverse[Transpose@rot[4]].Inverse[Transpose@evec[4]].HXXX[4].Transpose[evec[4]].Transpose[rot[4]],Inverse[Transpose@rot[4]].Inverse[Transpose@evec[4]].Qsr[3,4].Transpose[evec[4]].Transpose[rot[4]]}]

We can repeat the same computation at L=5

(* higher charge not diagonal *)
MF@QsrDiag[3,5]

Remove the zero eigenvalues in the last 6x6 block

{oval[5],ovec[5]}=Eigensystem[QsrDiag[3,5][[;;-7,;;-7]]];
rot[5]=ArrayFlatten[{{ovec[5],0},{0,IdentityMatrix[6]}}];
Length/@{evec[5],ovec[5],rot[5]}

Now diagonal

DiagonalMatrixQ[Inverse[Transpose@rot[5]].DiagonalMatrix[eval[5]].Transpose@rot[5]]
DiagonalMatrixQ[Inverse[Transpose@rot[5]].QsrDiag[3,5].Transpose@rot[5]]

Eigenvalues@QsrDiag[3,5]//FS