1

Guide.nb

A Guide to Mathematica Packages
for Physicists
Yi-Zhuang You, UCSD
2019-10-05

Abstract
This document provides a guide to the bundle of Mathematica packages for physicists.
The bundle is available in the Github repository.
Keywords: Mathematica

Introduction◼ Content
This repository contains Mathematica packages and stylesheets that are useful for me.
 Package PauliAlgebra: symbolic handling the algebra and representation of Pauli operators
 LoopIntegrate: performing loop integration in quantum ﬁeld theory (with dimension

regularization)
 MatsubaraSum: performing Matsubara summation analytically DiagramEditor: an interactive editor of Feynman diagrams (no diagrammatic evaluation) Themes: a self-made plot theme for Mathematica, called “Academic”
 Toolkit: miscellaneous functions, including BZPlot for plotting band structure, tTr for tensor
network contraction, ComplexMatrixPlot for complex matrix visualization, Pf for matrix

Pfaffian Stylesheet
 CMU Article: Mathematica style sheet based on Computer Modern Unicode fonts (thefonts need to be installed separately to the operating system) FrontEnd Conﬁguration

2

Guide.nb◼ Installation Instruction
The bundle of packages can be downloaded fromhttps://github.com/EverettYou/Mathematica-for-physics
To install everything:
1. unzip this repository in a folder,
2. open install.m in Mathematica,
3. click the Run Package button to the top right,
4. quit Mathematica and restart.

PauliAlgebra Package◼ Overview
The package PauliAlgebra provides the following functions:? PauliAlgebra`*
PauliAlgebra`
Abstract

LieAlgebraσExp

ActionSpacenTr

σHermitian

Anticommutator

OrthogonalTransform

σInverse

AnticommuteQ

Qubit

σLog

AntisymmetricQ

Represent

σPolynomialQ

C4

Swap

σPowerCl

SymmetricQ

σSelect

Commutator

UnitaryTransform

σSqrt

CommuteQσσTr

ConjugateTransform

σ0

σTranspose

Controlled

σConjugate

Hadamard

σDet

3

Guide.nb◼ Arithmetic◼ Symbolic Representation
Four basic Pauli matrices are deﬁned as
σ0 = 

1 0
0 1
0 -ⅈ
1 0
 , σ1 = 
 , σ2 =  , σ3 = 
.
0 1
1 0
0 -1
ⅈ 0

( 1)

Pauli operators are tensor products of series of Pauli matrices,
σabc… = σa ⊗ σb ⊗ σc ⊗ …,

( 2)

which can be input as σ[a,b,c,…] directly. The indices a, b, c, … only take values of 0, 1, 2, 3.
◻ Tensor ProductPauli operators can explicitly constructed by the tensor product ( ⊗ is entered as c*).
σ[a] ⊗ σ[b] ⊗ σ[c]σ[a, b, c]

The tensor product construction will be translated by Mathematica to the short-handed notation
automatically.
◻ Dot ProductThe composition of Pauli operators is denoted by a dot product ·, which can be entered as
.. Dot product is carried out according to the Pauli algebra. For example,
σ[2] · σ[3]ⅈ σ[1]
σ[1] · σ[1]σ[0]Following shows the multiplication table calculated by Mathematica.
TableFormTableσi · σj, i, 0, 3, j, 0, 3,TableHeadings → Tableσi, i, 0, 3, Tableσj, j, 0, 3
σ[0]
σ[1]
σ[2]
σ[3]

σ[0]
σ[0]
σ[1]
σ[2]
σ[3]

σ[1]
σ[1]
σ[0]
- ⅈ σ[3]ⅈ σ[2]σ[2]σ[2]ⅈ σ[3]σ[0]- ⅈ σ[1]σ[3]σ[3]- ⅈ σ[2]ⅈ σ[1]
σ[0]◼ Matrix Representation
The matrix representation of a Pauli operator can be constructed by Represent in the formof a sparse array.

4

Guide.nb

σ[1] // Represent�������� ��������� ������������ {�� �}

SparseArray

σ[1, 2] // Represent // MatrixForm
0 0
0 0
0 -ⅈ
ⅈ 0

0 -ⅈ
ⅈ 0
0 0
0 0

Symbolic expression can also represented
Represent[a σ[1, 1] + b σ[3, 3]]�������� ��������� ������������ {�� �}

SparseArray

