PHYS 130C

Part 1: Quantum Many-Body Physics

Qubits, Bosons and Fermions

Introduction

What are Quantum Many-Body Systems?

Electrons in a solid, (CMP)

Ultra-cold atoms in an optical lattice, (AMO)

Elementary particles in a collider, (HEP)

Qubits in a quantum computer. (QI)

What Make Them Interesting?

Everything can be Emergent

Emergent Particle: Dirac fermions in graphene

Emergent Force: Gauge force in quantum magnets

Emergent Spacetime: ER = EPR

Emergent Gravity: SYK model

Conclusion: Despite the seemingly simple building blocks, many-body systems can display phenomena that are rich, surprising, and deeply beautiful.

Connection to Real-World Technology

High-temperature superconductor: strongly-correlated collective phenomenon of interacting electrons

Quantum computation: store and process information by quantum many-body entanglement

Mathematical and Conceptual Beauty

Topological phases of matter

Renormalization group (RG) theory

Quantum entanglement and tensor network

Emergent gauge theories

Holography and AdS/CFT correspondence

Qubits

Single Qubit

Quantum States

A pure state of a qubit is represented by a vector in a 2-dimensional Hilbert space, spanned by the basis states (Ket[0]) and (Ket[1]), which can be represented as one-hot vectors:

(Ket[0])≏(({{1}, {0}})), (Ket[1])≏(({{0}, {1}})).

A generic state of a qubit is a linear combination of the basis states

where (Subscript[ψ, 0]),(Subscript[ψ, 1])∈ℂ and ((|(Subscript[ψ, 0])|)^2)+((|(Subscript[ψ, 1])|)^2)=1.

Physical Observables

In quantum mechanics, physical observables are represented by Hermitian operators. The only observable for a single qubit is its spin or polarization, often measured along different axis.

The Pauli matrices are basis of Hermitian operators, corresponding to the spin observable along three orthogonal axes.

It is also useful to introduce the identity operator

such that any Hermitian observable O can be written as

O=c I+x X+y Y+z Z,

with c,x,y,z∈ℝ.

⚠️ Disclaimer: different notations of Pauli operators

Quantum information: X,Y,Z and I;

Condensed matter physics: σ=((σ^x),(σ^y),(σ^z))=((σ^1),(σ^2),(σ^3)) and (σ^0).

Measuring a qubit in the computational basis (Z eigen basis) means measuring its Z observable.

If the qubit is in the state (Ket[ψ])=(Subscript[ψ, 0])(Ket[0])+(Subscript[ψ, 1])(Ket[1]), measuring in the {(Ket[0]),(Ket[1])} basis will:

yield (Ket[0]) with probability ((|(Subscript[ψ, 0])|)^2),

yield (Ket[1]) with probability ((|(Subscript[ψ, 1])|)^2).

Two Qubits

Quantum States

Each qubit has two basis states: (Ket[0]) and (Ket[1]), so two qubits have four basis states, arranged as follows:

({{ , , qubit 2, }, {, , (Ket[0]), (Ket[1])}, {qubit 1, (Ket[0]), (Ket[00]), (Ket[01])}, {, (Ket[1]), (Ket[10]), (Ket[11])}})

We can choose to represent them as one-hot vectors:

({{0, 1, 2, 3}, {(Ket[00]), (Ket[01]), (Ket[10]), (Ket[11])}, {≏, ≏, ≏, ≏}, {(({{1}, {0}, {0}, {0}})), (({{0}, {1}, {0}, {0}})), (({{0}, {0}, {1}, {0}})), (({{0}, {0}, {0}, {1}}))}})

This is also called the computational basis ((Subscript[Z, 1]) and (Subscript[Z, 2]) eigen basis).

The mathematical construction of (Ket[01]) is a tensor product of (Subscript[(Ket[0]), 1st]) and (Subscript[(Ket[1]), 2nd]) states, denoted as

(Ket[01])=(Ket[0])⊗(Ket[1]).

How to make sense of ⊗ ?

Practically, tensor product of two vectors means

({{(Ket[0]), ⊗, (Ket[1]),  , =,  , (Ket[01])}, {≏,  , ≏,  ,  ,  , ≏}, {(({{1}, {0}})), ⊗, (({{0}, {1}})), =, (({{1(({{0}, {1}}))}, {0(({{0}, {1}}))}})), =, (({{0}, {1}, {0}, {0}}))}})

Formally, tensor product is a binary operation that is linear in both factors:

((Subscript[α, 0])(Ket[0])+(Subscript[α, 1])(Ket[1]))⊗((Subscript[β, 0])(Ket[0])+(Subscript[β, 1])(Ket[1]))
=(Subscript[α, 0])(Subscript[β, 0])(Ket[0])⊗(Ket[0])+(Subscript[α, 0])(Subscript[β, 1])(Ket[0])⊗(Ket[1])+(Subscript[α, 1])(Subscript[β, 0])(Ket[1])⊗(Ket[0])+(Subscript[α, 1])(Subscript[β, 1])(Ket[1])⊗(Ket[1])
(Overscript[=, short])(Subscript[α, 0])(Subscript[β, 0])(Ket[00])+(Subscript[α, 0])(Subscript[β, 1])(Ket[01])+(Subscript[α, 1])(Subscript[β, 0])(Ket[10])+(Subscript[α, 1])(Subscript[β, 1])(Ket[11]).

The only way to implement this algebra correctly on the vector level is to require the tensor product be computed as (see also Eq. ())

Rules of tensor product on basis vectors

(Ket[a])⊗(Ket[b])=(Ket[a b]),
(Bra[a])⊗(Bra[b])=(Bra[a b]),

((Ket[a])(Bra[b]))⊗((Ket[c])(Bra[d]))=(Ket[a c])(Bra[b d]),
(BraKet[a c, b d])=((Bra[a])⊗(Bra[c]))((Ket[b])⊗(Ket[d]))=(BraKet[a, b])(BraKet[c, d]).

Physical Observables

Introduce the Pauli operators (Subscript[X, 1]), (Subscript[Y, 1]), (Subscript[Z, 1]) for the first qubit, and (Subscript[X, 2]), (Subscript[Y, 2]), (Subscript[Z, 2]) for the second qubit, any physical observable of a two-qubit system is described by a Hermitian operator as a linear combination of the followings:

Trivial observable:

𝟙=I⊗I.

Single-qubit observables:

qubit 1:(Subscript[X, 1])=X⊗I,(Subscript[Y, 1])=Y⊗I,(Subscript[Z, 1])=Z⊗I,
qubit 2:(Subscript[X, 2])=I⊗X,(Subscript[Y, 2])=I⊗Y,(Subscript[Z, 2])=I⊗Z.

Two-qubit observables (joint measurements):

({{(Subscript[X, 1])(Subscript[X, 2]), (Subscript[X, 1])(Subscript[Y, 2]), (Subscript[X, 1])(Subscript[Z, 2])}, {(Subscript[Y, 1])(Subscript[X, 2]), (Subscript[Y, 1])(Subscript[Y, 2]), (Subscript[Y, 1])(Subscript[Z, 2])}, {(Subscript[Z, 1])(Subscript[X, 2]), (Subscript[Z, 1])(Subscript[Y, 2]), (Subscript[Z, 1])(Subscript[Z, 2])}})

For example, by (Subscript[Z, 1])(Subscript[Y, 2]), one means:

(Subscript[Z, 1])(Subscript[Y, 2])=(Z⊗I)(I⊗Y)=Z⊗Y.

Altogether there are 16 basis observables.

The tensor product of operator is defined as

A=(∑_(a b))(Ket[a])(Subscript[A, a b])(Bra[b]),B=(∑_(c d))(Ket[c])(Subscript[B, c d])(Bra[d]),
=>A⊗B=(∑_(a b c d))(Ket[a])⊗(Ket[c])(Subscript[A, a b])(Subscript[B, c d])(Bra[b])⊗(Bra[d]).

Practically, it is computed as

Z⊗Y≏(({{1, 0}, {0, -1}}))⊗(({{0, -i}, {i, 0}}))
=(({{1(({{0, -i}, {i, 0}})), 0(({{0, -i}, {i, 0}}))}, {0(({{0, -i}, {i, 0}})), -1(({{0, -i}, {i, 0}}))}}))=(({{0, -i, 0, 0}, {i, 0, 0, 0}, {0, 0, 0, i}, {0, 0, -i, 0}})).

In Mathematica, you can use KroneckerProduct to calculate tensor product of matrices, such as:

KroneckerProduct[PauliMatrix[3],PauliMatrix[2]]//MatrixForm

Represent the Hamiltonian H=(Subscript[X, 1])(Subscript[X, 2])+(Subscript[Y, 1])(Subscript[Y, 2])+(Subscript[Z, 1])(Subscript[Z, 2]) as a 4×4 matrix in the computational basis. Diagonalize the matrix to find the eigen energies and corresponding eigen states.

H=KroneckerProduct[PauliMatrix[1],PauliMatrix[1]]+KroneckerProduct[PauliMatrix[2],PauliMatrix[2]]+KroneckerProduct[PauliMatrix[3],PauliMatrix[3]];
H//MatrixForm

Eigensystem[H]

({{energy, state}, {-3, (1/(Sqrt[2]))((Ket[01])-(Ket[10]))}, {1, (Ket[00])}, {, (1/(Sqrt[2]))((Ket[01])+(Ket[10]))}, {, (Ket[11])}})

By definition, operator tensor product has the following properties:

(A⊗B)(C⊗D)=(A B)⊗(C D),

(α A+β B)⊗C=α A⊗C+β B⊗C,
A⊗(β B+γ C)=β A⊗B+γ A⊗C,

Tr(A⊗B)=(Tr A) (Tr B).

Many Qubits

ℋ=(Subscript[ℋ, 1])⊗(Subscript[ℋ, 2])⊗(Subscript[ℋ, 3])⊗…⊗(Subscript[ℋ, N]).

Quantum States: a generic many-qubit quantum state

(Ket[ψ])=(∑_((Subscript[a, 1])(Subscript[a, 2])…))(Subscript[ψ, (Subscript[a, 1])(Subscript[a, 2])…])(Ket[(Subscript[a, 1])(Subscript[a, 2])…])=(∑_([a]))(Subscript[ψ, [a]])(Ket[[a]]).

Notation: bundled index [a]=(Subscript[a, 1])(Subscript[a, 2])… (Subscript[a, N]), with (Subscript[a, i])=0,1, and

(Ket[[a]])=(Ket[(Subscript[a, 1])])⊗(Ket[(Subscript[a, 2])])⊗…⊗(Ket[(Subscript[a, N])]).

Physical Observables

A generic operator takes the form of

O=(∑_((Subscript[a, 1])(Subscript[a, 2])…))(∑_((Subscript[b, 1])(Subscript[b, 2])…))(Ket[(Subscript[a, 1])(Subscript[a, 2])…])(Subscript[O, (Subscript[a, 1])(Subscript[a, 2])…,(Subscript[b, 1])(Subscript[b, 2])…])(Bra[(Subscript[b, 1])(Subscript[b, 2])…])
=(∑_([a],[b]))(Ket[[a]])(Subscript[O, [a][b]])(Bra[[b]]).

Any operator O can be decomposed to a linear combination of Pauli operators P with some coefficients (Subscript[o, P]),

O=(∑_P)(Subscript[o, P])P,

where P takes the form of

P=(Subscript[P, 1])⊗(Subscript[P, 2])⊗… ⊗(Subscript[P, N]),

with (Subscript[P, i])∈{I,X,Y,Z}.

There are all together (4^N) Pauli operators (including identity) for a N-qubit system.

They form a set of complete and orthogonal operator basis.

All Pauli operators can be generated by (Subscript[X, i]) and (Subscript[Z, i]), constructed as

or defined by the following algebraic definitions:

for all i,j∈{1,2,…,N}.

What about (Subscript[Y, i])? - It is simply defined as (Subscript[Y, i]):=i (Subscript[X, i])(Subscript[Z, i]).

Such that any Pauli operator can be written as