Use MatrixForm to convert the sparse array to its dense version for display.
% // MatrixFormb 0 0
0 -b a0 a -b
a 0 0

a
0
0
b

Use Abstract to recover the Pauli operator from the above matrix.
Abstract[%]a σ[1, 1] + b σ[3, 3]

Abstraction is just the inverse function of representation. Here is one more example.
AbstractDiagonalMatrix[{0, 0, 0, 1}]14

σ[0, 0] -1
4

σ[0, 3] -1
4

σ[3, 0] +1
4

σ[3, 3]◼ Dimension
Qubit gives the qubit number, i.e. the log 2 of the matrix dimension.
Qubit[σ[1, 2, 3]]3This means that the representation space of σ123 is spanned by 3 qubits or 23 = 8 states.
The identity operator of a given qubit number n can be constructed by σ0[n],5Guide.nb

σ0[3]σ[0, 0, 0]

or the qubit number can be automatically inferred from an operator
σ0[σ[1, 0, 0] + σ[2, 0, 0]]σ[0, 0, 0]◼ Distributive Properties
The tensor product and dot product automatically distributes over plus, such that the Paulialgebra expression is always expanded.
(σ[1] + 2 σ[2]) ⊗ (σ[1] - σ[2])σ[1, 1] - σ[1, 2] + 2 σ[2, 1] - 2 σ[2, 2](σ[1] + 2 σ[2]) · (σ[1] - σ[2])- σ[0] - 3 ⅈ σ[3]◼ Conjugation
Three types of conjugations are deﬁned:
 Complex conjugation A → A*σConjugate[σ[0] + σ[2] + ⅈ σ[3]]σ[0] - σ[2] - ⅈ σ[3] Transpose A → AσTranspose[σ[0] + σ[2] + ⅈ σ[3]]σ[0] - σ[2] + ⅈ σ[3] Hermitian conjugate A → A†
σHermitian[σ[0] + σ[2] + ⅈ σ[3]]σ[0] + σ[2] - ⅈ σ[3]◼ TraceσTr gives the trace of the Pauli operator.
σTr[σ[0, 0] + σ[3, 3]]2

6

Guide.nb◼ Operator Algebra◼ Action Space
A generic Hermitian operator can always be expanded as a superposition of Pauli operators
A = A[μ] σ[μ] ,
[μ]

( 3)

where [μ] = μ1 μ2 … labels the Pauli basis. The superposition coeﬃcient is given by
A[μ] =1
D

Tr A σ[μ] ,( 4)where D is the Hilbert space dimension. So each operator can be mapped to a super-state in theoperator space
A → A = A[μ] [μ]〉.
[μ]

( 5)

Composition of two operators can be calculated using operator product expansion,[μ][ν]
 c[λ] A[μ] B[ν] σ[λ] .

AB =[μ],[ν],[λ]( 6)Therefore each operator can also be represented as a super-operator in the operator spaceA→=[μ][ν] [λ]〉 c[λ] A[μ] 〈[ν],[μ],[ν],[λ]( 7)such that. The advantage of representing A in the operator space is to reduce
the representation dimension, because an operator acting on itself often only spans a subspace
whose dimension is smaller than the Hilbert space dimension.
Such subspace is called the action space of an operator. The matrix representation of an
operator in the action space and the corresponding basis can be found by ActionSpace.
MatrixForm /@ ActionSpace[a σ[1, 1, 0] + b σ[1, 0, 0] + c σ[0, 1, 0]]
0
c

b
a

c
0
a
b

b
a
0
c

a
b
,
c
0

σ[0, 0, 0]
σ[0, 1, 0]
σ[1, 0, 0]
σ[1, 1, 0]

Operator algebra can be carried out in the action space.

◼ Operator PowerσPower[A, n] returns the nth power of an operator A, i.e. An .
σPower[a σ[1, 1] + b σ[1, 0] + c σ[0, 1], 3]
6 a b c σ[0, 0] + 3 a2 c + 3 b2 c + c3  σ[0, 1] +
3 a2 b + b3 + 3 b c2  σ[1, 0] + a3 + 3 a b2 + 3 a c2  σ[1, 1]

7Guide.nb

σPower[a σ[0] + b σ[3], - 4]
4 a2 b2
- a2 + b2 4

a2 + b2 +2

- a2 + b2 

σ[0] -

4

4 a b a2 + b2  σ[3]- a2 + b2 

4

 The package uses divide-and-conquer algorithm for fast computation of high integer powers
Timing[σPower[a σ[1] + b σ[3], 99]]
2 16

2

0.003612, a3 a3 + a b2  + a2 b + b3  
22 16

a2 b a3 + a b2  + a2 b + b3  

2 16

2

+ a b2 a3 + a b2  + a2 b + b3    σ[1] +2

2 16+ b3 a3 + a b2  + a2 b + b3    σ[3] The inverse operator (A → A-1 ) can also be obtained via σInverse[A].
σInverse[a σ[0] + b σ[1] + c σ[2] + d σ[3]]- a σ[0] + b σ[1] + c σ[2] + d σ[3]- a2 + b 2 + c 2 + d 2

Singular matrix can not be inverted.
σInverse[σ[0] + σ[3]]������������ ������ �������� ����������� ���� ��� �� ���������
�������� ������ σ[�] + σ[�] �� ���������

σInverse[σ[0] + σ[3]] The square root operator (A → A1/2 ) can also be obtained via σSqrt[A].
σSqrt[a σ[0] + b σ[1]]
a-b
2

+

a+b
2

σ[0] + -

a-b
2

+

a+b
2

σ[1]

◼ Determinant
Determinant of an operator can be computed as
σDet[b σ[1, 3] + c σ[2, 2]]
b4 - 2 b2 c2 + c4

For large matrix, the algorithm is faster than the ordinary symbolic determinant algorithm, if thematrix has simple decomposition in terms of Pauli matrices.
σDet[b σ[0, 0, 0, 0, 0, 0, 0, 0, 0] + c σ[1, 3, 2, 2, 3, 1, 3, 2, 1]]b2 - c2 

256

8

Guide.nb◼ Operator Exp and Log
 Operator exponential (A → eA )σExp[ⅈ ϕ σ[3, 2]]
Cos[ϕ] σ[0, 0] + ⅈ Sin[ϕ] σ[3, 2]

 Operator logarithm (A → ln A)
σLog[(σ[0, 0] + ⅈ Sqrt[3] σ[1, 3]) / 2]13ⅈ π σ[1, 3]◼ Cliﬀord Algebra and Lie Algebra◼ Commutators
 Commutator A, B
Commutator[σ[1], σ[3]]- 2 ⅈ σ[2] Anticommutator A, B
Anticommutator[σ[3, 1], σ[2, 2]]
2 σ[1, 3]◼ Cliﬀord Algebra Cl[n] provides a choice of the generators of complex Cliﬀord algebra ℓn .
∀ i, j = 1, …, n : {γi , γj } = 2 δij .

( 8)

Cl[4]
{σ[1, 0], σ[2, 0], σ[3, 1], σ[3, 2]}

 Cl[p,q] provides a choice of the generators of real Cliﬀord algebra ℓp,q .
∀ i, j = 1, …, p + q : {γi , γj } = 2 δij ,∀ i = 1, …, p : γi  = γi ,∀ i = p + 1, …, p + q : γi  = -γi .
Cl[2, 3]{σ[1, 0, 0], σ[3, 1, 0], σ[2, 0, 0], σ[3, 2, 0], σ[3, 3, 2]}

( 9)

9

Guide.nb◼ Lie Algebra
LieAlgebra[{g1,g2,…}] completes the Lie algebra generators
LieAlgebra[{σ[1, 2], σ[3, 0]}]{σ[1, 2], σ[3, 0], σ[2, 2]}◼ Operator Transformations◼ Basis Gates
 C4[A] gives C4 rotation generated by a Pauli operator σ[μ]ⅈπ[μ]e4 σ ≡

1 + ⅈ σ[μ]
2

.

(10)

C4[σ[2, 3]]
σ[0, 0] + ⅈ σ[2, 3]
2

 Swap[σ[0,…,0,_,0,…,0,_,0,…,0]] gives the swap operator that exchange the two qubits

masked by _.
Swap[σ[0, _, 0, _]]12(σ[0, 0, 0, 0] + σ[0, 1, 0, 1] + σ[0, 2, 0, 2] + σ[0, 3, 0, 3])

 Hadamard[σ[0,…,0,_,0,…,0]] gives the Hadamard gate acting on the single qubit masked by_.
Hadamard[σ[0, 0, _, 0]]σ[0, 0, 1, 0] + σ[0, 0, 3, 0]
2 Controlled[σ[…,μ,_,ν,…]] gives the control gate, which implements σ[…μ0ν…] controlled by