(Subscript[P, ξ,ζ])=(i^(ξ·ζ))(∏_i)(Subscript[X, i]^(Subscript[ξ, i]))(∏_i)(Subscript[Z, i]^(Subscript[ζ, i])),

uniquely labeled by two binary vectors ξ=((Subscript[ξ, 1]),(Subscript[ξ, 2]),…,(Subscript[ξ, N])) and ζ=((Subscript[ζ, 1]),(Subscript[ζ, 2]),…,(Subscript[ζ, N])) with (Subscript[ξ, i]),(Subscript[ζ, i])=0,1. The factors of i are taken care by (i^(ξ·ζ)), where ξ·ζ means the dot product of the binary vectors.

Check that the following four operators all commute: 
(Subscript[X, 1])(Subscript[X, 2])(Subscript[X, 4]), (Subscript[Z, 1])(Subscript[Z, 2])(Subscript[X, 3])(Subscript[X, 4]), (Subscript[Y, 1])(Subscript[X, 2])(Subscript[Y, 3])(Subscript[Z, 4]), (Subscript[X, 1])(Subscript[Z, 2])(Subscript[Z, 4]),
using their operator algebra, without relying on explicit matrix representations.

({{(Subscript[X, 1]), (Subscript[X, 2]),  , (Subscript[X, 4])}, {(Subscript[Z, 1]), (Subscript[Z, 2]), (Subscript[X, 3]), (Subscript[X, 4])}})
=
({{(-), (-),  ,  }, {(Subscript[Z, 1]), (Subscript[Z, 2]), (Subscript[X, 3]), (Subscript[X, 4])}, {(Subscript[X, 1]), (Subscript[X, 2]),  , (Subscript[X, 4])}})

First align them by the qubit (each column corresponds to a qubit)

Then interchange the operators by swapping the two rows

Each time a pair of X and Z exchanges in the same column, a (-) sign is generated for that column

Count the total number of (-) signs. If even, then the two operator commutes.

Bosons

Single-Mode Bosons

Quantum states: Boson occupation states

({{(Ket[0]),  vacuum (empty)}, {(Ket[1]), one boson occupied}, {(Ket[2]), two bosons occupied}, {⋮, ⋮}, {(Ket[n]), n bosons occupied}, {⋮, ⋮}})

They are distinct eigenstates of the boson number operator (OverHat[n]),

where n∈ℕ is the corresponding eigenvalue.

They form a complete set of orthonormal basis,

They can be represented as one-hot vectors:

(Ket[0])≏(({{1}, {0}, {0}, {⋮}})), (Ket[1])≏(({{0}, {1}, {0}, {⋮}})), (Ket[2])≏(({{0}, {0}, {1}, {⋮}})), …

Comment: These vectors are infinite dimensional, because the Hilbert space dimension is infinite. In practice, we can (artificially) choose to truncate the Hilbert space dimension to some finite maximal boson number (Subscript[n, max]) (if we don’t care about states with that many bosons, as they are too high in energy).

Any quantum state of single-mode bosons must be a linear superposition of the boson occupation states

(Ket[ψ])=(∑_(n∈ℕ))(Subscript[ψ, n])(Ket[n]).

Boson Operators

Boson creation operator: (SuperDagger[b]) (Overscript[↔, †]) Boson annihilation operator: b
They are Hermitian conjugate to each other.

Actions on boson occupation states

Interpretation: creation operator raise the boson number by one, and annihilation operator lowers the boson number by one, BUT there are square root factors (also known as boson enhancement factors) in the front.

Boson number operator can be expressed as

Using Eq. () to show that Eq. () can reproduce Eq. () consistently.

(SuperDagger[b])b(Ket[n])=(SuperDagger[b])(Sqrt[n])(Ket[n-1])
=(Sqrt[n])(SuperDagger[b])(Ket[n-1])
=(Sqrt[n])(Sqrt[n])(Ket[n])
=n(Ket[n]).

Boson commutation relation

Note: commutator [A,B]:=A B-B A.

In fact, this should be treated as the algebraic definition of the boson creation and annihilation operators.

Use Eq. () to check that Eq. () holds for any state as operator identities.

[b,(SuperDagger[b])](Ket[n])=(b (SuperDagger[b])-(SuperDagger[b])b)(Ket[n])
=b (SuperDagger[b])(Ket[n])-(SuperDagger[b])b(Ket[n])
=b (Sqrt[n+1])(Ket[n+1])-(SuperDagger[b])(Sqrt[n])(Ket[n-1])
=(Sqrt[n+1])(Sqrt[n+1])(Ket[n])-(Sqrt[n])(Sqrt[n])(Ket[n])
=(n+1-n)(Ket[n])
=(Ket[n])
=𝟙(Ket[n]).

Boson operators can be represented as matrices

(SuperDagger[b])≏(({{0,  ,  ,  }, {(Sqrt[1]), 0,  ,  }, { , (Sqrt[2]), 0,  }, { ,  , ⋱, ⋱}})), b≏(({{0, (Sqrt[1]),  ,  }, { , 0, (Sqrt[2]),  }, { ,  , 0, ⋱}, { ,  ,  , ⋱}})),

such that the number operator is diagonal as expected

(SuperDagger[b])b≏(({{0,  ,  ,  }, { , 1,  ,  }, { ,  , 2,  }, { ,  ,  , ⋱}})).

Comment: These are infinite dimensional matrices, due to infinite dimensional Hilbert space. In practice, we can truncate them to finite dimension by setting a (artificial) maximal boson number (Subscript[n, max]).

Multi-Mode Bosons

a momentum eigenstate (a plane wave) for bosons to follow,

a position eigenstate (a site on an optical lattice) for atoms to stay,

a cavity eigenstate (a wave mode) for photons to populate,

a polarization mode (transverse/longitudinal) for phonon to take ...

Quantum States: Tensor product of boson occupation states

where (Subscript[n, i])∈ℕ is the number of bosons occupying the ith mode. 
[n]=(Subscript[n, 1])(Subscript[n, 2])(Subscript[n, 3])… denotes the sequence of boson occupation numbers.

They are joint eigenstates of boson number operators (Subscript[(OverHat[n]), i])

They form a complete set of orthonormal basis,

Any quantum many-body state of multi-mode bosons must be a linear superposition of these boson occupation states.

Boson Operators

For each mode i, define
Boson creation operator: (Subscript[b, i]^†) (Overscript[↔, †]) Boson annihilation operator: (Subscript[b, i])

Actions on boson occupation states

Comment: the boson creation/annihilation operator only acts on the tensor product Hilbert space of its own mode, and has no side effect on other modes.

Boson number operator of mode i can be express as

Boson commutation relations

These relations should be considered as the algebraic definition of the boson creation and annihilation operator for multiple modes.

Fermions

Single-Mode Fermions

Quantum states: Fermion occupation states

({{(Ket[0]),  vacuum (empty)}, {(Ket[1]), one fermion occupied}})

and no more state. The Hilbert space dimension is 2.

They are still distinct eigenstates of the fermion number operator (OverHat[n]),

Just now the eigenvalue can only be n=0,1.

They form a complete set of orthonormal basis, and can be represented as one-hot vectors:

(Ket[0])≏(({{1}, {0}})),(Ket[1])≏(({{0}, {1}})).

Question: This looks just the same as a single qubit, right?

Answer: Yes, indeed, a single-mode fermion and a single qubit are equivalent quantum systems. But multi-mode fermions and multi-qubit system will be different, as we will see.

Fermion Operators

Fermion creation operator: (SuperDagger[c]) (Overscript[↔, †]) Fermion annihilation operator: c
They are Hermitian conjugate to each other.

Actions on fermion occupation states

(SuperDagger[c])(Ket[0])=(Ket[1]), (SuperDagger[c])(Ket[1])=0;
c(Ket[0])=0, c(Ket[1])=(Ket[0]).

Note: be careful to distinguish the vacuum state (Ket[0]) (meaning that there is such a state, and the state contains no particle) and the zero vector 0 (meaning that there is no such a state). In quantum mechanics, a 0 probability amplitude indicates the state is impossible.

Empty state can not be empty: it is impossible to remove a fermion from the vacuum state, because there was no fermion there. Therefore, we have

c(Ket[0])=0.

Filled state can not be filled: it is impossible to add a second fermion to a mode when it is already occupied by one fermion — the Pauli exclusion principle. Therefore, we have

(SuperDagger[c])(Ket[1])=0.

Interestingly, Eq. () can be summarized in two lines

Interpretation: n↔1-n toggles between n=0 and n=1 state, as expected. The peculiar

Fermion number operator can be expressed as

Use Eq. () to show that Eq. () can reproduce Eq. () consistently.

(SuperDagger[c])c(Ket[n])=(SuperDagger[c])(Sqrt[n])(Ket[1-n])
=(Sqrt[n])(SuperDagger[c])(Ket[1-n])
=(Sqrt[n])(Sqrt[1-(1-n)])(Ket[1-(1-n)])
=n(Ket[n]).

Fermion anticommutation relation

Note: anticommutator {A,B}:=A B+B A.

In fact, this should be treated as the algebraic definition of the fermion creation and annihilation operators.

Use Eq. () to check that Eq. () holds for any state as operator identities.

{c,c}(Ket[n])=c c+c c(Ket[n])
=2 c c (Ket[n])
=2c (Sqrt[n])(Ket[1-n])
=2(Sqrt[n (1-n)])(Ket[n])

{c,(SuperDagger[c])}(Ket[n])=(c (SuperDagger[c])+(SuperDagger[c])c)(Ket[n])
=c (SuperDagger[c])(Ket[n])+(SuperDagger[c])c(Ket[n])
=c(Sqrt[1-n])(Ket[1-n])+(SuperDagger[c])(Sqrt[n])(Ket[1-n])
=(Sqrt[1-n])(Sqrt[1-n])(Ket[n])+(Sqrt[n])(Sqrt[n])(Ket[n])
=(1-n+n)(Ket[n])
=(Ket[n])
=𝟙(Ket[n]).

For a single mode, fermion operators have simple 2×2 matrix representations in the {(Ket[0]), |1⟩} basis

(SuperDagger[c])≏(({{0, 0}, {1, 0}})), c≏(({{0, 1}, {0, 0}})),

such that the number operator is diagonal as expected

(SuperDagger[c])c≏(({{0, 0}, {0, 1}})).

Jordan-Wigner transformation: If the fermion occupation states are mapped to the qubit states, fermion operators can also be represented as Pauli operators, as

(SuperDagger[c])=((X-i Y)/2), c=((X+i Y)/2),(SuperDagger[c])c=((I-Z)/2).

Multi-Mode Fermions

Quantum States: Tensor product of fermion occupation states

where (Subscript[n, i])∈ℕ is the number of fermions occupying the ith mode. 
[n]=(Subscript[n, 1])(Subscript[n, 2])(Subscript[n, 3])… denotes the sequence of fermion occupation numbers.

They are joint eigenstates of fermion number operators (Subscript[(OverHat[n]), i])

They form a complete set of orthonormal basis,

Any quantum many-body state of multi-mode fermions must be a linear superposition of these fermion occupation states.

Fermion Operators

For each mode i, define
Fermion creation operator: (Subscript[c, i]^†) (Overscript[↔, †]) Fermion annihilation operator: (Subscript[c, i])

Actions on fermion occupation states

Comment: the fermion creation/annihilation operator not only acts on its own mode. There is an important fermion sign that depends on the parity of the number of fermions occupying the proceeding modes (for this to make sense, a canonical ordering of single-particle mode is required), because each exchange of two fermions will change the sign of the whole many-body state — a defining feature of fermions.

Fermion number operator of mode i can be express as

Fermion commutation relations

These relations should be considered as the algebraic definition of the boson creation and annihilation operator for multiple modes.

Eq. () implies the following rules:

No double action: double creation or annihilation of fermions on the same mode i is forbidden

(Subscript[c, i])(Subscript[c, i])=0,
(Subscript[c, i]^†)(Subscript[c, i]^†)=0.

Exchange minus sign: exchanging two fermion operations of different modes (i≠j) will generate a minus sign (the fermion sign)