the qubit masked by _.
Controlled[σ[_, 1, 2]]12(σ[0, 0, 0] + σ[0, 1, 2] + σ[3, 0, 0] - σ[3, 1, 2])◼ Transformations
Three types of transformations are deﬁned:
 OrthogonalTransform[O] represents the orthogonal transformation

10

Guide.nb

A → O A O.

(11)

 UnitaryTransform[O] represents the unitary transformation

A → O† A O.

(12)

 ConjugateTransform[O] represents the conjugate transformation

A → O-1 A O.
They can be implemented as
OrthogonalTransform[C4[σ[2]]][σ[3]]σ[1]or can be applied to a list of operators
UnitaryTransform[Controlled[σ[_, 1]]] /@ {σ[1, 0], σ[3, 0], σ[0, 1], σ[0, 3]}{σ[1, 1], σ[3, 0], σ[0, 1], σ[3, 3]}

It is convenient to use AssociationMap to view the transformation
AssociationMapUnitaryTransform[Hadamard[σ[_, 0]]],{σ[3, 1], σ[3, 2], σ[3, 3], σ[1, 0], σ[2, 0]}
σ[3, 1] → σ[1, 1], σ[3, 2] → σ[1, 2],
σ[3, 3] → σ[1, 3], σ[1, 0] → σ[3, 0], σ[2, 0] → - σ[2, 0]

◼ Pauli Operator Selection◼ Boolean Functions
The following Boolean function are useful in setting selection criterion.
 Commutation relation
CommuteQ[σ[1, 2], σ[3, 1]]
True
AnticommuteQ[σ[1], σ[3]]True Symmetry condition
SymmetricQ[σ[2]]
False

(13)

11

Guide.nb

AntisymmetricQ[σ[2]]True◼ Pauli SelectσSelect[{criterion,…},n] selects the Pauli operators that satisfy given criterion. Number of

qubits can be speciﬁed by n.
σSelectAnticommuteQ /@ {σ[0, 2], σ[1, 1]}{σ[0, 3], σ[1, 3], σ[2, 1], σ[3, 1]}

Without any criterion, all Pauli basis are returned
σSelect[{}, 2]{σ[0, 0], σ[0, 1], σ[0, 2], σ[0, 3], σ[1, 0], σ[1, 1], σ[1, 2], σ[1, 3],
σ[2, 0], σ[2, 1], σ[2, 2], σ[2, 3], σ[3, 0], σ[3, 1], σ[3, 2], σ[3, 3]}

LoopIntegrate Package◼ Overview
The package LoopIntegrate provides the following functions:? LoopIntegrate`*
LoopIntegrate`
DimensionRegularize

LeviCivitaEpsilon

MomentumIntegrate

FeynmanParameterize

Loop

MomentumShift

Index

LoopIntegrate

ParameterReduce

IntegrandInformation

LoopReduce◼ Loop Integral and Dimensional Regularization◼ Loop Object
The central object of the package is called Loop. It is a symbolic container of the data thatdeﬁnes a loop integral in general dimensions. Its structure is likeLoop[expr,{p1,…},D,x,nx]12Guide.nb

 expr - the integrand of the loop integral.
 {p1,…} - a list of the integral variables, specifying the momenta to be integrated over.
 D - the symbol for the spacetime dimension. Better not specify an integer dimension here, unlessthe integral does not need to be regularize. (One should use DimensionRegularize to properly

calculate loop integrals in speciﬁc dimensions).
 x - the symbol for Feynman parameter.
 nx - the number of Feynman parameter.
D, x and nx are optional. To create a Loop object, typically one just need to specify the integrandand variables.
The Loop object is represented as a loop integral.
Loop[1 / (p ^ 2 + m ^ 2), p]
ⅆ p


(2 π)

1


m2 + p2

Loop[Index[k + 2 p, μ] / ((k + p) ^ 2 q ^ 3), {p, q}]
ⅆ p


(2 π)

ⅆ q
 

(2 π)

kμ + 2 pμ(k + p)2 q3◼ Momentum Indexing
The momentum can be indexed by Index.
Index[k, μ]kμIndex distributes into the linear combination of momenta.
Index[(k + 2 p) / 3, μ]kμ + 2 pμ3For a concrete vector with a speciﬁc index, it extracts the component.
Index[{k1, k2, k3}, 2]k2◼ Feynman Parameterization
FeynmanParameterize analyzes the denominator of the integrand and combines the denominatorfactors. The Feynman parameters are introduced automatically.

13

Guide.nb# → FeynmanParameterize[#] & @ Loop[Index[k + q, μ] / ((k + q) ^ 2 q ^ 3), q]
ⅆ q


(2 π)

kμ + q μ→

q3 (k + q)2

32
 ⅆ x

2

△

ⅆ q
(2 π)

(kμ + qμ )
x1

q2 x1 + (k + q)2 x2 

5/2

The small △ of ∫△ ⅆn x reminds that the integral must be performed under the constraint ∑ni=1 xi = 1(which is a codimension-1 simplex △). Demonstrate the Feynman parameterization formula for
multiple factors of the denominator.
FeynmanParameterize[Loop[1 / ((p ^ 2 + Δ1) ^ a1 (p ^ 2 + Δ2) ^ a2 (p ^ 2 + Δ3) ^ a3), p]]Gamma[a1 + a2 + a3]ⅆ p3

 ⅆ x
△(2 π)



x1 -1+a1 x2 -1+a2 x3 -1+a3 p2 + Δ1 x1 + p2 + Δ2 x2 + p2 + Δ3 x3 

-a1-a2-a3

(Gamma[a1] Gamma[a2] Gamma[a3])◼ Loop Reduction
LoopReduce applies the momentum integral formulas to reduce the loop integral. The momenta

are integrated out one by one. For each momentum integral, a momentum shift is ﬁrst performed
to transform the denominator to the standard form. Then the integral formulas are applied. Delta
symbols δμν … are generated automatically under Wick contraction.
LoopReduce[Loop[1 / (p ^ 2 + Δ) ^ a, p]]2- π-/2 Δ

-a+


2

Gammaa -  
2Gamma[a]LoopReduce[Loop[Index[p, μ] × Index[p, ν] / (p ^ 2 + Δ) ^ a, p]]
2-1- π-/2 Δ

1-a+


2

Gamma- 1 + a -   δμ,ν2

Gamma[a]LoopReduce[Loop[Index[p, μ] × Index[p, ν] × Index[p, λ] × Index[p, κ] / (p ^ 2 + Δ) ^ a, p]]
2-2- π-/2 Δ

2-a+


2

Gamma- 2 + a -   (δκ,ν δλ,μ + δκ,μ δλ,ν + δκ,λ δμ,ν )2

Gamma[a]14

Guide.nb

LoopReduce[Loop[Index[p, μ] × Index[p, ν] ×
Index[p, λ] × Index[p, κ] × Index[p, ρ] × Index[p, τ] / (p ^ 2 + Δ) ^ a, p]]1
Gamma[a]2-3- π-/2 Δ

3-a+



Gamma- 3 + a - (δκ,τ δλ,ρ δμ,ν + δκ,ρ δλ,τ δμ,ν + δκ,τ δλ,ν δμ,ρ + δκ,ν δλ,τ δμ,ρ +2δκ,ρ δλ,ν δμ,τ + δκ,ν δλ,ρ δμ,τ + δκ,τ δλ,μ δν,ρ + δκ,μ δλ,τ δν,ρ + δκ,λ δμ,τ δν,ρ +2δκ,ρ δλ,μ δν,τ + δκ,μ δλ,ρ δν,τ + δκ,λ δμ,ρ δν,τ + δκ,ν δλ,μ δρ,τ + δκ,μ δλ,ν δρ,τ + δκ,λ δμ,ν δρ,τ )Repeated indices are assumed to be summed over, so δμμ = D. For example we can show


ⅆD p

pμ pμ

(2 π)D p2 + Δ

=
a

p2

ⅆD p

(2 π)D p2 + Δa

.

(14)

LoopReduce[Loop[Index[p, μ] × Index[p, μ] / (p ^ 2 + Δ) ^ a, p]]
2-1- π-/2  Δ

1-a+


2

Gamma- 1 + a -  
2Gamma[a]LoopReduce[Loop[p ^ 2 / (p ^ 2 + Δ) ^ a, p]]2- π-/2 Δ

1-a+


2

Gamma- 1 + a -  
22- π-/2 Δ-

Gamma[- 1 + a]

1-a+


2

Gammaa -  
2Gamma[a]FullSimplify[% - %%]0

In general, Feynman parameters will be generated
LoopReduce[Loop[1 / ((k + p) ^ 2 p), p]]
2