(Subscript[c, i])(Subscript[c, j])=-(Subscript[c, j])(Subscript[c, i]),
(Subscript[c, i]^†)(Subscript[c, j]^†)=-(Subscript[c, j]^†)(Subscript[c, i]^†),
(Subscript[c, i]^†)(Subscript[c, j])=-(Subscript[c, j])(Subscript[c, i]^†),
(Subscript[c, i])(Subscript[c, j]^†)=-(Subscript[c, j]^†)(Subscript[c, i]).

Matter and anti-matter: (Subscript[c, i]^†)(Subscript[c, i]) counts the number of particles in the mode i, whereas (Subscript[c, i])(Subscript[c, i]^†) counts the number of anti-particles in the same mode:

(Subscript[c, i])(Subscript[c, i]^†)=𝟙-(Subscript[c, i]^†)(Subscript[c, i]).

({{state, particle, anti-particle, (Subscript[c, i]^†)(Subscript[c, i]), (Subscript[c, i])(Subscript[c, i]^†)}, {(Subscript[(Ket[0]), i]), empty, occupied, 0, 1}, {(Subscript[(Ket[1]), i]), occupied, empty, 1, 0}})

Slogan: where no particle is found, an anti-particle is bound.

({{operator, particle, anti-particle}, {creation, (Subscript[c, i]^†), (Subscript[c, i])}, {annihilation, (Subscript[c, i]), (Subscript[c, i]^†)}, {number, (Subscript[c, i]^†)(Subscript[c, i]), (Subscript[c, i])(Subscript[c, i]^†)}})

It turns out that all matter in our universe are made of fermions (leptons and quarks), so every matter has its anti-matter partner — a profound symmetry encoded in the quantum mechanical rules of fermion operators.

Consider a fermion system of 5 single-particle modes, the vacuum state is denoted as (modes are arranged from left to right)

(Ket[00000]):=(Subscript[(Ket[0]), 1])⊗(Subscript[(Ket[0]), 2])⊗(Subscript[(Ket[0]), 3])⊗(Subscript[(Ket[0]), 4])⊗(Subscript[(Ket[0]), 5]).

Any fermion occupation state can be constructed by applying fermion creation operators on the vacuum state

(Ket[10000])=(Subscript[c, 1]^†)(Ket[00000]),
(Ket[01000])=(Subscript[c, 2]^†)(Ket[00000]),
(Ket[00100])=(Subscript[c, 3]^†)(Ket[00000]),
(Ket[10100])=(Subscript[c, 1]^†)(Subscript[c, 3]^†)(Ket[00000]),
(Ket[10110])=(Subscript[c, 1]^†)(Subscript[c, 3]^†)(Subscript[c, 4]^†)(Ket[00000]).

If we add fermions one by one to the system, the order matters:

Starting from the vacuum state,

add fermion on mode 1 then mode 3 then mode 4:

(Ket[00000])
→(Subscript[c, 1]^†)(Ket[00000])=(Ket[10000])
→(Subscript[c, 3]^†)(Ket[10000])=(Subscript[c, 3]^†)(Subscript[c, 1]^†)(Ket[00000])=-(Subscript[c, 1]^†)(Subscript[c, 3]^†)(Ket[00000])=-(Ket[10100])
->-(Subscript[c, 4]^†)(Ket[10100])=-(Subscript[c, 4]^†)(Subscript[c, 1]^†)(Subscript[c, 3]^†)(Ket[00000])=-(Subscript[c, 1]^†)(Subscript[c, 3]^†)(Subscript[c, 4]^†)(Ket[00000])=-(Ket[10110])

add fermion on mode 1 then mode 4 then mode 3:

(Ket[00000])
→(Subscript[c, 1]^†)(Ket[00000])=(Ket[10000])
→(Subscript[c, 4]^†)(Ket[10000])=(Subscript[c, 4]^†)(Subscript[c, 1]^†)(Ket[00000])=-(Subscript[c, 1]^†)(Subscript[c, 4]^†)(Ket[00000])=-(Ket[10010])
->-(Subscript[c, 3]^†)(Ket[10010])=-(Subscript[c, 3]^†)(Subscript[c, 1]^†)(Subscript[c, 4]^†)(Ket[00000])=(Subscript[c, 1]^†)(Subscript[c, 3]^†)(Subscript[c, 4]^†)(Ket[00000])=(Ket[10110])

The resulting states are differed by a minus sign, due to the exchange of fermions between mode 3 and 4 effectively comparing the two processes.

Jordan-Wigner Transformation

State correspondence

Qubit: (Ket[(Subscript[a, 1])(Subscript[a, 2])…]) with (Subscript[a, i])=0,1 — the ith qubit in the 0 or 1 state.

Fermion: (Ket[(Subscript[n, 1])(Subscript[n, 2])…]) with (Subscript[n, i])=0,1 — the ith mode being empty or occupied.

Operator correspondence: The fermion operators can be written in terms of Pauli operators as

Verify that the construction Eq. () satisfies the fermion anticommutation relation Eq. (), therefore justified to be called fermion operators.

First check Eq. (),

(Subscript[c, i])(Subscript[c, i])=(∏_(j<i))(Subscript[Z, j])(((Subscript[X, i])+i (Subscript[Y, i]))/2)(∏_(j<i))(Subscript[Z, j])(((Subscript[X, i])+i (Subscript[Y, i]))/2)
=(∏_(j<i))((Subscript[Z, j])(Subscript[Z, j]))((((Subscript[X, i])+i (Subscript[Y, i]))/2)(((Subscript[X, i])+i (Subscript[Y, i]))/2))
=(((Subscript[I, i])-(Subscript[Z, i])+(Subscript[Z, i])-(Subscript[I, i]))/4)=0.

Then (Subscript[c, i]^†)(Subscript[c, i]^†)=(SuperDagger[((Subscript[c, i])(Subscript[c, i]))])=(SuperDagger[0])=0. So Eq. () is verified.

Check Eq. (). With out lost of generality, we may assume i<j.

(Subscript[c, i]^((†)))(Subscript[c, j]^((†)))=(∏_(k<i))(Subscript[Z, k])(((Subscript[X, i])±i (Subscript[Y, i]))/2)(∏_(l<j))(Subscript[Z, l])(((Subscript[X, j])±i (Subscript[Y, j]))/2).

Due to the setting of i<j,

(((Subscript[X, i])±i (Subscript[Y, i]))/2)(∏_(l<j))(Subscript[Z, l])=-(∏_(l<j))(Subscript[Z, l])(((Subscript[X, i])±i (Subscript[Y, i]))/2),
(∏_(k<i))(Subscript[Z, k])(((Subscript[X, j])±i (Subscript[Y, j]))/2)=(((Subscript[X, j])±i (Subscript[Y, j]))/2)(∏_(k<i))(Subscript[Z, k]),

i.e., the (Subscript[X, i])±i (Subscript[Y, i]) operator anticommute with the Jordan-Wigner string from j, but (Subscript[X, j])±i (Subscript[Y, j]) operator commute with the Jordan-Wigner string from i. So only one minus sign will be produced in the commutation process

(Subscript[c, i]^((†)))(Subscript[c, j]^((†)))=-(∏_(k<i))(Subscript[Z, k])(∏_(l<j))(Subscript[Z, l])(((Subscript[X, j])±i (Subscript[Y, j]))/2)(((Subscript[X, i])±i (Subscript[Y, i]))/2)
=-(∏_(l<j))(Subscript[Z, l])(((Subscript[X, j])±i (Subscript[Y, j]))/2)(∏_(k<i))(Subscript[Z, k])(((Subscript[X, i])±i (Subscript[Y, i]))/2)
=-(Subscript[c, j]^((†)))(Subscript[c, i]^((†))).

Check Eq. ().

(Subscript[c, i])(Subscript[c, i]^†)=(∏_(j<i))(Subscript[Z, j])(((Subscript[X, i])+i (Subscript[Y, i]))/2)(∏_(j<i))(Subscript[Z, j])(((Subscript[X, i])-i (Subscript[Y, i]))/2)
=(∏_(j<i))((Subscript[Z, j])(Subscript[Z, j]))((((Subscript[X, i])+i (Subscript[Y, i]))/2)(((Subscript[X, i])-i (Subscript[Y, i]))/2))
=(((Subscript[I, i])+(Subscript[Z, i])+(Subscript[Z, i])+(Subscript[I, i]))/4)=((𝟙+(Subscript[Z, i]))/2).

(Subscript[c, i]^†)(Subscript[c, i])=(∏_(j<i))(Subscript[Z, j])(((Subscript[X, i])-i (Subscript[Y, i]))/2)(∏_(j<i))(Subscript[Z, j])(((Subscript[X, i])+i (Subscript[Y, i]))/2)
=(∏_(j<i))((Subscript[Z, j])(Subscript[Z, j]))((((Subscript[X, i])-i (Subscript[Y, i]))/2)(((Subscript[X, i])+i (Subscript[Y, i]))/2))
=(((Subscript[I, i])-(Subscript[Z, i])-(Subscript[Z, i])+(Subscript[I, i]))/4)=((𝟙-(Subscript[Z, i]))/2).

Therefore,

(Subscript[c, i])(Subscript[c, i]^†)+(Subscript[c, i]^†)(Subscript[c, i])=((𝟙+(Subscript[Z, i]))/2)+((𝟙-(Subscript[Z, i]))/2)=𝟙.

The product (∏_(j<i))(Subscript[Z, j]) is called the Jordan-Wigner string, a non-local operator necessary to ensure the correct fermion anticommutation relation is reproduced.

Comment: The non-locality of the Jordan-Wigner string has sparked deep reflection among some physicists.

If we simply accept the anticommutation relations of fermionic operators as the defining property of fermions, then we can build a quantum many-body theory of fermions, without questioning their origin.

However, if we consider the universe to be, at the fundamental level, a quantum computer—where every physical phenomenon ultimately arises from some computable quantum algorithm on foundational qubits—then it becomes difficult to imagine why all matter particles in the universe, as fermions, would have to be realized in such a complex, non-local manner. This line of thought has led to the discovery of topological order, a theoretical framework that provides a local physical origin for the non-local nature of fermions.

Matrix representation: Despite of the complexity, Eq. () actually provides a systematic method to represent fermion operators as matrices, given that we know how to represent each Pauli matrix and how to tensor product them together.

Quantum Many-Body Problems

Overview

What are the Problems?

Input: Hamiltonian H(g), parametrized by g=((Subscript[g, 1]),(Subscript[g, 2]),…),

Output: ground state (Ket[(Subscript[Ψ, 0])]) of H(g) for any g, i.e.,

such that given any observable O, one can predict its ground state expectation value ⟨O⟩=(Bra[(Subscript[Ψ, 0])])O(Ket[(Subscript[Ψ, 0])]).

Why care about ground state? — The ground state reveals properties of many-body systems at low temperature, which are crucial for understanding quantum phases (quantum magnet, superconductor, superfluid, topological insulator) and quantum phase transitions.

What are the Challenges?

Quantum Entanglement: When particle interact with each other, they become quantum entangled with each other. Entangled quantum systems are fundamentally non-decomposable due to non-local quantum information sharing among all the entangled parties.

Curse of Dimensionality: The Hilbert space dimension of quantum many-body system grows exponentially with the number of particles (or qubits) [Recall: N-qubit system Hilbert space dimension is (2^N)], making it exponentially hard to solve for the ground state (or even just to write down the solution).

Exact Diagonalization: A Brute-Force Approach

Algorithm

Input: Hamiltonian H (we will omit its parameters for now).

Step 1: represent the Hamiltonian H as a matrix.

Step 2: solve the eigen problem to find the eigenvalues (energy levels) and eigenvectors (corresponding states) of H

H(Ket[(Subscript[Ψ, n])])=(Subscript[E, n])(Ket[(Subscript[Ψ, n])]).

Step 3: find the lowest energy eigenvalue (Subscript[E, 0]), and pick out its corresponding eigenstate (Ket[(Subscript[Ψ, 0])]).

Output: (Ket[(Subscript[Ψ, 0])]) will be a ground state of H.

Explanation

(Ket[Ψ])=(∑_n)(Subscript[c, n])(Ket[(Subscript[Ψ, n])]),

(Bra[Ψ])H(Ket[Ψ])=(∑_n)((|(Subscript[c, n])|)^2)(Subscript[E, n])≥(∑_n)((|(Subscript[c, n])|)^2)(Subscript[E, 0])=(Subscript[E, 0])(∑_n)((|(Subscript[c, n])|)^2)=(Subscript[E, 0]).

Transverse-Field Ising Model

Model Hamiltonian

Graphics[{Circle[],Translate[{EdgeForm[Black],White,Disk[{0,0},0.06]},#]&/@CirclePoints[{1,π/6},12],Text[Style[StringTemplate["\!\(i=`1`\)"][#1],12],ReIm@(I Exp[-I π/6 #2]),{-Sign[#2],-1}]&@@@{{"N",0},{"1",1},{"2",2},{"N-1",-1},{"N-2",-2}},{AbsolutePointSize[1.5],Point[ReIm[# 1.2Exp[I 0.04π Range[-2,2]]]]&/@{1,-1}}},ImageSize->{Automatic,120}]

Hamiltonian

with the understanding that (Subscript[Z, i+N])=(Subscript[Z, i]) under periodic boundary condition.

Parameters:

Ising coupling strength J: large J favors all qubits to be either all 0 or all 1, such that (Subscript[Z, i])(Subscript[Z, i+1])=+1 for all neighbors, such that the Ising energy is minimized.

Transverse field strength h: large h favors all qubits to be in superposition states (Ket[+])∝(Ket[0])+(Ket[1]), such that (Subscript[X, i])=+1 for all qubits independently, such that the transverse field energy is minimized.

However, only their ratio J/h matters, as the overall energy scale can always be redefined.

(Subscript[ℤ, 2]) Ising Symmetry

Under (Subscript[ℤ, 2]) symmetry transformation,

(Subscript[Z, i])->U (Subscript[Z, i])(SuperDagger[U])=-(Subscript[Z, i]),
(Subscript[X, i])->U (Subscript[X, i])(SuperDagger[U])=(Subscript[X, i]),

the Ising spin operator (Subscript[Z, i]) gets a minus sign (spin flipped), but the Hamiltonian H in Eq. () remains unchanged.

In other words, the system is symmetric ⇔ the Hamiltonian H commutes with the symmetry operator U

for any parameters J and h in H.

The symmetry is named by the mathematical structure of a (Subscript[ℤ, 2]) group, which consists of two group elements (Subscript[ℤ, 2])={𝟙,U}, with the following group multiplication rules

({{·, 𝟙, U}, {𝟙, 𝟙, U}, {U, U, 𝟙}})  ≅  ({{+, 0, 1}, {0, 0, 1}, {1, 1, 0}}),

which is isomorphic to the mod-2 addition of (Subscript[ℤ, 2]) integers {0,1} (hence the name).

Two Limits

Strong coupling limit J/h→∞: ferromagnetic, ordered phase.

Ground states: two fold degenerated

(Ket[(Subscript[Ψ, 0])])∈span{(Ket[000…]),(Ket[111…])}.

No matter which state we choose, the ground state will spontaneously break the (Subscript[ℤ, 2]) symmetry.

Spontaneous symmetry breaking: the Hamiltonian H has a symmetry, but all its short-range-entangled ground states necessarily breaks the symmetry.

If a discrete symmetry group G is spontaneously broken to its subgroup G', the ground state degeneracy (GSD) will be given by

GSD=((|G|)/(|G'|)),

where |G| denotes the order (number of elements) of a group G.

=> In our case, (Subscript[ℤ, 2])→(Subscript[ℤ, 1]), GSD=2/1=2.

Weak coupling limit J/h→0: paramagnetic, disordered phase.

Ground state: unique

(Ket[(Subscript[Ψ, 0])])=(Ket[+++…])=(1/(2^(N/2)))(∑_([a]))(Ket[[a]]),

as an equal-amplitude superposition of all possible Ising (bit-string) configurations [a]=(Subscript[a, 1])(Subscript[a, 2])… ((Subscript[a, i])=0,1).

Symmetric: the (Subscript[ℤ, 2]) symmetry is preserved (no symmetry breaking).

Representing the Hamiltonian

Let us first pick a finite system size N, and enumerate terms in our Hamiltonian.

Block[{Ν=4,J,h,H},
H=-J Sum[RotateLeft[PadLeft[P["Z","Z"],Ν,"I"],i],{i,Ν}]-h Sum[RotateLeft[PadLeft[P["X"],Ν,"I"],i],{i,Ν}]]

As a symbolic programing language, Mathematica allows us to manipulate mathematical symbols directly, which makes it well-suited for algebraic manipulation, equation solving, and calculus.

The key mechanism for the symbolic processing is its rule-based replacement system, which works by defining rules that determines how symbols are matched and replaced according to specific patterns.

Here is how it works:

First define a rule set called “represent”, that replaces I,X,Y,Z by their corresponding Pauli matrices.

Block[{represent},
represent={
"I"->SparseArray@PauliMatrix[0],
"X"->SparseArray@PauliMatrix[1],
"Y"->SparseArray@PauliMatrix[2],
"Z"->SparseArray@PauliMatrix[3]};
MatrixForm/@{"I","X","Y","Z"}//.represent]

Then extend the rule set by rules for Pauli operators

Block[{represent},
represent={,
P[a_]:>a,
P[a_,b__]:>KroneckerProduct[a,P[b]]};
MatrixForm/@{P["X"],P["Z","Z"],P["X","I","Z"]}//.represent]

Block[{Ν=4,J,h,H,represent},
represent=;
H=-J Sum[RotateLeft[PadLeft[P["Z","Z"],Ν,"I"],i],{i,Ν}]-h Sum[RotateLeft[PadLeft[P["X"],Ν,"I"],i],{i,Ν}];
MatrixForm[H]//.represent]

Solving the Eigen Problem

Eigen solver. For numerical matrix, the eigen problems can be solved by an algorithm called Arnoldi iteration. One can request the number of eigen values and states to be calculated.

Parametrization. For this purpose, we substitute numerical values for J and h. Since only J/h matters, we can choose to parametrize them as

J=g, h=1-g,

with g∈[0,1], such that there is only a single parameter to tune, and the overall energy scale is bounded between ±N.

Block[{Ν=4,g=0.,H,represent,Es,Ψs},
represent=;
H[J_,h_]:=;
{Es,Ψs}=Eigensystem[N@H[g,1-g],2,Method->{"Arnoldi","Shift"->-Ν-0.1}]]

Manipulate[DynamicModule[{Ν=4,H,represent,Es,Ψs,basis},
represent={"I"->SparseArray[PauliMatrix[0]],"X"->SparseArray[PauliMatrix[1]],"Y"->SparseArray[PauliMatrix[2]],"Z"->SparseArray[PauliMatrix[3]],P[a_]:>a,P[a_,b__]:>KroneckerProduct[a,P[b]]};
H[J_,h_]:=-J Sum[RotateLeft[PadLeft[P["Z","Z"],Ν,"I"],i],{i,Ν}]-h Sum[RotateLeft[PadLeft[P["X"],Ν,"I"],i],{i,Ν}]//. represent;
{Es,Ψs}=Eigensystem[N@H[g,1-g],2,Method->{"Arnoldi","Shift"->-Ν-0.1}];
basis=Row[Join[{"\[LeftBracketingBar]"},#,{"⟩"}]]&/@Tuples[{Style[0,Blue],Style[1,Red]},Ν];
Style[Column[{StringTemplate["\!\(J/h=`1`\)"][If[PossibleZeroQ[1-g],∞,g/1-g]],Grid[Thread[{Es,Chop[#,0.01].basis&/@Ψs}],
ItemStyle->{Automatic,Automatic,{{{2,2},{1,2}}->Opacity[Exp[2(Es[[1]]-Es[[2]])]]}}]},Dividers->Center],"DisplayFormula"]],{{g,1},0,1},Paneled->False]

Strictly speaking, the ground state degeneracy splits as long as the system is tuned away from the J/h->∞ limit. However, the splitting is small in the strong coupling regime.

It seems that the ground state interpolates between the two limits continuously.

Energy Gap Scaling

Block[{Ν=12,H,represent,gs,Es},
represent={"I"->SparseArray[PauliMatrix[0]],"X"->SparseArray[PauliMatrix[1]],"Y"->SparseArray[PauliMatrix[2]],"Z"->SparseArray[PauliMatrix[3]],
P[a_]:>a,P[a_,b__]:>KroneckerProduct[a,P[b]]};
H[J_,h_]:=-J Sum[RotateLeft[PadLeft[P["Z","Z"],Ν,"I"],i],{i,Ν}]-h Sum[RotateLeft[PadLeft[P["X"],Ν,"I"],i],{i,Ν}]//. represent;
gs=Range[0,1,0.05];
Es=Table[Eigenvalues[N[H[g,1-g]],2,Method->{"Arnoldi","Shift"->-Ν-0.1}]/Ν,{g,gs}];
ListLinePlot[Thread@{gs,#}&/@Thread@Es,Epilog->{{Red,Text["\!\(E\_0\)",{1/6,-0.93}]},{Yellow,Text["\!\(E\_1\)",{1/6,-0.68}]}},PlotStyle->{Red,Directive[Yellow,Dashed]},PlotLabel->StringTemplate["\!\(N=`1`\)"][Ν],FrameTicks->{{Automatic,Automatic},{{{0,0},{1/6,0.2},{1/3,0.5},{1/2,1},{2/3,2},{5/6,5},{1,∞}},None}},FrameLabel->{"\!\(J/h\)","\!\(E\_n/N\)"},ImageSize->200]]

Δ=(Subscript[E, 1])-(Subscript[E, 0]).

In the symmetry breaking phase (strong coupling), e.g. J/h=2, the energy gap decays to zero with system size N exponentially

Block[{Ns,Δs},
Ns=Range[4,12];
Δs=Table[Block[{g=2/3,represent,H},
represent={"I"->SparseArray[PauliMatrix[0]],"X"->SparseArray[PauliMatrix[1]],"Y"->SparseArray[PauliMatrix[2]],"Z"->SparseArray[PauliMatrix[3]],
P[a_]:>a,P[a_,b__]:>KroneckerProduct[a,P[b]]};
H[J_,h_]:=-J Sum[RotateLeft[PadLeft[P["Z","Z"],Ν,"I"],i],{i,Ν}]-h Sum[RotateLeft[PadLeft[P["X"],Ν,"I"],i],{i,Ν}]//. represent;
First@Differences@Eigenvalues[N[H[g,1-g]],2,Method->{"Arnoldi","Shift"->-Ν-0.1}]],{Ν,Ns}];
ListLogPlot[Thread[{Ns,Δs}],PlotStyle->Black,Joined->True,Mesh->All,GridLines->Automatic,FrameLabel->{"\!\( N\)","\!\(Δ\)"},ImageSize->{Automatic,140}]
]

In the symmetric phase (weak coupling), e.g. J/h=1/2, the energy gap tends to a constant in the large N limit.

Block[{Ns,Δs},
Ns=Range[4,12];
Δs=Table[Block[{g=1/3,represent,H},
represent={"I"->SparseArray[PauliMatrix[0]],"X"->SparseArray[PauliMatrix[1]],"Y"->SparseArray[PauliMatrix[2]],"Z"->SparseArray[PauliMatrix[3]],
P[a_]:>a,P[a_,b__]:>KroneckerProduct[a,P[b]]};
H[J_,h_]:=-J Sum[RotateLeft[PadLeft[P["Z","Z"],Ν,"I"],i],{i,Ν}]-h Sum[RotateLeft[PadLeft[P["X"],Ν,"I"],i],{i,Ν}]//. represent;
First@Differences@Eigenvalues[N[H[g,1-g]],2,Method->{"Arnoldi","Shift"->-Ν-0.1}]],{Ν,Ns}];
ListPlot[Thread[{Ns,Δs}],PlotStyle->Black,Joined->True,Mesh->All,GridLines->Automatic,FrameLabel->{"\!\( N\)","\!\(Δ\)"},ImageSize->{Automatic,140}]
]

Physical Observables

⟨O⟩:=(Bra[(Subscript[Ψ, 0])])O(Ket[(Subscript[Ψ, 0])]).

Correlation function: used to indicate the ordered phase.

It characterizes the correlation between Ising spins (Subscript[Z, i]) and (Subscript[Z, i+d]) separated by distance d on the lattice. Strong correlation indicates the spins are ordered.

Long-distance correlation function v.s. J/h ratio, indicating phase transition (transition gets sharper as N→∞)

Block[{Ν=12,represent,H,ZZ,gs,Ψs,Cs},
represent={"I"->SparseArray[PauliMatrix[0]],"X"->SparseArray[PauliMatrix[1]],"Y"->SparseArray[PauliMatrix[2]],"Z"->SparseArray[PauliMatrix[3]],
P[a_]:>a,P[a_,b__]:>KroneckerProduct[a,P[b]]};
H[J_,h_]:=-J Sum[RotateLeft[PadLeft[P["Z","Z"],Ν,"I"],i],{i,Ν}]-h Sum[RotateLeft[PadLeft[P["X"],Ν,"I"],i],{i,Ν}]//. represent;
ZZ=Join[#,#]&@PadLeft[P["Z"],Ν/2,"I"]//.represent;
gs=Range[0,1,0.05];
Ψs=Table[First@Eigenvectors[N[H[g,1-g]],1,Method->{"Arnoldi","Shift"->-Ν-0.1}],{g,gs}];
Cs=Conjugate[#].ZZ.#&/@Ψs;
ListLinePlot[Thread@{gs,Cs},PlotRange->{0,1},PlotRangePadding->0.05,GridLines->{{1/2},{0,1}},PlotStyle->Black,PlotLabel->StringTemplate["\!\(N=`1`\)"][Ν],FrameTicks->{{Automatic,Automatic},{{{0,0},{1/6,0.2},{1/3,0.5},{1/2,1},{2/3,2},{5/6,5},{1,∞}},None}},FrameLabel->{"\!\(J/h\)","\!\(C(N/2)\)"},ImageSize->180]
]

In the symmetry breaking phase (strong coupling), e.g. J/h=2, the correlation function approaches to constant at long range

C(d)∼const.

Block[{Ns,Cs},
Ns=Range[4,12,2];
Cs=Table[Block[{g=2/3,represent,H,ZZ},
represent={"I"->SparseArray[PauliMatrix[0]],"X"->SparseArray[PauliMatrix[1]],"Y"->SparseArray[PauliMatrix[2]],"Z"->SparseArray[PauliMatrix[3]],
P[a_]:>a,P[a_,b__]:>KroneckerProduct[a,P[b]]};
H[J_,h_]:=-J Sum[RotateLeft[PadLeft[P["Z","Z"],Ν,"I"],i],{i,Ν}]-h Sum[RotateLeft[PadLeft[P["X"],Ν,"I"],i],{i,Ν}]//. represent;
ZZ=Join[#,#]&@PadLeft[P["Z"],Ν/2,"I"]//.represent;
Conjugate[#].ZZ.#&@First@Eigenvectors[N[H[g,1-g]],1,Method->{"Arnoldi","Shift"->-Ν-0.1}]],{Ν,Ns}];
ListPlot[Thread[{Ns,Cs}],PlotStyle->Black,Joined->True,Mesh->All,GridLines->Automatic,FrameLabel->{"\!\( N\)","\!\(C(N/2)\)"},ImageSize->{Automatic,140}]]

In the symmetric phase (weak coupling), e.g. J/h=1/2, the correlation function decays to zero exponentially with the distance d.

C(d)∼(e^(-d/ξ)),

the decay length scale ξ is called the correlation length.

Block[{Ns,Cs},
Ns=Range[4,12,2];
Cs=Table[Block[{g=1/3,represent,H,ZZ},
represent={"I"->SparseArray[PauliMatrix[0]],"X"->SparseArray[PauliMatrix[1]],"Y"->SparseArray[PauliMatrix[2]],"Z"->SparseArray[PauliMatrix[3]],
P[a_]:>a,P[a_,b__]:>KroneckerProduct[a,P[b]]};
H[J_,h_]:=-J Sum[RotateLeft[PadLeft[P["Z","Z"],Ν,"I"],i],{i,Ν}]-h Sum[RotateLeft[PadLeft[P["X"],Ν,"I"],i],{i,Ν}]//. represent;
ZZ=Join[#,#]&@PadLeft[P["Z"],Ν/2,"I"]//.represent;
Conjugate[#].ZZ.#&@First@Eigenvectors[N[H[g,1-g]],1,Method->{"Arnoldi","Shift"->-Ν-0.1}]],{Ν,Ns}];
ListLogPlot[Thread[{Ns,Cs}],PlotStyle->Black,Joined->True,Mesh->All,GridLines->Automatic,FrameLabel->{"\!\( N\)","\!\(C(N/2)\)"},ImageSize->{Automatic,140}]]

String operator: used to indicate the disordered phase.

String operator is a segment of symmetry operator of length d, which flips all the Ising spins in the segment. Its expectation value measure the probability to see the segment of spin being flipped under quantum fluctuation in the ground state. In the symmetric (disordered) phase, such probability should be high.

Long-range string operator v.s. J/h ratio, indicating phase transition (transition gets sharper as N→∞)

Block[{Ν=12,represent,H,ΠX,gs,Ψs,Ss},
represent={"I"->SparseArray[PauliMatrix[0]],"X"->SparseArray[PauliMatrix[1]],"Y"->SparseArray[PauliMatrix[2]],"Z"->SparseArray[PauliMatrix[3]],
P[a_]:>a,P[a_,b__]:>KroneckerProduct[a,P[b]]};
H[J_,h_]:=-J Sum[RotateLeft[PadLeft[P["Z","Z"],Ν,"I"],i],{i,Ν}]-h Sum[RotateLeft[PadLeft[P["X"],Ν,"I"],i],{i,Ν}]//. represent;
ΠX=PadLeft[P@@ConstantArray["X",Ν/2],Ν,"I"]//.represent;
gs=Range[0,1,0.05];
Ψs=Table[First@Eigenvectors[N[H[g,1-g]],1,Method->{"Arnoldi","Shift"->-Ν-0.1}],{g,gs}];
Ss=Conjugate[#].ΠX.#&/@Ψs;
ListLinePlot[Thread@{gs,Ss},PlotRange->{0,1},PlotRangePadding->0.05,GridLines->{{1/2},{0,1}},PlotStyle->Black,PlotLabel->StringTemplate["\!\(N=`1`\)"][Ν],FrameTicks->{{Automatic,Automatic},{{{0,0},{1/6,0.2},{1/3,0.5},{1/2,1},{2/3,2},{5/6,5},{1,∞}},None}},FrameLabel->{"\!\(J/h\)","\!\(S(N/2)\)"},ImageSize->180]
]

In the symmetry breaking phase (strong coupling), e.g. J/h=2, the string operator decays to zero exponentially with the distance d.

S(d)∼(e^(-d/ξ)).

Block[{Ns,Ss},
Ns=Range[4,12,2];
Ss=Table[Block[{g=2/3,represent,H,ΠX},
represent={"I"->SparseArray[PauliMatrix[0]],"X"->SparseArray[PauliMatrix[1]],"Y"->SparseArray[PauliMatrix[2]],"Z"->SparseArray[PauliMatrix[3]],
P[a_]:>a,P[a_,b__]:>KroneckerProduct[a,P[b]]};
H[J_,h_]:=-J Sum[RotateLeft[PadLeft[P["Z","Z"],Ν,"I"],i],{i,Ν}]-h Sum[RotateLeft[PadLeft[P["X"],Ν,"I"],i],{i,Ν}]//. represent;
ΠX=PadLeft[P@@ConstantArray["X",Ν/2],Ν,"I"]//.represent;
Conjugate[#].ΠX.#&@First@Eigenvectors[N[H[g,1-g]],1,Method->{"Arnoldi","Shift"->-Ν-0.1}]],{Ν,Ns}];
ListLogPlot[Thread[{Ns,Ss}],PlotStyle->Black,Joined->True,Mesh->All,GridLines->Automatic,FrameLabel->{"\!\( N\)","\!\(S(N/2)\)"},ImageSize->{Automatic,140}]]

In the symmetric phase (weak coupling), e.g. J/h=1/2, the string operator approaches to constant at long range.

S(d)∼const.

Block[{Ns,Ss},
Ns=Range[4,12,2];
Ss=Table[Block[{g=1/3,represent,H,ΠX},
represent={"I"->SparseArray[PauliMatrix[0]],"X"->SparseArray[PauliMatrix[1]],"Y"->SparseArray[PauliMatrix[2]],"Z"->SparseArray[PauliMatrix[3]],
P[a_]:>a,P[a_,b__]:>KroneckerProduct[a,P[b]]};
H[J_,h_]:=-J Sum[RotateLeft[PadLeft[P["Z","Z"],Ν,"I"],i],{i,Ν}]-h Sum[RotateLeft[PadLeft[P["X"],Ν,"I"],i],{i,Ν}]//. represent;
ΠX=PadLeft[P@@ConstantArray["X",Ν/2],Ν,"I"]//.represent;
Conjugate[#].ΠX.#&@First@Eigenvectors[N[H[g,1-g]],1,Method->{"Arnoldi","Shift"->-Ν-0.1}]],{Ν,Ns}];
ListPlot[Thread[{Ns,Ss}],PlotStyle->Black,Joined->True,Mesh->All,GridLines->Automatic,FrameLabel->{"\!\( N\)","\!\(S(N/2)\)"},ImageSize->{Automatic,140}]]

Kramers-Wannier Duality

Block[{R=2.25,h1=0.5,h2=0.5,a1=0.1,a2=0.2,lt=0.85,latt,ah},
latt[col_,s_]:={{col,Line[{#,0}&/@{-R,R}]},Translate[{EdgeForm[col],White,Disk[{0,0},a1]},{#,0}]&/@Cases[Range[-3,3]+s,x_/;-R<x<R]};
ah={Arrowheads[{0.03}],Arrow[{0,0.6#}&/@{h2,-h2}]};
Graphics[{Translate[{Translate[{{Lighter[Red,lt],AbsoluteThickness[16],Line[{#,h1}&/@{-a2,1+a2}]},latt[Red,0],{Red,Text["\!\(Z\_i\)",{0,h1}],Text["\!\(Z\_\(i+1\)\)",{1,h1}]}},{0,h2}],Translate[ah,{1/2,0}],Translate[{{Lighter[Blue,lt],AbsoluteThickness[16],Line[{#,-h1}&/@{1/2-a2,1/2+a2}]},latt[Blue,1/2],{Blue,Text["\!\(X\_\(i+1/2\)\)",{1/2,-h1}]}},{0,-h2}]},{-3,0}],
Translate[{Translate[{{Lighter[Red,lt],AbsoluteThickness[16],Line[{#,h1}&/@{-a2,a2}]},latt[Red,0],{Red,Text["\!\(X\_i\)",{0,h1}]}},{0,h2}],
Translate[ah,{0,0}],Translate[{{Lighter[Blue,lt],AbsoluteThickness[16],Line[{#,-h1}&/@{-1/2-a2,1/2+a2}]},latt[Blue,1/2],{Blue,Text["\!\(Z\_\(i-1/2\)\)",{-1/2,-h1},{0.2,0}],Text["\!\(Z\_\(i+1/2\)\)",{1/2,-h1},{-0.2,0}]}},{0,-h2}]},{3,0}]},ImageSize->300]]

It maps the Pauli operators acting the original qubits (in red) to the Pauli operators acting a new set of qubits (in blue) on the dual lattice (shifted by 1/2 lattice constant).

Requirements:

The algebraic relations (such as commutation relations) among these operators are preserved under the mapping.

The locality of the Hamiltonian is also preserved.

({{J/h, ↔, h/J}, {lattice, ↔, dual-lattice}, {order, ↔, disorder}, {C(d), , S(d)}, {disorder, ↔, order}, {S(d), , C(d)}})

Show that the correlation function C(d) and string operator S(d) maps to each other under Kramers-Wannier duality.

C(d)=⟨(Subscript[Z, i])(Subscript[Z, i+d])⟩=⟨(∏_(j=i)^(i+d-1))(Subscript[Z, j])(Subscript[Z, j+1])⟩
→⟨(∏_(j=i)^(i+d-1))(Subscript[X, j+1/2])⟩=S(d),

S(d)=⟨(∏_(j=i)^(i+d-1))(Subscript[X, j])⟩
→⟨(∏_(j=i)^(i+d-1))(Subscript[Z, j-1/2])(Subscript[Z, j+1/2])⟩=⟨(Subscript[Z, i-1/2])(Subscript[Z, i+d-1/2])⟩=C(d),

Implication: under Kramers-Wannier duality, the ordered (disordered) phase of the original Ising model maps to the disordered (ordered) phase of the dual Ising model, but both models are identical, so if there is a single phase transition separating the ordered and disordered phases, it must happen at

This is the critical point of the transverse-field Ising model.

Bose-Hubbard Model

Model Hamiltonian

Hamiltonian

where

(Subscript[b, i]^†) and (Subscript[b, i]) are boson creation and annihilation operators on site-i, with (Subscript[b, i+N])=(Subscript[b, i]) given the periodic boundary condition,

(Subscript[n, i])=(Subscript[b, i]^†)(Subscript[b, i]) denotes the boson number operator (we will omit the hat of (OverHat[n]) for simplicity).

Parameters:

Hopping amplitude t: measures boson mobility on the lattice. (Subscript[b, i+1]^†)(Subscript[b, i]) describes the process of moving a boson from site-i to the neighboring site-(i+1), and the Hermitian conjugate term (Subscript[b, i]^†)(Subscript[b, i+1]) describes the reverse process.

On-site interaction strength U: also known as the Hubbard interaction, can be attractive (U<0) or reclusive (U>0). (Subscript[n, i])((Subscript[n, i])-1)/2 counts the number of pairs that can form among (Subscript[n, i]) bosons on a site, as interaction is pair-wise.

Chemical potential μ: adjust to set the number of particles.

👉 Remark: If we fix the total boson number (Subscript[N, b]) (i.e., assuming a canonical ensemble), we can drop the chemical potential term.

U(1) Symmetry

Under the U(1) symmetry transformation,

(Subscript[b, i])→U(θ)(Subscript[b, i])(SuperDagger[U(θ)])=(e^(-i θ))(Subscript[b, i]),
(Subscript[b, i]^†)→U(θ)(Subscript[b, i]^†)(SuperDagger[U(θ)])=(e^(i θ))(Subscript[b, i]^†).

Prove Eq. () given the definition Eq. ().

(Subscript[n, i])(Subscript[b, i])=(Subscript[b, i]^†)(Subscript[b, i])(Subscript[b, i])=((Subscript[b, i])(Subscript[b, i]^†)-𝟙)(Subscript[b, i])=(Subscript[b, i])(Subscript[b, i]^†)(Subscript[b, i])-(Subscript[b, i])=(Subscript[b, i])((Subscript[b, i]^†)(Subscript[b, i])-𝟙)=(Subscript[b, i])((Subscript[n, i])-𝟙),

(Subscript[N, b]) (Subscript[b, i])=((∑_i)(Subscript[n, i]))(Subscript[b, i])=(Subscript[b, i])((∑_i)(Subscript[n, i])-𝟙)=(Subscript[b, i])((Subscript[N, b])-𝟙).

U(θ)(Subscript[b, i])(SuperDagger[U(θ)])=(e^(i θ (Subscript[N, b])))(Subscript[b, i])(e^(-i θ (Subscript[N, b])))
=(Subscript[b, i]) (e^(i θ ((Subscript[N, b])-𝟙)))(e^(-i θ (Subscript[N, b])))
=(Subscript[b, i]) (e^(i θ ((Subscript[N, b])-𝟙)))(e^(-i θ (Subscript[N, b])))
=(Subscript[b, i]) (e^(i θ 𝟙))
=(e^(-i θ))(Subscript[b, i]) .

The Bose-Hubbard model Hamiltonian H in Eq. () is invariant under such opposite phase rotations of the boson creation and annihilation operators (as they always appear in pairs in H).

The U(1) symmetry is associated with the conservation law of the boson number conservation, as the symmetry implies

the charge operator (Subscript[N, b]) will not evolve in time, and is therefore conserved.

Because (Subscript[N, b]) is conserved, we can always restrict the Hilbert space to a fixed eigen-subspace of (Subscript[N, b]), such that (Subscript[N, b]) becomes a number (the eigenvalue) other than an operator.

Hilbert Space

(Subscript[D, b])=((((Subscript[N, b])+N-1)!)/((Subscript[N, b])!(N-1)!)),

DynamicModule[{Νb=5,Ν=4,ball,div,sample,set},
ball=Graphics[{Red,Disk[]},ImageSize->{12,12}];
div=Graphics[{Gray,AbsoluteThickness[2],Line[{0,#}&/@{0,1}]},PlotRange->All,ImageSize->{12,12}];
sample[]:=RandomSample[Join[ConstantArray[ball,Νb],ConstantArray[div,Ν-1]]];
set=sample[];
EventHandler[Deploy[Dynamic[Grid[{set,Style[ToString[#],10,"Text"]&/@Range[Νb+Ν-1]}]]],
{"MouseClicked":>(set=sample[])}]]

Following the above idea, all basis states can be enumerated given N and (Subscript[N, b]). Here is a simple example:

Block[{Ν=3,Νb=3,basis},
basis=Ket[Differences[Join[{0},#,{Νb+Ν}]]-1]&/@Subsets[Range[Νb+Ν-1],{Ν-1}]]

The Hilbert space dimension grows exponentially with system size N at any fixed density (Subscript[N, b])/N. — Exactly solving the problem is hard for large systems.

Representing the Hamiltonian

We can use Mathematica to quickly build a symbolic system to implement the action of boson creation and annihilation operators on boson occupation states.

Block[{act,b,b†,𝟙,n},
(* Rules *)
act[b[i_],Ket[ns_]]:=Sqrt[ns[[i]]]Ket[ReplacePart[ns,i->ns[[i]]-1]];
act[b†[i_],Ket[ns_]]:=Sqrt[ns[[i]]+1]Ket[ReplacePart[ns,i->ns[[i]]+1]];
act[𝟙,st_]:=st;
act[op_,x_ st:(_Ket|_Plus)]:=x act[op,st];
act[x_ op:(𝟙|_b|_b†|_Composition|_Plus),st_]:=x act[op,st];
act[op1_@*op2_,st_]:=act[op1,act[op2,st]];
act[op1_+op2_,st_]:=act[op1,st]+act[op2,st];
act[op_,st1_+st2_]:=act[op,st1]+act[op,st2];
act[op_,0]:=0;
n[i_]:=b†[i]@*b[i];
(* Examples *)
Column[Row[{#1,"@",#2}]->act[#1,#2]&@@@{
{b[1],Ket[{0,0,0}]},
{b†[1],Ket[{0,0,0}]},
{b†[2]@*b[3],Ket[{0,4,3}]},
{n[2],Ket[{0,2,1}]+2Ket[{0,3,0}]}}]]

This enables us to apply the Hamiltonian Eq. () to each basis state and compute the resulting state.

Block[{Ν=3,Νb=3,basis,act,b,b†,𝟙,n,H},
basis=(Ket[Differences[Join[{0},#1,{Νb+Ν}]]-1]&)/@Subsets[Range[Νb+Ν-1],{Ν-1}];
act[b[i_],Ket[ns_]]:=Sqrt[ns[[i]]]Ket[ReplacePart[ns,i->ns[[i]]-1]];
act[b†[i_],Ket[ns_]]:=Sqrt[ns[[i]]+1]Ket[ReplacePart[ns,i->ns[[i]]+1]];
act[𝟙,st_]:=st;
act[op_,x_ st:(_Ket|_Plus)]:=x act[op,st];
act[x_ op:(𝟙|_b|_b†|_Composition|_Plus),st_]:=x act[op,st];
act[op1_@*op2_,st_]:=act[op1,act[op2,st]];
act[op1_+op2_,st_]:=act[op1,st]+act[op2,st];
act[op_,st1_+st2_]:=act[op,st1]+act[op,st2];
act[op_,0]:=0;
n[i_]:=b†[i]@*b[i];
H=-t Sum[(b†[Mod[i+1,Ν,1]]@*b[i]+b†[i]@*b[Mod[i+1,Ν,1]]),{i,Ν}]+(U/2)Sum[n[i]@*(n[i]-𝟙),{i,Ν}];
Column[Row[{"H","@",#}]->Expand[act[H,#]]&/@basis]
]

We can collect coefficients in front of the basis states and organize them into the matrix representation.

Block[{Ν=3,Νb=3,basis,id,act,b,b†,𝟙,n,H},
basis=(Ket[Differences[Join[{0},#1,{Νb+Ν}]]-1]&)/@Subsets[Range[Νb+Ν-1],{Ν-1}];
id=AssociationThread[basis->Range@Length@basis];
(* Rules *)
act[b[i_],Ket[ns_]]:=Sqrt[ns[[i]]]Ket[ReplacePart[ns,i->ns[[i]]-1]];
act[b†[i_],Ket[ns_]]:=Sqrt[ns[[i]]+1]Ket[ReplacePart[ns,i->ns[[i]]+1]];
act[𝟙,st_]:=st;
act[op_,x_ st:(_Ket|_Plus)]:=x act[op,st];
act[x_ op:(𝟙|_b|_b†|_Composition|_Plus),st_]:=x act[op,st];
act[op1_@*op2_,st_]:=act[op1,act[op2,st]];
act[op1_+op2_,st_]:=act[op1,st]+act[op2,st];
act[op_,st1_+st2_]:=act[op,st1]+act[op,st2];
act[op_,0]:=0;
n[i_]:=b†[i]@*b[i];
(* operator form *)
H=-t Sum[(b†[Mod[i+1,Ν,1]]@*b[i]+b†[i]@*b[Mod[i+1,Ν,1]]),{i,Ν}]+(U/2)Sum[n[i]@*(n[i]-𝟙),{i,Ν}];
(* matrix form *)
H=SparseArray[Flatten[Function[{i,dat},{i,#1}->#2&@@@dat]@@@Thread[{Values@id,Cases[Expand[act[H,#]],x_ st_Ket:>{id[st],x}]&/@basis}]],{#,#}&@Length@basis];
MatrixForm@H]

Solving the Eigen Problem

DynamicModule[{Ν=5,Νb=5,basis,id,act,b,b†,𝟙,n,rep,Ht,HU},
basis=(Ket[Differences[Join[{0},#1,{Νb+Ν}]]-1]&)/@Subsets[Range[Νb+Ν-1],{Ν-1}];
id=AssociationThread[basis->Range@Length@basis];
(* Rules *)
act[b[i_],Ket[ns_]]:=Sqrt[ns[[i]]]Ket[ReplacePart[ns,i->ns[[i]]-1]];
act[b†[i_],Ket[ns_]]:=Sqrt[ns[[i]]+1]Ket[ReplacePart[ns,i->ns[[i]]+1]];
act[𝟙,st_]:=st;
act[op_,x_ st:(_Ket|_Plus)]:=x act[op,st];
act[x_ op:(𝟙|_b|_b†|_Composition|_Plus),st_]:=x act[op,st];
act[op1_@*op2_,st_]:=act[op1,act[op2,st]];
act[op1_+op2_,st_]:=act[op1,st]+act[op2,st];
act[op_,st1_+st2_]:=act[op,st1]+act[op,st2];
act[op_,0]:=0;
n[i_]:=b†[i]@*b[i];
rep[op_]:=SparseArray[Flatten[Function[{i,dat},{i,#1}->#2&@@@dat]@@@Thread[{Values@id,Block[{expr,kets,cs},expr=Expand[act[op,#]];kets=DeleteDuplicates[Cases[expr,_Ket,{0,∞}]];If[Length[kets]==0,{},cs=Flatten@Normal@Rest@CoefficientArrays[expr,kets];
Thread@{id/@kets,cs}]]&/@basis}]],{#,#}&@Length@basis];
Ht=rep[-Sum[(b†[Mod[i+1,Ν,1]]@*b[i]+b†[i]@*b[Mod[i+1,Ν,1]]),{i,Ν}]];
HU=rep[Sum[n[i]@*(n[i]-𝟙),{i,Ν}]/2];
Manipulate[DynamicModule[{t=0.17(1-g),U=g,E0,Ψ0},
(* Eigen solver *)
{{E0},{Ψ0}}=Eigensystem[N@(t Ht+U HU),1,Method->{"Arnoldi","Shift"->-2Νb-0.1}];
Ψ0=Ψ0 Sign[Total[Ψ0]];(* gauging *)
(* Visualization of Ψ0 *)
Block[{terms},terms=KeySort[GroupBy[DeleteCases[Thread[{Round[Ψ0,0.001],basis}],{0.,_}],First->Last],Greater];
Style[Column[{StringTemplate["\!\(U/t=`1`\)"][If[PossibleZeroQ[t],∞,U/t]],
Column[Table[Row[{"+ ",coeff First@terms[coeff]/.Ket[ns_]:>Block[{a=3,site},site[n_]:={{EdgeForm[Directive[AbsoluteThickness[1],LightGray]],White,Disk[]},{Disk[(n-1)/n#,0.35]&/@CirclePoints[{1,π/2},n]}};Graphics[{Text["\[LeftBracketingBar]",{0.5a,0}],Text["⟩",{(Ν+0.8)a,0}],{LightGray,Line[{a#,0}&/@{1,Ν}]},Translate[site[#2],{a#1,0}]&@@@Thread@{Range@Length@ns,ns}},ImageSize->{Automatic,14}]],If[#==0,Nothing,StringTemplate[" + (`1` terms)"][#]]&@Length@Rest@terms[coeff]}],{coeff,Keys@terms}]]},Dividers->Center],"InlineFormula"]]],{g,0,1},Paneled->False]]

Weak Interaction regime (small U/t)

Bosons are delocalized across the entire lattice. Each boson can be found anywhere in the system, independent of the other bosons.

Strong Interaction regime (large U/t)

Bosons are localized. Each boson is confined to a lattice site, due to the strong on-site repulsive interaction that prevents bosons from hopping to neighboring sites.

There is an energy gap to boson excitations.

There is a quantum phase transition separating the weak and strong interacting phases, known as the Kosterlitz-Thouless (KT) transition. [2016 Nobel Prize in Physics]

Physical Observables

Boson Correlation

Block[{Ν=7,Νb=7,basis,id,act,b,b†,𝟙,n,rep,Ht,HU,bb},
basis=(Ket[Differences[Join[{0},#1,{Νb+Ν}]]-1]&)/@Subsets[Range[Νb+Ν-1],{Ν-1}];
id=AssociationThread[basis->Range@Length@basis];
(* Rules *)
act[b[i_],Ket[ns_]]:=Sqrt[ns[[i]]]Ket[ReplacePart[ns,i->ns[[i]]-1]];
act[b†[i_],Ket[ns_]]:=Sqrt[ns[[i]]+1]Ket[ReplacePart[ns,i->ns[[i]]+1]];
act[𝟙,st_]:=st;
act[op_,x_ st:(_Ket|_Plus)]:=x act[op,st];
act[x_ op:(𝟙|_b|_b†|_Composition|_Plus),st_]:=x act[op,st];
act[op1_@*op2_,st_]:=act[op1,act[op2,st]];
act[op1_+op2_,st_]:=act[op1,st]+act[op2,st];
act[op_,st1_+st2_]:=act[op,st1]+act[op,st2];
act[_,0.]:=0;
act[0.,_]:=0;
n[i_]:=b†[i]@*b[i];
rep[op_]:=SparseArray[Flatten[Function[{i,dat},{i,#1}->#2&@@@dat]@@@Thread[{Values@id,Block[{expr,kets,cs},expr=Expand[act[op,#]];kets=DeleteDuplicates[Cases[expr,_Ket,{0,∞}]];If[Length[kets]==0,{},cs=Flatten@Normal@Rest@CoefficientArrays[expr,kets];
Thread@{id/@kets,cs}]]&/@basis}]],{#,#}&@Length@basis];
Ht=rep[Sum[(b†[Mod[i+1,Ν,1]]@*b[i]+b†[i]@*b[Mod[i+1,Ν,1]]),{i,Ν}]];
HU=rep[Sum[n[i]@*(n[i]-𝟙),{i,Ν}]];
bb=rep[b†[Ceiling[Ν/2]]@*b[1]];
Block[{gs,Uts,Cs},
gs=Range[0,1,0.05];
Cs=Table[Block[{t=0.17(1-g),U=g,H,E0,Ψ0},
H=-t Ht+(U/2)HU;
{{E0},{Ψ0}}=Eigensystem[N@H,1,Method->{"Arnoldi","Shift"->-2Νb-0.1}];
Conjugate[Ψ0].bb.Ψ0],{g,gs}];
ListLinePlot[Thread@{gs,Cs},PlotRange->{0,1},PlotRangePadding->0.05,GridLines->{{1/2},{0,1}},PlotStyle->Black,PlotLabel->StringTemplate["\!\(N=`1`, N\_b=`2`\)"][Ν,Νb],FrameTicks->{{Automatic,Automatic},{{#,If[#==1,∞,NumberForm[#/(0.17(1-#)),2]]}&/@Range[0,1,0.2],None}},FrameLabel->{"\!\(U/t\)","\!\(G(d=3)\)"},ImageSize->180]
]]

Weak interacting phase: boson correlation G(d) decays in a power law (algebraically) with the distance d

G(d)∼(d^(-α)).

Strong interacting phase: boson correlation G(d) decays exponentially with the distance d

G(d)∼(e^(-d/ξ)).

Our current exact diagonalization approach will not be able to further study these behaviors and pin down the KT transition, because the system size we can achieve is too small. Yet, it still provides us qualitative understanding of helpful physical pictures.

Sachdev-Ye-Kitaev Model

Majorana Fermion

Real and Imaginary Decomposition: Majorana fermions can be thought as the “real part” and “imaginary part” of the traditional fermion operator

({{{, ({{(Subscript[c, i])=(1/2)((Subscript[χ, 2i-1])+i (Subscript[χ, 2i]))}, {(Subscript[c, i]^†)=(1/2)((Subscript[χ, 2i-1])-i (Subscript[χ, 2i]))}})}})  ({{{, ({{(Subscript[χ, 2i-1])=(Subscript[c, i]^†)+(Subscript[c, i])}, {(Subscript[χ, 2i])=i ((Subscript[c, i]^†)-(Subscript[c, i]))}})}})

These Majorana fermion operators still acts on the states in the same fermion many-body Hilbert space. They are just linear combination of (Subscript[c, i]) and (Subscript[c, i]^†).

Anti-commutation Relations: Majorana fermion operators (Subscript[χ, i]) are defined by the following algebraic property:

This is also the defining property of the Clifford algebra generators.

Jordan-Wigner Transformation: Majorana fermion operator can be systematically represented as Pauli matrices

Model Hamiltonian

(Subscript[χ, i]) are Majorana fermion operators,

(Subscript[J, i j k l]) are random coupling constants drawn from independent identical Gaussian distributions, with

(OverBar[(Subscript[J, i j k l])])=0, (OverBar[(Subscript[J, i j k l]^2)])=(2/(N^3)).

(Subscript[ℤ, 2]^F) Symmetry

(Subscript[n, i])=(Subscript[c, i]^†)(Subscript[c, i]) is the fermion number operator,

(∑_i)(Subscript[n, i]) counts the total number of fermions in the system,

(((-))^F)=+1 on those states with even number of fermions and -1 on those states with odd number of fermions, hence the name fermion parity.

Fermion parity operator is a two-fold operator (like an Ising symmetry)

(((-))^F)(((-))^F)=𝟙.

therefore, the symmetry group is (Subscript[ℤ, 2]). (A superscript F is usually put on (Subscript[ℤ, 2]^F) to indicate the Fermion parity symmetry — a fundamental symmetry unique to fermion systems).

Under fermion parity transformation: all fermion operators changes sign

(Subscript[c, i])->(((-))^F)(Subscript[c, i]) (((-))^F)=-(Subscript[c, i]),
(Subscript[c, i]^†)->(((-))^F)(Subscript[c, i]^†) (((-))^F)=-(Subscript[c, i]^†),
(Subscript[χ, i])->(((-))^F)(Subscript[χ, i]) (((-))^F)=-(Subscript[χ, i]).

But since each term in the Hamiltonian always contains even number of fermion operators, the fermion parity symmetry is always respected.

Prove Eq. () given the definition Eq. ().

(Subscript[c, i])(Subscript[n, i])=((Subscript[n, i])+𝟙)(Subscript[c, i]),

lhs=(Subscript[c, i])(Subscript[n, i])=(Subscript[c, i])(Subscript[c, i]^†)(Subscript[c, i])=(𝟙-(Subscript[c, i]^†)(Subscript[c, i]))(Subscript[c, i])=(Subscript[c, i]),
rhs=((Subscript[n, i])+𝟙)(Subscript[c, i])=((Subscript[c, i]^†)(Subscript[c, i])+𝟙)(Subscript[c, i])=(Subscript[c, i]).

(Subscript[c, i]) (((-))^F)=(Subscript[c, i]) (((-1))^((∑_i)(Subscript[n, i])))=(((-1))^((∑_i)(Subscript[n, i])+𝟙))(Subscript[c, i])=-(((-1))^((∑_i)(Subscript[n, i])))(Subscript[c, i])=-(((-))^F)(Subscript[c, i]),

(((-))^F)(Subscript[c, i]) (((-))^F)=-(((-))^F)(((-))^F)(Subscript[c, i])=-(Subscript[c, i]),

Fermion parity conservation: as a result, the Hamiltonian commutes with the fermion parity operator

such that the fermion parity is always conserved in fermion systems.

Jordan-Wigner Transformation: the fermion parity operator can be represented as

It indeed anti-commutes with all fermion operators, as expected.

Matrix Representations

First we can construct the matrix representation for Majorana fermion operators using Jordan-Wigner transformation Eq. ().

Block[{Ν=8,χs},
χs=Block[{n=Ceiling[Ν/2]},
SparseArray[Fold[KroneckerProduct,SparseArray/@PauliMatrix/@#]&/@Take[Flatten[Table[PadRight[Append[ConstantArray[3,m],a],n],{m,0,n-1},{a,{1,2}}],1],Ν]]];
Row[ComplexMatrixPlot[#1,PlotLabel->StringTemplate["\!\(χ\_`1`\)"][#2],FrameStyle->Black,FrameTicks->False,PlotRangePadding->0,Mesh->All,MeshStyle->Directive[AbsoluteThickness[1],Opacity[0.1]],ImageSize->60]&@@@Thread[{χs,Range@Ν}]]]

👉 Legend: each pixel corresponds to a matrix element. 
Color code: =1, =i, =-1, =-i, =0.

Then we generate random coupling constant tensor (Subscript[J, i j k l]) by random sampling, and construct the SYK Hamiltonian H by tensor contraction.

DynamicModule[{Ν=8,χs,FP,getH,H},
χs=Block[{n=Ceiling[Ν/2]},
SparseArray[Fold[KroneckerProduct,SparseArray/@PauliMatrix/@#]&/@Take[Flatten[Table[PadRight[Append[ConstantArray[3,m],a],n],{m,0,n-1},{a,{1,2}}],1],Ν]]];
FP=Fold[KroneckerProduct,SparseArray/@PauliMatrix/@ConstantArray[3,Ceiling[Ν/2]]];
getH[]:=Block[{J},
J=Sqrt[2/Ν^3]Symmetrize[RandomVariate[NormalDistribution[],{Ν,Ν,Ν,Ν}],Antisymmetric[{1,2,3,4}]];
Activate[TensorContract[Inactive[TensorProduct]@@{J,χs,χs,χs,χs},{{1,5},{2,8},{3,11},{4,14},{7,9},{10,12},{13,15}}]]];
H=getH[];
EventHandler[Dynamic[Deploy[Row[ComplexMatrixPlot[#1,PlotLabel->#2,FrameStyle->Black,FrameTicks->False,PlotRangePadding->0,Mesh->All,MeshStyle->Directive[AbsoluteThickness[1],Opacity[0.1]],ImageSize->120]&@@@{{H,"\!\( H\)"},{FP,"\!\((-)\^F\)"}},Spacer[1]]],TrackedSymbols:>{H}],{"MouseClicked":>(H=getH[])}]]

Some matrix elements in H are always 0, why?

Inspect the fermion parity operator (((-))^F). It becomes clear that those matrix elements in H are forbidden by the (Subscript[ℤ, 2]^F) symmetry, as they correspond to transitions between states of different fermion parities (thereby violating the fermion parity conservation).

If we rearrange the basis state by fermion parity into even and odd parity subspaces, the Hamiltonian will be transformed to the block diagonal form.

DynamicModule[{Ν=8,χs,FP,getH,ord,H},
χs=Block[{n=Ceiling[Ν/2]},
SparseArray[Fold[KroneckerProduct,SparseArray/@PauliMatrix/@#]&/@Take[Flatten[Table[PadRight[Append[ConstantArray[3,m],a],n],{m,0,n-1},{a,{1,2}}],1],Ν]]];
FP=Fold[KroneckerProduct,SparseArray/@PauliMatrix/@ConstantArray[3,Ceiling[Ν/2]]];
ord=Ordering@Normal@Diagonal@FP;
FP=FP[[ord,ord]];
getH[]:=Block[{J},
J=Sqrt[2/Ν^3]Symmetrize[RandomVariate[NormalDistribution[],{Ν,Ν,Ν,Ν}],Antisymmetric[{1,2,3,4}]];
Activate[TensorContract[Inactive[TensorProduct]@@{J,χs,χs,χs,χs},{{1,5},{2,8},{3,11},{4,14},{7,9},{10,12},{13,15}}]][[ord,ord]]];
H=getH[];
EventHandler[Dynamic[Deploy[Row[ComplexMatrixPlot[#1,PlotLabel->#2,FrameStyle->Black,FrameTicks->False,PlotRangePadding->0,Mesh->All,MeshStyle->Directive[AbsoluteThickness[1],Opacity[0.1]],ImageSize->120]&@@@{{H,"\!\( H\)"},{FP,"\!\((-)\^F\)"}},Spacer[1]]],TrackedSymbols:>{H}],{"MouseClicked":>(H=getH[])}]]

The block diagonal structure reduces the computational complexity for exact diagonalization, as it enables us to diagonalize each block independently. — The complexity is reduced to 1/8+1/8=1/4 of the original task.

Quantum Chaos

Examples of Classical Chaos:

Weather Systems: long-term weather forecast is almost impossible.

Double Pendulum: a simple pendulum is periodic, but you can’t predict the behavior of a pendulum hanging below another one.

Three-Body Problem (in classical gravity): there is a famous Chinese science fiction [by Cixin Liu] about it ...

Key Takeaway: unpredictable and random behavior can emerge in deterministic systems due to sensitivity to initial conditions.

Naively, quantum system should not exhibit chaos: because the time-evolution equation is linear, small changes in the initial state remain small.

Instead, quantum chaos manifests through other features, such as

Energy-level statistics, [][]

Eigenstate thermalization hypothesis (ETH) and volume-law entanglement entropy,

Local operator growth and quantum information scrambling.

E. Wigner, Characteristic Vectors of Bordered Matrices
with Infinite Dimensions, Ann. Math. 62, 548 (1955).

F. J. Dyson, Statistical Theory of the Energy Levels of
Complex Systems. I, J. Math. Phys. (N.Y.) 3, 140 (1962).

Energy-Level Statistics (Billiard)

Rectangle shape: the motion will be very regular (integrable). (the horizontal and vertical motions are decoupled, and both are periodic)

Stadium shape: the motion will be chaotic.

Block[{Ω1,Ω2},
Ω1=Rectangle[{0,0},{5^(1/4),1}];
Ω2=DiscretizeRegion[RegionUnion[{Disk[{0,1/2},1/2],Rectangle[{0,0},{1/7^(1/3),1}],Disk[{1/7^(1/3),1/2},1/2]}]];
Row[RegionPlot[#1,PlotLabel->#2,BoundaryStyle->Black,PlotStyle->LightGray,Frame->False,AspectRatio->Automatic,ImageSize->{Automatic,80}]&@@@{{Ω1,"Rectangle"},{Ω2,"Stadium"}},Spacer[1]]
]

(Subscript[ψ, n])(x) is the eigen wave function, restricted in x∈Ω (with the Dirichlet boundary condition (Subscript[ψ, n])(x∈∂Ω)=0),

(Subscript[E, n]) is the corresponding eigen energy (energy level).

(Subscript[E, 0])<=(Subscript[E, 1])<=(Subscript[E, 2])<=…<=(Subscript[E, n])≤(Subscript[E, n+1])≤…,

the energy level spacing

the log level spacing ratio (to construct a dimensionless scale)

Block[{Ω1,Ω2,En1,En2,sn1,sn2,λn1,λn2,pWD,PWD},
Ω1=Rectangle[{0,0},{5^(1/4),1}];
Ω2=DiscretizeRegion[RegionUnion[{Disk[{0,1/2},1/2],Rectangle[{0,0},{1/7^(1/3),1}],Disk[{1/7^(1/3),1/2},1/2]}]];
{En1,En2}=NDEigenvalues[{-Laplacian[u[x,y],{x,y}],DirichletCondition[u[x,y]==0,True]},u[x,y],{x,y}∈#,1000,Method->{"PDEDiscretization"->{"FiniteElement",{"MeshOptions"->{"MaxCellMeasure"->0.002}}}}]&/@{Ω1,Ω2};
{sn1,sn2}=Differences/@{En1,En2};
{λn1,λn2}=Log[Rest[#]/Most[#]]&/@{sn1,sn2};
pWD[β_][s_]:=If[β==0,Exp[-s],Block[{B,A},
B=Switch[β,1,π/4,2,4/π,4,64/9π];
A=2B^(β/2+1)Gamma[β/2+1];
A s^β Exp[-B s^2]]];
PWD[β_][λ_]:=If[β==0,1/(2 Cosh[λ/2])^2,Block[{Z,r},
Z=Switch[β,1,8/27,2,4π/(81 Sqrt[3]),4,4π/(729Sqrt[3])];
r=Exp[λ];
r(r+r^2)^β/(1+r+r^2)^(1+3β/2)/Z]];
Grid[{Style[#,"Text"]&/@{"Rectangle","Stadium"},Histogram[#1/Mean[#1],{0.1},"PDF",Epilog->{Line[Table[{s,#2[s]},{s,0,3,0.1}]]},ChartStyle->LightGray,PlotRange->{{0,3},{0,1}},PlotRangeClipping->True,FrameLabel->{"\!\( s/s\&_\)","\!\(p(s)\)"},ImageSize->180]&@@@{{sn1,pWD[0]},{sn2,pWD[1]}},
Histogram[#1,{0.2},"PDF",Epilog->{Line[Table[{s,#2[s]},{s,-3,3,0.1}]]},ChartStyle->LightGray,PlotRange->{{-3,3},{0,0.5}},PlotRangeClipping->True,FrameLabel->{"\!\( λ\)","\!\(p(λ)\)"},ImageSize->180]&@@@{{λn1,PWD[0]},{λn2,PWD[1]}}}]
]

Rectangle case: the classical dynamics is integrable, the energy levels will be uncorrelated, and the level statistics p(s) follows Poisson distribution, correspondingly the distribution p(λ) will be boarder.

What is Poisson distribution? Sample N independent, uniformly distributed random real numbers from [0,N], sort them, and calculate the spacings s between consecutive values. In the large-N limit, the spacings s will follow the Poisson distribution

p(s)=(e^(-s)).

Stadium case: the classical dynamics is chaotic, the energy levels are correlated with level repulsion, and the level statistics p(s) follows Wigner-Dyson distribution (more specifically the GOE distribution here), correspondingly the distribution p(λ) will be more concentrated.

What are Wigner-Dyson distributions? Sample independent Gaussian random numbers to construct N×N Hermitian matrices, and calculate the spacings s between consecutive eigenvalues. In the large-N limit, the spacing s will follow the Wigner-Dyson distribution

p(s)∝(s^β)(e^(-α (s^2))),

where β depends on whether the matrix elements are:

real - Gaussian Orthogonal Ensemble (GOE): β=1,

complex - Gaussian Unitary Ensemble (GUE): β=2,

quaternion - Gaussian Symplectic Ensemble (GSE): β=4.

Energy-Level Statistics (SYK)

Legends:({{({{}, { }}), (Underscript[({{, , }, {GOE, GUE, GSE}}), ⏟])}, {Poisson, Wigner-Dyson}})

Grid[Partition[Block[{Ν=#1,M,χs,FP,ie,io,λs,PWD,col},
M=2^(13-Ceiling[Ν/2]);
χs=Block[{n=Ceiling[Ν/2]},
SparseArray[Fold[KroneckerProduct,SparseArray/@PauliMatrix/@#]&/@Take[Flatten[Table[PadRight[Append[ConstantArray[3,m],a],n],{m,0,n-1},{a,{1,2}}],1],Ν]]];
FP=Fold[KroneckerProduct,SparseArray/@PauliMatrix/@ConstantArray[3,Ceiling[Ν/2]]];
{ie,io}=PositionIndex[Normal@Diagonal@FP]/@{1,-1};
λs=Flatten[Table[Block[{J,Hfull},
J=Sqrt[2/Ν^3]Symmetrize[RandomVariate[NormalDistribution[],{Ν,Ν,Ν,Ν}],Antisymmetric[{1,2,3,4}]];
Hfull=Activate[TensorContract[Inactive[TensorProduct]@@{J,χs,χs,χs,χs},{{1,5},{2,8},{3,11},{4,14},{7,9},{10,12},{13,15}}]];
Flatten[Log[Rest[#]/Most[#]]&@DeleteCases[Differences@Sort@Eigenvalues[#],x_/;Abs[x]<1*^-10]&/@{Hfull[[ie,ie]],Hfull[[io,io]]}]],M]];
PWD[β_][λ_]:=If[β==0,1/(2 Cosh[λ/2])^2,Block[{Z,r},
Z=Switch[β,1,8/27,2,4π/(81 Sqrt[3]),4,4π/(729Sqrt[3])];
r=Exp[λ];
r(r+r^2)^β/(1+r+r^2)^(1+3β/2)/Z]];
col=<|0->Blue,1->Red,2->Orange,4->Brown|>;
Histogram[λs,{0.2},"PDF",Epilog->{{Opacity[0.5],Table[{col[β],Line[Table[{s,PWD[β][s]},{s,-3,3,0.1}]]},{β,{0,1,2,4}}]},Text[StringTemplate["\!\(N=`1`\)"][Ν],Scaled@{0,1},{-1.2,1.2}],{col[#2],Text[<|1->"GOE",2->"GUE",4->"GSE"|>[#2],Scaled@{1,1},{1.2,1.2}]}},ChartStyle->LightGray,PlotRange->{{-3,3},{0,0.8}},PlotRangeClipping->True,FrameLabel->{"\!\( λ\)","\!\(p(λ)\)"},ImageSize->180]]&@@@{{8,1},{9,1},{10,2},{11,4},{12,4},{13,4},{14,2},{15,1}},2]]

For various different Majorana mode number N, the SYK model always exhibits the Wigner-Dyson level statistics, indicating quantum chaotic behavior.

Interestingly, the distribution cycles through the three Gaussian ensembles (GOE, GUE, GSE) with respect to N with a period 8, a phenomenon that is profoundly related to the quantum anomaly associated with the (Subscript[ℤ, 2]^F) and time-reversal symmetry of the SYK model. []

Yi-Zhuang You, Andreas W. W. Ludwig, Cenke Xu. Sachdev-Ye-Kitaev Model and Thermalization on the Boundary of Many-Body Localized Fermionic Symmetry Protected Topological States. arXiv:1602.06964