-

1



π 2 2 Gamma
- -

3
2

3


-

2

  ⅆ x

- k2 (- 1 + x2 ) x2  2 2
- +

2

△

x1

◼ Dimensional Regularization
Explicitly specifying the dimension in Loop may leads to divergent result.
LoopReduce[Loop[1 / (p ^ 2 + m ^ 2), p, 2]]ComplexInfinity

One should ﬁrst perform the loop integral with a symbolic dimension, and then use the dimensional regularization to ﬁnd the regular part of the integral.

15

Guide.nb

DimensionRegularize[LoopReduce[Loop[1 / (p ^ 2 + m ^ 2), p, D]], D → 2]
Log m Λ-2π◼ Parameter Reduction
Finally Feynman parameters can be integrated over by parameter reduction. For example,
ParameterReduce[Loop[Indexed[x, 1] ^ 2 × Indexed[x, 2], {}, , x, 2]]
1
12

It is equivalent to the following integral,
IntegrateIndexed[x, 1] ^ 2 × Indexed[x, 2] × DiracDeltaSumIndexedx, i, i, 2 - 1,x ∈ Rectangle[{0, 0}, {1, 1}]1
12

Apply to loop integral results.
ParameterReduceDimensionRegularize[LoopReduce[Loop[Index[k + p, μ] / ((k + p) ^ 2 p), p, D]], D → 3]
kμ Log k 
-

Λ

6 π2

◼ Function LoopIntegrate
The package also provide the high-level function LoopIntegrate that automates the above
procedures of Feynman parameterization, momentum shift and integration, dimensional regularization and parameter reduction.
LoopIntegrate[expr,p]LoopIntegrate[expr,p,D]
LoopIntegrate[expr,{p1,p2,…},D]It can be used with either general or speciﬁc dimensions.
FullSimplify @ LoopIntegrate[Index[k + p, μ] / ((k + p) ^ 2 p), p]
141- k2  2(-3+)

π-/2 Gamma 3-  Gamma[- 1 + ] kμ2Gamma- 1 + 2

If a speciﬁc dimension is given, dimensional regularization will be applied to obtain the regularpart of the integral.

16

Guide.nb

LoopIntegrate[Index[k + p, μ] / ((k + p) ^ 2 p), p, 3]
kμ Log k 
Λ

-

6 π2

MatsubaraSum Package◼ Overview
The package MatsubaraSum provides the following functions:? MatsubaraSum`*
MatsubaraSum`
Bosonic

Fermionic

ControlledPlane

MatsubaraSum

DistributionFunction

StatisticalSign

ZeroTemperatureLimit◼ Single Frequency Summation◼ Basic Summation
MatsubaraSum[f(z), z] evaluates the following summation

1β f (z ).
zwhere f (z ) is a function of the Matsubara frequency z = ⅈ ωn and ωn is taken from either one of thefollowing setsn ∈ ℤ : ωn =

2π
n
β
2π
n + 12 βBosonic,
Fermionic.

In the current version of the package, it is assumed that f (z ) is a fraction, whose denominator is apolynomial of z and whose numerator can be arbitrary.
Without specifying the type of the Matsubara frequency z, the general result is returned.
MatsubaraSum[1 / (z - ϵ), z]- nηz [ϵ] ηzηz = ± 1 is the statistical sign of the frequency z:

17Guide.nb

ηz = 

+ 1 if z ∈ Bosonic,- 1 if z ∈ Fermionic.

nη (ϵ) represents the distribution functionnη (ϵ) =1eβ ϵ - η= nB (ϵ) η = + 1,
nF (ϵ) η = - 1.

MatsubaraSum[1 / (z - ϵ) ^ 2, z]- ηz n′ηz [ϵ]MatsubaraSum[1 / (z - ϵ) ^ 3, z]-1
2

ηz n″ηz [ϵ]MatsubaraSum[1 / (z - ϵ) ^ 4, z]-1
6

ηz n(3)ηz [ϵ]n′η (ϵ), n″η (ϵ) and n(k)
η (ϵ) represent the 1st, 2nd and kth derivatives of the distribution function
respectively.

◼ Specify Frequency TypeOne can specify the Matsubara frequency type byz ∈ Bosonic : asserts z = ⅈ ωn to be bosonic,z ∈ Fermionic : asserts z = ⅈ ωn to be fermionic.

The symbol ∈ can be entered as el.
MatsubaraSum1 / (z - ϵ), z ∈ Fermionic
nF [ϵ]
MatsubaraSum1 / (z - ϵ), z ∈ Bosonic
- nB [ϵ]Another way to specify the frequency type is to use the Assumptions option.
MatsubaraSum1 / (z - ϵ), z, Assumptions → z ∈ Fermionic
nF [ϵ]
MatsubaraSum1 / (z - ϵ), z, Assumptions → z ∈ Bosonic
- nB [ϵ]The Assumptions option can also be use to specify the type of external frequencies (i.e. the Matsubara frequencies that will not be summed over).

18

Guide.nb

MatsubaraSum1 / (z1 - ϵ1) / (z1 - z2 - ϵ2), z1 ∈ Fermionic, Assumptions → z2 ∈ Bosonic
-nF [ϵ1]z2 - ϵ1 + ϵ2+

nF [ϵ2]z2 - ϵ1 + ϵ2

◼ Multiple Frequency Summation
Summation over multiple frequencies can be calculated by specifying more frequency variables
to be summed over.
MatsubaraSum[1 / (z1 - ϵ1) / (z2 - ϵ2), z1, z2]nηz1 [ϵ1] nηz2 [ϵ2] ηz1 ηz2Each frequency variable can be assigned a frequency type (either Bosonic or Fermionic).
MatsubaraSum1 / (z1 - ϵ1) / (z2 - ϵ2), z1 ∈ Bosonic, z2 ∈ Fermionic
- nF [ϵ2] nB [ϵ1]The external frequency type can be speciﬁed by Assumptions.
MatsubaraSum1 / (z1 - ϵ1) / (z2 - ϵ2) / (z1 + z2 - z - ϵ3),z1 ∈ Fermionic, z2 ∈ Fermionic, Assumptions → z ∈ Fermionic
-nF [ϵ2] (nF [ϵ1] - nF [ϵ3])z - ϵ1 - ϵ2 + ϵ3-

(nF [ϵ1] - nF [ϵ3]) nB [- ϵ1 + ϵ3]
z - ϵ1 - ϵ2 + ϵ3

◼ Miscellaneous◼ Simplify Distribution Functions
The package knowns about the properties of distribution functions, such asnη (- x) = -η - nη (x),
k (k)
n(k)
η (- x) = -(- 1) nη (x),2 n B (x ) nF (x ) = nB (x ) - nF (x ).
Using these, expressions of distribution functions can be simpliﬁed.
FullSimplify[1 - nF [- ϵ]]nF [ϵ]
FullSimplify[(1 + 2 nB [ϵ]) (1 - 2 nF [ϵ])]1Generic distribution functions and statistical signs can be reduced by making explicit assumptions about the types of frequency variables.

19

Guide.nb

Assumingz1 ∈ Bosonic && z2 ∈ Fermionic, nηz1 [ϵ1] nηz2 [ϵ2] ηz1 ηz2 - nF [ϵ2] nB [ϵ1]Assumingz ∈ Fermionic, nF [z + ϵ]
- nB [ϵ]◼ Zero Temperature Limit
The distribution function can further be expressed in terms of Heaviside Θ function in the zerotemperature limit.
Θ(ϵ) = 

1 ϵ > 0,0 ϵ < 0.

(15)

ZeroTemperatureLimit[nF [ϵ]]1 - HeavisideTheta[ϵ]ZeroTemperatureLimit[nB [ϵ]]- 1 + HeavisideTheta[ϵ]◼ Convergence Control
If the Matsubara summation does not converge, the result will depend on the regularization.
We regularize the summation by e-δ z factor with δ → 0,1βf (z ) e-δ z .
z

If δ = 0+ (δ = 0- ) the convergence is controlled on the right(left)-half plane. This choice can be set
by the option ControlledPlane, which can take Right (default) or Left.
TableFullSimplify @ MatsubaraSum1 / (z - ϵ), z ∈ Fermionic, ControlledPlane → cp,cp, Right, Left, AllnF [ϵ], - 1 + nF [ϵ], -1
2+ nF [ϵ]

For convergent summation, the result does not depend on regularization.
TableFullSimplify @ MatsubaraSum1 / (z - ϵ) ^ 2, z ∈ Fermionic, ControlledPlane → cp,cp, Right, Left, All{n′F [ϵ], n′F [ϵ], n′F [ϵ]}