Condensed Matter Physics

Part 2. Electronic System

Electrons in Solids

Quantum Mechanics of Electrons

Schrödinger Equation of Electron

(-(((ℏ^2)(∇^2))/(2m))+V(x))ψ(x)=E ψ(x).

m - electron mass,

V(x) - electronic potential at position x,

ψ(x) - electron wave function (eigenstate), s.t. the probability (density) to observe an electron at position x is ((|ψ(x)|)^2),

E - the corresponding eigen energy of the eigenstate.

A Toy Atom in 1D

Block[{γ=0.2},Row[{Graphics[{Arrowheads[0.5],AbsoluteThickness[1],Arrow[{0,#}&/@{-3,1}],Line[Offset[{#,0}]&/@{-2,2}],Text["\!\(V(x)\)",{0,1},{1.2,1}],Text["0",{0,0},{2,0.5}]},AspectRatio->6,ImageSize->{Automatic,80}],Row[Append[Table[Graphics[{{LightGray,Polygon[{{-1,0},{-a,0},{-a,-γ/a},{a,-γ/a},{a,0},{1,0},{1,-3},{-1,-3}}]},{Line[{{-1,0},{-a,0},{-a,-γ/a},{a,-γ/a},{a,0},{1,0}}]}},PlotRange->{{-1,1},{-3,1}},ImageSize->{Automatic,80}],{a,{0.4,0.2,0.1}}],Graphics[{{LightGray,Polygon[{{-1,0},{1,0},{1,-3},{-1,-3}}]},{Line[{{-1,0},{1,0}}]},{Arrowheads[0.3],Arrow[{0,#}&/@{0,-2.8}]}},PlotRange->{{-1,1},{-3,1}},ImageSize->{Automatic,80}]],Style["→",FontFamily->"LaTeX"]]}]]

-((ℏ^2)/(2m))(∂_x^2)ψ(x)-γ δ(x)ψ(x)=E ψ(x).

Bound state (E<0)

Wave function decays exponentially near the origin,

with a decay length

Graphics[{{Arrowheads[0.05],AbsoluteThickness[1],Arrow[{#,0}&/@{-3,3}],Arrow[{0,#}&/@{-0.5,1.5}],Text["\!\(,x\)",{3,0},{-1,0}],Text["\!\(ψ(x)\)",{0,1.5},{-1.5,0}]},{AbsoluteThickness[1],Dashed,Line[{{-1,0},{0.2,1.2}}],Line[{{1,0},{-0.2,1.2}}],Text["\!\(0\)",{0,0},{1,1}],Text["\!\(ξ\)",{1,0},{0,1}],Text["\!\(-ξ\)",{-1,0},{0,1}]},
First@Plot[Exp[-Abs[x]],{x,-2.5,2.5},PlotStyle->Black]}]

In the bound state, the electron is trapped around the potential within the length scale ξ. => The model mimics an atom (in 1D), where ξ can be viewed as the atomic radius (Bohr radius).

The eigen energy is negative (with respect to the potential zero)

This is the unique bound state in the system, which is also the lowest-energy ground state of the electron.

Verify that Eq. () is a normalized solution of Eq. ().

(∂_x)ψ(x)=(ξ^(-1/2))(∂_x)(e^(-|x|/ξ))
=(ξ^(-1/2))(-(ξ^(-1)))(e^(-|x|/ξ))(∂_x)|x|
=-(ξ^(-3/2))(e^(-|x|/ξ))sgn(x),

(∂_x^2)ψ(x)=-(ξ^(-3/2))(∂_x)((e^(-|x|/ξ))sgn(x))
=-(ξ^(-3/2))(((∂_x)(e^(-|x|/ξ)))sgn(x)+(e^(-|x|/ξ))((∂_x)sgn(x)))
=-(ξ^(-3/2))((-(ξ^(-1)))(e^(-|x|/ξ))((sgn(x))^2)+2(e^(-|x|/ξ))δ(x))
=((ξ^(-2))-2(ξ^(-1))δ(x))(ξ^(-1/2))(e^(-|x|/ξ))
=((ξ^(-2))-2(ξ^(-1))δ(x))ψ(x).

-((ℏ^2)/(2m))(∂_x^2)ψ(x)-γ δ(x)ψ(x)
=-((ℏ^2)/(2m))((ξ^(-2))-2(ξ^(-1))δ(x))ψ(x)-γ δ(x)ψ(x)
=-((ℏ^2)/(2m))((1/(ξ^2))-2(((m γ)/ℏ)-(1/ξ))δ(x))ψ(x).

-((ℏ^2)/(2m))((1/(ξ^2))-2(((m γ)/ℏ)-(1/ξ))δ(x))ψ(x)=E ψ(x),

((m γ)/ℏ)-(1/ξ)=0, -((ℏ^2)/(2m (ξ^2)))=E.

ξ=((ℏ^2)/(m γ)), E=-((ℏ^2)/(2m (ξ^2)))=-((m (γ^2))/(2(ℏ^2))).

(∫_(-∞)^∞)((|ψ(x)|)^2)dx=(∫_(-∞)^∞)(ξ^(-1))(e^(-2|x|/ξ))dx
=(2/ξ)(∫_0^∞)(e^(-2x /ξ))dx
=(∫_0^∞)(e^(-2x /ξ))d(2x/ξ)
=1.

Scattering states (E>0)

There are infinitely many scattering state solutions (forming a continuum in the spectrum). Each scattering state is associated with an positive eigen energy E>0, which also sets a wave number

k=((Sqrt[2m E])/ℏ).

Wave function

ψ(x)=({{{, ({{(e^(i k x))+(Subscript[A, r]) (e^(-i k x)), x<0}, {(Subscript[A, t]) (e^(i k x)), x>0}})}}),

with the reflection and transmission amplitudes

(Subscript[A, r])=-(1/(1+i k ξ)),
(Subscript[A, t])=(1/(1-i /(k ξ))).

Manipulate[Plot[Block[{Ar,At},{Ar,At}={-1/(1+I k),1/(1-I/k)};Re[Exp[-I t]Piecewise[{{Exp[I k x]+Ar Exp[-I k x],x<0},{At Exp[I k x],x>0}}]]],{x,-10π,10π},PlotStyle->Black,PlotRange->{-2,2},FrameLabel->{"\!\(,x\)","\!\(ψ(x)\)"},ImageSize->180],{{k,1},0.01,2,ImageSize->Small},{t,0,2π,ImageSize->Small},LabelStyle->{"InlineFormula",12},Paneled->False,SaveDefinitions->True]

In the scattering state, the electron has high enough energy to scatter off the potential. => The atom ionize into dissociated electron and nucleus.

In the following, we will be interested in bound states of electrons (in atomic orbitals), which are relevant to the solid state physics.

Two Toy Atoms in 1D

Bound state (E<0) wave functions

Even parity solution

ψ(x)=({{{, ({{A (e^(κ x)), x<-d/2,}, {cosh κ x, -d/2<x<d/2,}, {A (e^(-κ x)), x>d/2.}})}})

Odd parity solution

ψ(x)=({{{, ({{-A (e^(κ x)), x<-d/2,}, {sinh κ x, -d/2<x<d/2,}, {A (e^(-κ x)), x>d/2.}})}})

The amplitude A is given by (continuity of wave function)

A=({{{, ({{(1/2)((e^(κ d))+1), even parity,}, {(1/2)((e^(κ d))-1), odd parity.}})}})

Substitute Eq. (), Eq. () into the Schrödinger equation, the imaginary wave number κ is determined by

(e^(κ d))(κ ξ-1)=({{{, ({{+1, even parity,}, {-1, odd parity,}})}})

Derive Eq. ().

-((ℏ^2)/(2m))(∂_x^2)ψ(x)-γ δ(x-d/2)ψ(x)-γ δ(x+d/2)ψ(x)=E ψ(x).

(∂_x^2)ψ(x)=(κ^2)ψ(x),

(∂_x^2)ψ(±d/2)=((∂_x)ψ(±d/2+(SubPlus[0]))-(∂_x)ψ(±d/2+(SubMinus[0])))δ(x∓d/2).

For even parity solution

(∂_x)ψ(d/2+(SubPlus[0]))-(∂_x)ψ(d/2+(SubMinus[0]))
=(Subscript[((∂_x)A (e^(-κ x))-(∂_x)cosh κ x), x=d/2])
=-(1/2)((e^(κ d))+1) (e^(-κ d/2))-κ sinh(κ d/2)
=-κ (e^(κ d/2)),

(∂_x)ψ(-d/2+(SubPlus[0]))-(∂_x)ψ(-d/2+(SubMinus[0]))
=(Subscript[((∂_x)cosh κ x-(∂_x)A (e^(κ x))), x=-d/2])
=κ sinh(-κ d/2)-(κ/2)((e^(κ d))+1) (e^(-κ d/2))
=-κ (e^(κ d/2)).

For odd parity solution

(∂_x)ψ(d/2+(SubPlus[0]))-(∂_x)ψ(d/2+(SubMinus[0]))
=(Subscript[((∂_x)A (e^(-κ x))-(∂_x)sinh κ x), x=d/2])
=-(κ/2)((e^(κ d))-1) (e^(-κ d/2))-κ cosh(κ d/2)
=-κ (e^(κ d/2)),

(∂_x)ψ(-d/2+(SubPlus[0]))-(∂_x)ψ(-d/2+(SubMinus[0]))
=(Subscript[((∂_x)sinh κ x-(∂_x)(-A (e^(κ x)))), x=-d/2])
=κ cosh(-κ d/2)+(κ/2)((e^(κ d))-1)(e^(-κ d/2))
=κ (e^(κ d/2)),

(∂_x^2)ψ(±d/2)=-κ (e^(κ d/2))δ(x∓d/2)(((±1))^p),

-((ℏ^2)/(2m))((κ^2)ψ(x)-κ (e^(κ d/2))δ(x-d/2)-(((-1))^p)κ (e^(κ d/2))δ(x+d/2))-γ δ(x-d/2)ψ(x)-γ δ(x+d/2)ψ(x)=E ψ(x).

-((ℏ^2)/(2m))(κ^2)=E,

-((ℏ^2)/(2m))(-κ (e^(κ d/2)))δ(x-d/2)-γ δ(x-d/2)ψ(x)=0,

-((ℏ^2)/(2m))(-(((-1))^p)κ (e^(κ d/2)))δ(x+d/2)-γ δ(x+d/2)ψ(x)=0,

(((ℏ^2)κ (e^(κ d/2)))/(2m)) δ(x-d/2)
=γ δ(x-d/2)ψ(x)
=γ ψ(d/2)δ(x-d/2)
=γ A (e^(-κ d/2))δ(x-d/2),

(((-1))^p)(((ℏ^2)κ (e^(κ d/2)))/(2m))δ(x+d/2)
=γ δ(x+d/2)ψ(x)
=γ ψ(-d/2)δ(x+d/2)
=(((-1))^p)γ A (e^(-κ d/2))δ(x+d/2).

(((ℏ^2)κ (e^(κ d/2)))/(2m))=γ A (e^(-κ d/2))=(γ/2)((e^(κ d))+(((-1))^p))(e^(-κ d/2)),

(e^(κ d))((((ℏ^2)κ)/(m γ))-1)=(((-1))^p).

The number of solutions of Eq. () depends on the inter-atomic distance d and the atomic radius ξ,

d>ξ: two bound states of both parities,

d<ξ: one bound state of even parity.

Manipulate[Block[{L=5},
Graphics[{{AbsoluteThickness[1],Arrowheads[0.05],
Arrow[{#,0}&/@{-L,L}],Line[Offset[{0,#}]&/@{-2,2}],Arrow[{-L,#}&/@{-4,0.5}]},Translate[{{LightGray,Line[{0,#}&/@{-4,0.5}]},{White,EdgeForm[Black],Disk[{0,0},Offset[{5,5}]]},{Text["+"]}},{# d/2,0}]&/@{1,-1},
{Text["\!\(,x\)",{L,0},{0,-1}],Text["\!\(,E\)",{-L,0.5},{-2,0.8}]},
Table[Block[{col,κ,A,Ε},
col=Switch[p,0,Blue,1,Red];
If[d<=1&&p==1,
{},
κ=If[d==0,2,1+ProductLog[(-1)^p d Exp[-d]]/d];
Ε=-κ^2;
A=(Exp[κ d]+(-1)^p)/2;
Translate[{{AbsoluteThickness[1],Line[{#,0}&/@{-L,L}]},First@Plot[Piecewise[{{(Exp[κ x]+(-1)^p Exp[-κ x])/2,Abs[x]<d/2},{A Exp[-κ x],x>d/2},{(-1)^p  A Exp[κ x],x<-d/2}}]/(2A Exp[-κ d/2]/κ),{x,-L,L},PlotRange->All,PlotStyle->col,Filling->Axis,Exclusions->None]},{0,Ε}]]],{p,{0,1}}]},AspectRatio->0.8,PlotRange->{-4,0.5},PlotRangePadding->Scaled[.02],PlotLabel->Style["\!\(d/ξ="<>ToString[d]<>"\)",14,Black],ImageSize->180]],{{d,1},0,5,ImageSize->Small},Paneled->False,LabelStyle->"InlineFormula"]

: even parity solution,
: odd parity solution.

Based on the solutions of κ, the eigen energies are given by

Plot[{-If[d>0,1+ProductLog[d Exp[-d]]/d,2]^2,-If[d>1,1+ProductLog[-d Exp[-d]]/d,Indeterminate]^2},{d,0,5},GridLines->{None,{-1}},FrameLabel->{"\!\(d/ξ\)","\!\(E/\[LeftBracketingBar]E\_0\[RightBracketingBar]\)"},PlotRange->{-4,0},ImageSize->180]

The energy E is plotted in unit of the ground state energy (Subscript[E, 0]) of a single atom given in Eq. ()

(Subscript[E, 0])=-((ℏ^2)/(2m (ξ^2))).

The even parity state is always the ground state of the two-atom system.

Assuming electrons always occupy the ground state, as the atoms merge together (from d→∞ to d->0), the binding energy (Subscript[E, b])=(Subscript[E, d→∞])-(Subscript[E, d=0])=3|(Subscript[E, 0])| will be released (per electron).

Force between atoms are given by

F=-((∂E)/(∂d))

Inter-atomic force can be mediated by exchanging electron => valence bond.

Molecular orbitals of valence σ-bond

Block[{orbit},
orbit[col_,dat_]:=Block[{L=2},
{{AbsoluteThickness[1],Line[{#,0}&/@{-L,L}]},
First@Plot[Total[#1 Exp[-2Abs[x-#2]]&@@@dat],{x,-L,L},PlotStyle->col,PlotRange->{-1,1},Filling->Axis],
Translate[{{White,EdgeForm[Black],Disk[{0,0},Offset[{4,4}]]},{Text["+"]}},{#2,0}]&@@@dat}];
Graphics[{
Translate[{orbit[Red,{{-1/Sqrt[2],-1},{1/Sqrt[2],1}}],Text["↑↑",{0,0.5}],Text["\!\(σ\^*\)",{0,-0.8}]},{0,3}],
Translate[{orbit[Blue,{{1/Sqrt[2],-1},{1/Sqrt[2],1}}],Text["↑↓",{0,0.5}],Text["\!\(σ\)",{0,-0.8}]},{0,-3}],
Translate[{orbit[Black,{{1,0}}],Text["↑",{-0.5,0.8}]},{-5,0}],
Translate[{orbit[Black,{{1,0}}],Text["↑",{0.5,0.8}]},{5,0}],
{Dashed,Scale[Line[{{-3,0},{-2,3}}],#,{0,0}]&/@{{1,1},{1,-1},{-1,1},{-1,-1}}},
Translate[{Arrow[{0,#}&/@{-4,4}],Text["\!\(,E\)",{0,3},{-2,0}]},{-7.5,0}]},ImageSize->200]]

Bonding orbital: electron in even parity state, atoms attract each other.

Anti-bonding orbital: electron in odd parity state, atoms repel each other.

Tight Binding Approximation

Atomic orbitals

where (Subscript[x, 1]) and (Subscript[x, 2]) are the positions of two nuclei respectively.

(Ket[1]) and (Ket[2]) are not orthogonal, but close to orthogonal in the tight binding limit

(BraKet[1, 2])=(∫_(-∞)^∞)(Subscript[ψ, 1]^( * ))(x)(Subscript[ψ, 2])(x)dx
=(ξ^(-1))(∫_(-∞)^∞)(e^(-|x-(Subscript[x, 1])|/ξ))(e^(-|x-(Subscript[x, 2])|/ξ))dx
=(e^(-d/ξ))(1+d/ξ)(Overscript[→, d≫ξ])0,

where d=|(Subscript[x, 1])-(Subscript[x, 2])|.

Tight binding basis: orthogonalized atomic orbitals

such that (BraKet[(OverTilde[i]), (OverTilde[j])])=(Subscript[δ, i,j]).

The transformation matrix is given by

S=(({{(Subscript[S, (OverTilde[1])1]), (Subscript[S, (OverTilde[1])2])}, {(Subscript[S, (OverTilde[2])1]), (Subscript[S, (OverTilde[2])2])}}))=(Sqrt[cosh 2ε])(({{cosh ε, -sinh ε}, {-sinh ε, cosh ε}})).

where ε is defined by

tanh 2ε=(e^(-d/ξ))(1+d/ξ).

Block[{ε,S},
S=Sqrt[Cosh[2ε]]{{Cosh[ε],-Sinh[ε]},{-Sinh[ε],Cosh[ε]}};
Simplify[S.{{1,Tanh[2ε]},{Tanh[2ε],1}}.Transpose[S]]]

Verification.

ε is expected to vanish in the tight biding limit. To the leading order in ε,

(Ket[(OverTilde[1])])≃(Ket[1])-ε(Ket[2]),
(Ket[(OverTilde[2])])≃(Ket[2])-ε(Ket[1]),

with ε≃(1/2)(e^(-d/ξ))(1+d/ξ).

Hamiltonian operator

H=-((ℏ^2)/(2m))(∂_x^2)+V(x)
=-((ℏ^2)/(2m))(∂_x^2)-γ δ(x-(Subscript[x, 1]))-γ δ(x-(Subscript[x, 2])).

It is assumed that γ and ξ are always related by ξ=(ℏ^2)/(m γ), such that (Ket[1]) and (Ket[2]) are eigenstates of the Hamiltonian when only one potential is present, i.e.

(-((ℏ^2)/(2m))(∂_x^2)-γ δ(x-(Subscript[x, i])))(Subscript[ψ, i])(x)=(Subscript[E, 0])(Subscript[ψ, i])(x),

with the atomic orbital eigen energy

Tight binding Hamiltonian: the effective Hamiltonian in the tight binding basis

(Subscript[H, eff])=(∑_(i,j=1,2))(Ket[(OverTilde[i])])(Subscript[(OverTilde[H]), i,j])(Bra[(OverTilde[j])]),

where (Subscript[(OverTilde[H]), i,j])≡(Bra[(OverTilde[i])])H(Ket[(OverTilde[j])]).

The matrix element (Subscript[(OverTilde[H]), i,j]) can be obtained from

(Subscript[H, i,j])≡(Bra[i])H(Ket[j])
=(∫_(-∞)^∞)(Subscript[ψ, i]^( * ))(x)(-((ℏ^2)/(2m))(∂_x^2)-γ δ(x-(Subscript[x, 1]))-γ δ(x-(Subscript[x, 2]))) (Subscript[ψ, j])(x)dx,

via basis transformation Eq. ()

(Subscript[(OverTilde[H]), i,j])=(∑_(i,j))(Subscript[S, (OverTilde[i])i])(Subscript[H, i,j])(Subscript[S, (OverTilde[j])j]),

or written in terms of matrix multiplication (OverTilde[H])=S H Sᵀ.

Substitute Eq. () into Eq. (), one finds

H=(({{(Subscript[H, 11]), (Subscript[H, 12])}, {(Subscript[H, 21]), (Subscript[H, 22])}}))=(Subscript[E, 0])(({{1+ 2(e^(-2d/ξ)), (e^(-d/ξ))(3+d/ξ)}, {(e^(-d/ξ))(3+d/ξ), 1+ 2(e^(-2d/ξ))}})),

where (Subscript[E, 0]) is the atomic orbital energy.

Derive Eq. ().

(Subscript[H, 11])=(ξ^(-1))(∫_(-∞)^∞)(e^(-|x-(Subscript[x, 1])|/ξ))(-((ℏ^2)/(2m))(∂_x^2)-γ δ(x-(Subscript[x, 1]))-γ δ(x-(Subscript[x, 2]))) (e^(-|x-(Subscript[x, 1])|/ξ))dx
=(ξ^(-1))(∫_(-∞)^∞)(e^(-|x-(Subscript[x, 1])|/ξ))(-((ℏ^2)/(2m))(∂_x^2)-γ δ(x-(Subscript[x, 1]))) (e^(-|x-(Subscript[x, 1])|/ξ))dx-γ (ξ^(-1))(∫_(-∞)^∞)(e^(-|x-(Subscript[x, 1])|/ξ)) δ(x-(Subscript[x, 2])) (e^(-|x-(Subscript[x, 1])|/ξ))dx.

(-((ℏ^2)/(2m))(∂_x^2)-γ δ(x-(Subscript[x, 1]))) (e^(-|x-(Subscript[x, 1])|/ξ))=(Subscript[E, 0])(e^(-|x-(Subscript[x, 1])|/ξ)),

(Subscript[H, 11])=(Subscript[E, 0])(ξ^(-1))(∫_(-∞)^∞)(e^(-2|x-(Subscript[x, 1])|/ξ))dx-γ (ξ^(-1))(∫_(-∞)^∞)δ(x-(Subscript[x, 2])) (e^(-2|x-(Subscript[x, 1])|/ξ))dx
=(Subscript[E, 0])-γ (ξ^(-1)) (e^(-2|(Subscript[x, 2])-(Subscript[x, 1])|/ξ))
=(Subscript[E, 0])(1+ 2(e^(-2d/ξ))),

(Subscript[H, 21])=(ξ^(-1))(∫_(-∞)^∞)(e^(-|x-(Subscript[x, 2])|/ξ))(-((ℏ^2)/(2m))(∂_x^2)-γ δ(x-(Subscript[x, 1]))-γ δ(x-(Subscript[x, 2]))) (e^(-|x-(Subscript[x, 1])|/ξ))dx
=(ξ^(-1))(∫_(-∞)^∞)(e^(-|x-(Subscript[x, 2])|/ξ))(-((ℏ^2)/(2m))(∂_x^2)-γ δ(x-(Subscript[x, 1]))) (e^(-|x-(Subscript[x, 1])|/ξ))dx-γ (ξ^(-1))(∫_(-∞)^∞)(e^(-|x-(Subscript[x, 2])|/ξ)) δ(x-(Subscript[x, 2])) (e^(-|x-(Subscript[x, 1])|/ξ))dx.

(Subscript[H, 21])=(Subscript[E, 0])(ξ^(-1))(∫_(-∞)^∞)(e^(-|x-(Subscript[x, 2])|/ξ))(e^(-|x-(Subscript[x, 1])|/ξ))dx-γ (ξ^(-1))(∫_(-∞)^∞)(e^(-|x-(Subscript[x, 2])|/ξ)) δ(x-(Subscript[x, 2])) (e^(-|x-(Subscript[x, 1])|/ξ))dx
=(Subscript[E, 0])(e^(-d/ξ))(1+d/ξ)-γ (ξ^(-1)) (e^(-d/ξ))
=(Subscript[E, 0])(e^(-d/ξ))(3+d/ξ).

Substitute Eq. () and Eq. () into Eq. (), to the leading order of (e^(-d/ξ)) (in the tight binding limit), the effective Hamiltonian reads

(OverTilde[H])=S H Sᵀ
=(Subscript[E, 0])(({{1, 2(e^(-d/ξ))}, {2(e^(-d/ξ)), 1}}))+𝒪( (e^(-2d/ξ))).

Block[{ε,S,H},
ε=1/2 ArcTanh[Exp[-x](1+x)];
S=Sqrt[Cosh[2ε]]{{Cosh[ε],-Sinh[ε]},{-Sinh[ε],Cosh[ε]}};
H={{1+2Exp[-2x],Exp[-x](3+x)},{Exp[-x](3+x),1+2Exp[-2x]}};
Simplify[S.H.Transpose[S]]]

Matrix multiplication.

or written in the tight binding basis as

The diagonal term (Ket[(OverTilde[1])])(Bra[(OverTilde[1])])+(Ket[(OverTilde[2])])(Bra[(OverTilde[2])]) describes that the electron trapped around an atom will experience the on-site energy (Subscript[E, 0]).

The off-diagonal term (Ket[(OverTilde[1])])(Bra[(OverTilde[2])])+(Ket[(OverTilde[2])])(Bra[(OverTilde[1])]) describes that the electron can tunnel from one atom to another, and the process is associated with an energy scale 2(Subscript[E, 0])(e^(-d/ξ)), which is also called the hopping energy.

Molecular orbitals (tight binding model)

In the tight biding basis, the Schrödinger equation reduces from a differential equation to a matrix equation

(Subscript[H, eff])(Ket[ψ])=E(Ket[ψ]),

or explicitly

(Subscript[E, 0])(({{1, 2(e^(-d/ξ))}, {2(e^(-d/ξ)), 1}}))(({{(Subscript[ψ, (OverTilde[1])])}, {(Subscript[ψ, (OverTilde[2])])}}))=E(({{(Subscript[ψ, (OverTilde[1])])}, {(Subscript[ψ, (OverTilde[2])])}})).

Solving the eigen problem, one finds two eigen states (Ket[(Subscript[ψ, ±])]) with the eigen energy (Subscript[E, ±])

Block[{orbit},
orbit[col_,dat_]:=Block[{L=2},
{{AbsoluteThickness[1],Line[{#,0}&/@{-L,L}]},
First@Plot[Total[#1 Exp[-2Abs[x-#2]]&@@@dat],{x,-L,L},PlotStyle->col,PlotRange->{-1,1},Filling->Axis],
Translate[{{White,EdgeForm[Black],Disk[{0,0},Offset[{4,4}]]},{Text["+"]}},{#2,0}]&@@@dat}];Plot[{-If[d>0,1+ProductLog[d Exp[-d]]/d,2]^2,-If[d>1,1+ProductLog[-d Exp[-d]]/d,Indeterminate]^2,-(1+2Exp[-d]),-(1-2Exp[-d])},{d,2,6},PlotStyle->{Blue,Red,Directive[Dashed,Blue],Directive[Dashed,Red]},Epilog->{Translate[Scale[orbit[Red,{{1/Sqrt[2],1},{-1/Sqrt[2],-1}}],{1/3,1/9},{0,0}],{4,-0.8}],
Translate[Scale[orbit[Blue,{{1/Sqrt[2],1},{1/Sqrt[2],-1}}],{1/3,1/9},{0,0}],{4,-1.2}]},GridLines->{None,{-1}},FrameLabel->{"\!\(d/ξ\)","\!\(E/\[LeftBracketingBar]E\_0\[RightBracketingBar]\)"},PlotRange->{-1.3,-0.7},ImageSize->200]]

- Solid curve: exact energy. 
- Dashed curve: approximate energy to the 1st order in (e^(-d/ξ)) (given by the tight binding model).

: even parity, bonding orbital (Ket[(SubPlus[ψ])]),
: odd parity, anti-bonding orbital (Ket[(SubMinus[ψ])])

1D Electronic Systems

Electron on 1D Lattice

H=-((ℏ^2)/(2m))(∂_x^2)+V(x)

Periodic potential: V(x+a)=V(x), i.e. the potential is invariant under  lattice translation x→x+a.

Block[{ξ=1.,a=0.1},
Graphics[{
{AbsoluteThickness[1],Arrow[{#,0}&/@{-4.5,5}],Text["\!\(,x\)",{5,0},{1,1}],Arrow[{-4.5,#}&/@{-1,1}],Text["\!\(V(x)\)",{-4.5,1},{-1,0}]},Translate[{{Line[{{-0.5,0},{-a,0},{-a,-1},{a,-1},{a,0},{0.5,0}}]},Translate[{{White,EdgeForm[Black],Disk[{0,0},Offset[{4,4}]]},{Text["+"]}},{0,0.5}]},{#,0}]&/@Range[-4,4]
},PlotRange->{{-5,5},All},ImageSize->250]
]

(BraKet[x, i])=(Subscript[ψ, i])(x)∝(e^(-|x-(Subscript[x, i])|/ξ))+(orthogonality corrections...),

Block[{ξ=1.,basis},
basis[op_]:=Block[{L=4,ψ,x},ψ=MatrixPower[Table[Exp[-#](1+#)&[Abs[x1-x2]/ξ],{x1,-L,L},{x2,-L,L}],-1/2,KroneckerDelta/@Range[-L,L]].(Exp[-Abs[x-#]/ξ]&/@Range[-L,L]);
First@Plot[ψ,{x,-L,L},Filling->Axis,FillingStyle->Directive[Black,Opacity[0.2op]],PlotStyle->Directive[Black,Opacity[op]],PlotRange->All]];
Graphics[{Translate[{basis[If[#==0,1,0.2]]},{#,0}]&/@Range[-4,4],
{Arrow[{#,0}&/@{-4.5,5}],Text["\!\(,x\)",{5,0},{1,-1}]},
Text["\!\(,i\)",{0,0},{0,1.5}]},PlotRange->{{-5,5},All},ImageSize->250]
]

The lattice translation correspond to (Ket[i])→(Ket[i+1]) for tight binding basis.

H=(∑_(i,j))(Ket[i])(Subscript[H, i,j])(Bra[j]),

The lattice translation symmetry implies

(Bra[i+1])H(Ket[j+1])=(Bra[i])H(Ket[j]),

i.e. (Subscript[H, i,j]) only depends on the index difference i-j.

Hermiticity of the Hamiltonian requires

(Subscript[H, i,j])=(Subscript[H, j,i]^( * )).

The Hamiltonian matrix must take the form of

(Subscript[H, i,j])=({{{, ({{-(Subscript[t, i-j]), i>j}, {(Subscript[ϵ, 0]), i==j}, {-(Subscript[t, j-i]^( * )), i<j}})}}),

or in the matrix form as

H=(({{⋱, ⋱, ⋱,  , ⋱}, {⋱, (Subscript[ϵ, 0]), -(Subscript[t, 1]^( * )), -(Subscript[t, 2]^( * )),  }, {⋱, -(Subscript[t, 1]), (Subscript[ϵ, 0]), -(Subscript[t, 1]^( * )), ⋱}, { , -(Subscript[t, 2]), -(Subscript[t, 1]), (Subscript[ϵ, 0]), ⋱}, {⋱,  , ⋱, ⋱, ⋱}})).

(Subscript[ϵ, 0]) - on-site potential energy.

(Subscript[t, n]) - nth neighbor hopping parameter (kinetic energy). (Subscript[t, n]) is expected to decay exponentially

(Subscript[t, n])∼(e^(-|(Subscript[x, n])-(Subscript[x, 0])|/ξ))=(e^(-n a /ξ)),

therefore (Subscript[t, n]) is often truncated to the first few leading ones (the few nearest neighbors hopping parameters).

Tight Binding Model and Energy Band

The quasi-momentum basis and the tight binding basis are related Fourier / inverse Fourier transforms

where N is the total number of sites in the lattice.

Assuming N is even, (∑_k) sums over

k=0,±((2π)/(N a)),…,±((2π m)/(N a)),…,±((2π (N/2-1))/(N a)),(π/a).

(∑_i) sums over i=1,2,...,N-1,N.

Substitute Eq. () into Eq. (), the Hamiltonian is diagonalized

Derive Eq. ().

(∑_i)(Ket[i+n])(Bra[i])=(∑_i)((N^(-1/2))(∑_k)(e^(-i k (Subscript[x, i+n])))(Ket[k]))((N^(-1/2))(∑_(k'))(e^(i k' (Subscript[x, i])))(Bra[k']))
=(1/N)(∑_(k,k'))(∑_i)(e^(-i (k (Subscript[x, i+n])-k'(Subscript[x, i]))))(Ket[k])(Bra[k'])
=(1/N)(∑_(k,k'))(∑_i)(e^(-i (k (Subscript[x, i+n])-k (Subscript[x, i])+k (Subscript[x, i])-k'(Subscript[x, i]))))(Ket[k])(Bra[k'])
=(1/N)(∑_(k,k'))(∑_i)(e^(-i k((Subscript[x, i+n])-(Subscript[x, i]))))(e^(-i (k-k')(Subscript[x, i])))(Ket[k])(Bra[k'])
=(∑_(k,k'))(e^(-i n k a))((1/N)(∑_i)(e^(-i (k-k')(Subscript[x, i]))))(Ket[k])(Bra[k'])
=(∑_(k,k'))(e^(-i n k a))(Subscript[δ, k,k'])(Ket[k])(Bra[k'])
=(∑_k)(e^(-i n k a))(Ket[k])(Bra[k]).

H=(Subscript[ϵ, 0])(∑_i)(Ket[i])(Bra[i])-(Subscript[t, 1])((∑_i)(Ket[i+1])(Bra[i])+(∑_i)(Ket[i-1])(Bra[i]))-(Subscript[t, 2])((∑_i)(Ket[i+2])(Bra[i])+(∑_i)(Ket[i-2])(Bra[i]))
=(∑_k)((Subscript[ϵ, 0])-(Subscript[t, 1])((e^(-i k a))+(e^(i k a)))-(Subscript[t, 2])((e^(-i 2k a))+(e^(i 2k a))))(Ket[k])(Bra[k])
=(∑_k)((Subscript[ϵ, 0])-2(Subscript[t, 1])cos(k a)-2(Subscript[t, 2])cos(2k a))(Ket[k])(Bra[k]).

Dispersion relation (Subscript[ϵ, k]) of the electron (also called band structure, band dispersion)

Manipulate[Plot[ϵ0-2t1 Cos[k]-2t2 Cos[2k],{k,-π,π},PlotStyle->Black,PlotRange->{-4,4},FrameLabel->{"\!\(k a\)","\!\(ϵ\_k\)"},FrameTicks->{{Automatic,Automatic},{Range[-1,1,1/2]π,Automatic}},ImageSize->180],{{ϵ0,0,"\!\(ϵ\_0\)"},-1,1,ImageSize->Small},{{t1,1,"\!\(t\_1\)"},0,1,ImageSize->Small},{{t2,0,"\!\(t\_2\)"},0,1,ImageSize->Small},Paneled->False,LabelStyle->"InlineFormula"]

The dispersion relation (Subscript[ϵ, k]) is periodic in the quasi-momentum, i.e.

(Subscript[ϵ, k+((2π)/a)])=(Subscript[ϵ, k]),

so the quasi-momentum k of electron can be taken to be in the first Brillouin zone k∈(-π/a,π/a].

Energy band: the allowed energy range of the electron. The electron can not take an energy outside the energy band (as there is no state outside the energy band for the electron to occupied).

Eigenstates (quasi-momentum basis states) in real space representation

(BraKet[x, k])=(N^(-1/2))(∑_i)(e^(i k (Subscript[x, i])))(BraKet[x, i])
=(N^(-1/2))(∑_i)(e^(i k a i))(Subscript[ψ, i])(x).

Tight binding basis (Subscript[ψ, i]) enveloped by a plane wave.

Manipulate[Plot[Block[{ξ=0.2},Re[Exp[-I t]Sum[Exp[I k i] Exp[-Abs[x-i]/ξ],{i,-8,8}]]],{x,-5,5},PlotRange->{-2,2},PlotStyle->Black,FrameLabel->{"\!\(,x\)","\!\((BraKet[x, k])\)"},ImageSize->180],{{k,0},-π,π,ImageSize->Small},{t,0,2π,ImageSize->Small},LabelStyle->"InlineFormula",Paneled->False]

Group velocity of the electron

Take the dispersion relation (Subscript[ϵ, k]) in Eq. (),

(Subscript[v, k])=((2a (Subscript[t, 1]))/ℏ)sin(k a)+((4a (Subscript[t, 2]))/ℏ)sin(2k a).

Manipulate[Plot[2 t1 Sin[k]+4 t2 Sin[2 k],{k,-π,π},PlotStyle->Black,PlotRange->{-4,4},FrameLabel->{"\!\(k a\)","\!\(ℏ v\_k/a\)"},FrameTicks->{{Automatic,Automatic},{Range[-1,1,1/2]π,Automatic}},ImageSize->180],{{t1,1,"\!\(t\_1\)"},0,1,ImageSize->Small},{{t2,0,"\!\(t\_2\)"},0,1,ImageSize->Small},Paneled->False,LabelStyle->"InlineFormula"]

The group velocity always vanishes ((Subscript[v, k])=0) at both the Brillouin zone center (k=0) and zone boundary (k=π/a), because they are inversion symmetric points (i.e. k≅-k) in the Brillouin zone. Since velocity is odd under inversion

(Subscript[v, -k])=-(Subscript[v, k]),

it must vanish at inversion symmetric momentum points.

The quasi-momentum basis states (Ket[k]) are standing wave states at k=0 and k=π/a, i.e. the electron wave does not propagate and the group velocity vanishes. The standing wave at k=π/a is caused by the strong back-scattering (Bragg scattering) of electron by the lattice.

Effective mass (Subscript[m, eff]) of the electron

Take the dispersion relation (Subscript[ϵ, k]) in Eq. (),

(Subscript[m, eff]^(-1))=((2(a^2)(Subscript[t, 1]))/(ℏ^2))cos(k a)+((8(a^2)(Subscript[t, 2]))/(ℏ^2))cos(2k a).

At band bottom (k=0), the effective mass is positive

(Subscript[m, eff])=((ℏ^2)/(2 (a^2)((Subscript[t, 1])+4(Subscript[t, 2])))).

At band top, the effective mass is negative. For example, for (Subscript[t, 2])<(Subscript[t, 1])/4, the band top is at k=π/a, where the effective mass is

(Subscript[m, eff])=((ℏ^2)/(2 (a^2)(-(Subscript[t, 1])+4(Subscript[t, 2])))).

When the effective mass is negative, the electron velocity decreases with increasing momentum. The negative effective mass is a peculiar feature of electron on the lattice.

Thermodynamic Limit and Density of State

(∑_k)=N a(∫_(-π)^π)((dk)/(2π))=V(∫_(-π)^π)((dk)/(2π)),

Extensive property that involves summation (∑_k) of states will be proportional to the volume V.

Intensive property (extensive property per volume) eliminates the volume dependence, and has well-defined thermodynamic limit.

(∑_k)=V∫dϵ (∑_((Subscript[ϵ, k])=ϵ))(1/(2π |(∂_k)(Subscript[ϵ, k])|))=V∫dϵ (∑_((Subscript[ϵ, k])=ϵ))(1/(2π ℏ |(Subscript[v, k])|)).

g(ϵ)=(∑_((Subscript[ϵ, k])=ϵ))(1/(2π ℏ |(Subscript[v, k])|)).

For (Subscript[ϵ, k])=-2t cos(k a), we have ℏ (Subscript[v, k])=2a t sin(k a),

g(ϵ)=(∑_(-2t cos(k a)=ϵ))(1/(4π a t |sin(k a)|))
=(1/(2π a t |sin(arccos(ϵ/(2t)))|))
=(1/(π a (Sqrt[(((2t))^2)-(ϵ^2)]))).

Plot[Re[1/(π Sqrt[4-ϵ^2])],{ϵ,-3,3},PlotStyle->Black,Filling->0,Exclusions->None,FrameLabel->{"\!\(ϵ\)","\!\(g(ϵ)\)"},ImageSize->180]

The divergences of density of state near band edges are examples of van Hove singularities.

The mid-band density of state saturates to g(0)=1/(2π a t).

Fermi Statistics and Band Filling

β=1/(Subscript[k, B])T, where T is the temperature.

μ is the chemical potential (Fermi energy).

Manipulate[Plot[1/(Exp[ϵ/T]+1),{ϵ,-5,5},PlotRange->{0,1},PlotRangePadding->Scaled[.05],PlotStyle->Black,FrameLabel->{"\!\(ϵ-μ\)","\!\(n\_F(ϵ)\)"},ImageSize->180],{{T,0.1},0.01,1,ImageSize->Small},Paneled->False,LabelStyle->"InlineFormula"]

The chemical potential is set by the total number (Subscript[N, e]) of electrons

The factor 2 in front of the summation takes into account the fact that each electron has two degenerated spin states.

Then the total internal energy E of the electronic system is

(Subscript[N, e])=2(∑_((Subscript[ϵ, k])<μ))1, E=2(∑_((Subscript[ϵ, k])<μ))(Subscript[ϵ, k]).

For (Subscript[ϵ, k])=-2t cos(k a), use the density of state g(ϵ) in Eq. (), the total number of electrons increases with μ from 0 to 2 as

(Subscript[N, e])=2V(∫_(-2t)^μ)dϵ g(ϵ)
=((2V)/(π a))(∫_(-2t)^μ)dϵ (1/(Sqrt[(((2t))^2)-(ϵ^2)]))
=((2N)/π)(arcsin((μ/(2 t)))+(π/2)).

Plot[2/π(ArcSin[μ/2]+π/2),{μ,-2,2},PlotStyle->Black,FrameLabel->{"\!\(μ/t\)","\!\(N\_e/N\)"},ImageSize->180]

and the total internal energy varies with μ as

E=2V(∫_(-2t)^μ)dϵ ϵ g(ϵ)
=((2V)/(π a))(∫_(-2t)^μ)dϵ (ϵ/(Sqrt[(((2t))^2)-(ϵ^2)]))
=-((2N)/π)(Sqrt[(((2t))^2)-(μ^2)]).

Plot[-(2/π)Sqrt[4-μ^2],{μ,-2,2},PlotStyle->Black,FrameLabel->{"\!\(μ/t\)","\!\(E/N\)"},ImageSize->180]

Filling fraction: the fraction of a energy band filled by electrons

ν=((Subscript[N, e])/(2N)),

where 2N is the number of states (including 2-fold spin degeneracy) in the energy band. N is the number of sites in the lattice.

When 0<ν<1, the many-body system of electrons forms a metal (in the non-interacting limit).

Row[Table[Graphics[{{AbsoluteThickness[1],Arrowheads[0.06],Arrow[{#,0}&/@{-π,π}],Arrow[{0,#}&/@{-2.5,2.5}],Text["\!\(,k\)",{π,0},{1,-1}],Text["\!\(ϵ\_k\)",{0,2.5},{1.5,1}]},
Line@Table[{k,-2Cos[k]},{k,-π,π,0.01π}],{AbsoluteThickness[5],Line@Table[{k,-2Cos[k]},{k,-π/2+δk,π/2+δk,0.01π}]},Text["\!\(ν=1/2\), "<>If[δk==0,"\!\(⟨j⟩=0\)","\!\(⟨j⟩≠0\)"],{0,2.5},{0,-1}],
Text[If[δk==0,"ground state","current-carrying"],{0,-2.5},{0,1}]},ImageSize->120],{δk,{0,π/16}}]]

In the many-body ground state: electrons occupy all states from the band bottom up to the designated filling fraction => total energy E is minimized.

A current-carrying excited state can be created by an imbalanced occupation of electrons, such that there are more electrons on the (Subscript[v, k])>0 side than those on the (Subscript[v, k])<0 side. => ⟨j⟩=e(∑_k)⟨(Subscript[v, k])⟩≠0. => Metal can conduct electricity.

However, if the band if empty (ν=0) or fully filled (ν=1), the current-carrying state can not be constructed. => The system can not conduct electricity, hence an insulator (also called band insulator).

Row[Table[Graphics[{{AbsoluteThickness[1],Arrowheads[0.06],Arrow[{#,0}&/@{-π,π}],Arrow[{0,#}&/@{-2.5,2.5}],Text["\!\(,k\)",{π,0},{1,-1}],Text["\!\(ϵ\_k\)",{0,2.5},{1.5,1}]},
Line@Table[{k,-2Cos[k]},{k,-π,π,0.01π}],If[kF!=0,{AbsoluteThickness[5],Line@Table[{k,-2Cos[k]},{k,-kF,kF,0.01π}]}],Text[If[kF==0,"\!\(ν=0\)","\!\(ν=1\)"]<>", \!\(⟨j⟩=0\)",{0,2.5},{0,-1}],
Text[If[kF==0,"empty band","fully-filled band"],{0,-2.5},{0,1}]},ImageSize->120],{kF,{0,π}}]]

Metallic bond: N atoms can gain energy jointly by allowing their electrons to delocalize and hybridize across the whole system. This energy gain binds all atoms together forming a metallic state.

Block[{β=10.,ϵ0=0.,t1=1.,t2=0.,ϵ,NE},
ϵ[k_]:=ϵ0-2t1 Cos[k]-2t2 Cos[2k];
NE[μ_]:=NIntegrate[{1,2ϵ[k]}/(Exp[β(ϵ[k]-μ)/2]+1),{k,-π,π}]/(2π);
ListLinePlot[Table[NE[μ],{μ,-4,4,0.2}],PlotStyle->Black,FrameLabel->{"\!\(ν=N\_e/2N\)","\!\(E/N\)"},ImageSize->180]]

Given (Subscript[ϵ, k]), one can calculate E v.s. (Subscript[N, e]) by tuning μ in Eq. ().

The total energy E is typically minimized at a some ν∈(0,1) => electrons forming a metallic state.

Su-Schrieffer-Heeger Model

Block[{h=0.3,d=0.08},
Graphics[{{LightGray,AbsoluteThickness[25],Line[{{1,-h},{2,h}}]},Table[If[EvenQ[i],#,Table[Translate[#,d s {2h(-1)^i,1}],{s,{-1,1}}]]&@Line[{{i,h(-1)^i},{i+1,h(-1)^(i+1)}}],{i,11}],Table[Translate[{EdgeForm[Black],Lighter[If[EvenQ[i],Blue,Red],0.6],Disk[{0,0},0.18]},{i,h(-1)^i}],{i,12}],
{Red,Text["A",{1,-h},{0,1.5}]},{Blue,Text["B",{2,h},{0,-1}]},Text["\!\(,v\)",{3.5,0},{-1,0.8}],Text["\!\(,w\)",{4.5,0},{-1,-0.8}]},ImageSize->250]]

Each unit cell (gray shaded) now contains two atoms (labeled by A and B). Each lattice site is labeled jointly by the unit cell index i and the sublattice label A/B, such that a set of tight binding basis can be defined

…,(Ket[i,A]), (Ket[i,B]),(Ket[i+1,A]), (Ket[i+1,B])…

The hopping parameters are different on the double bond and the single bond (double bond is shorter => larger hopping).

v: intra-unit-cell hopping,

w: inter-unit-cell hopping.

Column[Function[{F,gl},Block[{h=0.3,d=0.08},
Graphics[{{LightGray,AbsoluteThickness[25],Line[{{3,-h},{4,h}}]},Table[If[F[i],#,Table[Translate[#,d s {2h(-1)^i,1}],{s,{-1,1}}]]&@Line[{{i,h(-1)^i},{i+1,h(-1)^(i+1)}}],{i,11}],Table[Translate[{EdgeForm[Black],Lighter[If[EvenQ[i],Blue,Red],0.6],Disk[{0,0},0.18]},{i,h(-1)^i}],{i,12}],
Text["\!\(,v\)",{3.5,0},{-1,0.8}],Text["\!\(,w\)",{4.5,0},{-1,-0.8}],Text["\!\((v"<>gl<>"w)\)",{14,0}]},ImageSize->280]]]@@@{{EvenQ,">"},{OddQ,"<"}}]

(Ket[i,A])=(N^(-1/2))(∑_k)(e^(-i k (Subscript[x, i])))(Ket[k,A]),
(Ket[i,B])=(N^(-1/2))(∑_k)(e^(-i k (Subscript[x, i])))(Ket[k,B]),

Derive Eq. ().

H=-v(∑_i)((Ket[i,B])(Bra[i,A])+h.c.)-w(∑_i)((Ket[i+1,A])(Bra[i,B])+h.c.)
=-(v/N)(∑_i)(∑_(k,k'))((e^(-i k (Subscript[x, i])))(Ket[k,B])(Bra[k',A])(e^(i k' (Subscript[x, i])))+h.c.)-(w/N)(∑_i)(∑_(k,k'))((e^(-i k (Subscript[x, i+1])))(Ket[k,A])(Bra[k',B])(e^(i k' (Subscript[x, i])))+h.c.)
=-(v/N)(∑_(k,k'))(∑_i)((e^(-i (k-k')i ))(Ket[k,B])(Bra[k',A])+h.c.)-(w/N)(∑_(k,k'))(∑_i)((e^(-i (k-k')i))(e^(-i k))(Ket[k,A])(Bra[k',B])+h.c.)
=-v(∑_(k,k'))((Subscript[δ, k,k'])(Ket[k,B])(Bra[k',A])+h.c.)-w(∑_(k,k'))((Subscript[δ, k,k'])(e^(-i k))(Ket[k,A])(Bra[k',B])+h.c.)
=-v(∑_k)((Ket[k,B])(Bra[k,A])+h.c.)-w(∑_k)((e^(-i k))(Ket[k,A])(Bra[k,B])+h.c.)
=-v(∑_k)((Ket[k,B])(Bra[k,A])+h.c.)-w(∑_k)((e^(i k))(Ket[k,B])(Bra[k,A])+h.c.)
=(∑_k)((-v-w (e^(i k)))(Ket[k,B])(Bra[k,A])+h.c.).

Represent in the {(Ket[k,A]),(Ket[k,B])} basis, H takes the matrix form

H≏(Underscript[⊕, k])(Subscript[H, k]),(Subscript[H, k])≏(({{0, (Subscript[f, k]^( * ))}, {(Subscript[f, k]), 0}})),

No mixing between different momentum sectors => the Hamiltonian can be diagonalized in each momentum sector independently.

Goal: find eigen vectors

(Ket[k,m])=(Subscript[u, k,m,A])(Ket[k,A])+(Subscript[u, k,m,B])(Ket[k,B])
≏(({{(Subscript[u, k,m,A])}, {(Subscript[u, k,m,B])}})).

such that

(Subscript[H, k])(Ket[k,m])=(Subscript[ϵ, k,m])(Ket[k,m]),

or more explicitly

(({{0, (Subscript[f, k]^( * ))}, {(Subscript[f, k]), 0}}))(({{(Subscript[u, k,m,A])}, {(Subscript[u, k,m,B])}}))=(Subscript[ϵ, k,m])(({{(Subscript[u, k,m,A])}, {(Subscript[u, k,m,B])}})).

Manipulate[Plot[{-Sqrt[v^2+w^2+2v w Cos[k]],Sqrt[v^2+w^2+2v w Cos[k]]},{k,-π,π},PlotRange->{-3,3},PlotStyle->Black,PlotLabel->StringTemplate["\!\((v,w)=(`1`,`2`)\)"][v,w],FrameLabel->{"\!\( k\)","\!\(ϵ\_\(k,±\)\)"},FrameTicks->{{Automatic,Automatic},{Range[-π,π,π/2],Automatic}},ImageSize->180],{{v,1.5},0,2,ImageSize->Small},{{w,0.5},0,2,ImageSize->Small},Paneled->False,LabelStyle->"InlineFormula"]

There are two bands (labeled by ±), matching the number of sites in the unit cell.

Band gap: the range of energy between two bands. Electrons are forbidden to appear in the band gap (there is no state for the electron to occupy within the gap).

Band gap opens at the Brillouin zone boundary, with the gap size

2Δ=2|v-w|.

In the low-temperature (T→0) limit, electrons will fully occupy the lower band, forming a band insulator.

Row[Function[{v,w},Plot[{-Sqrt[v^2+w^2+2v w Cos[k]],Sqrt[v^2+w^2+2v w Cos[k]]},{k,-π,π},PlotRange->{-3,3},PlotStyle->{Directive[AbsoluteThickness[5],Black],Black},FrameLabel->{"\!\( k\)","\!\(ϵ\_\(k,±\)\)"},FrameTicks->{{Automatic,Automatic},{Range[-π,π,π/2],Automatic}},PlotLabel->StringTemplate["\!\((v,w)=(`1`,`2`)\)"][v,w],AspectRatio->0.8,ImageSize->150]]@@@{{0.5,1.5},{1.5,0.5}}]

The band structure does not distinguish v>w from v<w, more information is contained in the eigenstates.

(Subscript[f, k])/|(Subscript[f, k])|∈U(1)≅(S^1) is a phase factor, that can wind around as k∈(S^1) cycles over the first Brillouin zone.

The map f:(S^1)→(S^1) is classified by the homotopy group (Subscript[π, 1])((S^1))≅ℤ. Distinct homotopy classes are indexed by different winding number

v>w (intra > inter): ω=0 (trivial insulator),

v<w (intra < inter): ω=1 (topological insulator).

Manipulate[Block[{f},
f[k_]:=-v -w Exp[I k];
Graphics[{{AbsoluteThickness[1],Arrow[{#,0}&/@{-2,2}],Arrow[{0,#}&/@{-2,2}],Text["Re\!\(f\_k\)",{2,0},{1,-1}],Text["Im\!\(f\_k\)",{0,2},{-1.4,1}]},
{Red,Arrowheads[0.06],Arrow@ReIm[f/@Range[-π,π,0.01π]]},{Point[{0,0}],Point[{-v,0}],{AbsoluteThickness[1],Dashed,Line[ReIm/@{-v,-v+w Exp[2π I/3]}]},Text["\!\(,v\)",ReIm[-v/2],{0,0.7}],
Text["\!\(,w\)",ReIm[-v+w Exp[2π I/3]/2],{-1,0}]},Text[StringTemplate["\!\(ω=`1`\)"][Chop[NIntegrate[I  w/(v Exp[-I k]+w),{k,-π,π},AccuracyGoal->5]/(2π I)]],{2,2},{1,1}]},
PlotRange->2,PlotLabel->Style[StringTemplate["\!\((v,w)=(`1`,`2`)\)"][v,w],Black,"InlineFormula"],ImageSize->160]],{{v,0.5},0,2,ImageSize->Small},{{w,0.5},0,2,ImageSize->Small},Paneled->False,LabelStyle->"InlineFormula"]

A topological phase transition happens at v=w, where the system transition between the trivial insulator and the topological insulator, via gap closing and reopening.

(Subscript[𝒜, k,±]) is the Berry connection of the vector bundle,

2π (Subscript[ω, ±]) is the Berry phase accumulated along the Brillouin zone cycle, which is also called the Zak phase for the 1D case.

Column[Function[{F,gl},Block[{h=0.3,d=0.08},
Graphics[{{LightGray,AbsoluteThickness[25],Line[{{3,-h},{4,h}}]},Table[If[F[i],#,Table[Translate[#,d s {2h(-1)^i,1}],{s,{-1,1}}]]&@Line[{{i,h(-1)^i},{i+1,h(-1)^(i+1)}}],{i,11}],Table[Translate[{EdgeForm[Black],Lighter[If[EvenQ[i],Blue,Red],0.6],Disk[{0,0},0.18]},{i,h(-1)^i}],{i,12}],
Text["\!\(,v\)",{3.5,0},{-1,0.8}],Text["\!\(,w\)",{4.5,0},{-1,-0.8}],Text["\!\((v"<>gl<>"w)\)",{14,0}]},ImageSize->280]]]@@@{{EvenQ,">"},{OddQ,"<"}}]

H≏(({{0, v,  ,  ,  }, {v, 0, w,  ,  }, { , w, ⋱, ⋱,  }, { ,  , ⋱, ⋱, v}, { ,  ,  , v, 0}})).

Energy spectrum:

Manipulate[Block[{Ν=32,H},
H=2.Symmetrize[SparseArray[{Band[{1,2},{Ν-1,Ν}]->{v,w}},{Ν,Ν}]];
ListPlot[Sort@Eigenvalues@H,PlotStyle->Black,PlotLabel->StringTemplate["\!\((v,w)=(`1`,`2`)\)"][v,w],PlotRange->{{0,Ν+1},{-3,3}},FrameLabel->{"mode index \!\(,n\)","\!\(ϵ\_n\)"},ImageSize->180]],
{{v,1},0,2,ImageSize->Small},{{w,2},0,2,ImageSize->Small},Paneled->False,LabelStyle->"InlineFormula"]

When v<w (topological), two zero-energy modes are left in the energy gap.

Edge states (zero-energy model): assuming the states take the form of

(Ket[L/R])=(∑_i)((Subscript[u, i,A]^(L/R))(Ket[i,A])+(Subscript[u, i,B]^(L/R))(Ket[i,B])).

The wave amplitude (Subscript[u, i,A/B]^±) ( on A sublattice,  on B sublattice) are plotted for both zero modes:

Block[{Ν=32,v=1.,w=2.,H,S,us},
H=2Symmetrize[SparseArray[{Band[{1,2},{Ν-1,Ν}]->{v,w}},{Ν,Ν}]];
S=SparseArray[{Band[{1,1},{Ν,Ν}]->{1,-1}},{Ν,Ν}];
us=Cases[Thread@Eigensystem@(H+0.01S),{e_,u_}/;Abs[e]<0.1:>{u,u.S.u}];
Column[Append[Graphics[{If[EvenQ[#1],Blue,Red],Point[{#1,#2}],Opacity[0.5],Line[{{#1,#2},{#1,0}}]}&@@@Thread@{Range[Ν],#1},PlotRange->{{0,Ν+1},{-1,1}},Frame->True,FrameLabel->{None,"\!\(u\_\(i,A/B\)\%"<>If[#2>0,"L","R"]<>"\)"},AspectRatio->0.4,ImageSize->220]&@@@us,Style["          site index","Text"]],Alignment->Center]]

The wave amplitude decays exponentially from the edge with an localization length ξ, e.g.

|(Subscript[u, i,A]^±)|∼(e^(-(Subscript[x, i])/ξ)).

Chiral Symmetry

Sublattice projection operators

(Subscript[P, A])=(∑_i)(Ket[i,A])(Bra[i,A]),
(Subscript[P, B])=(∑_i)(Ket[i,B])(Bra[i,B]).

For any state (Ket[ψ]) of an electron, 
(Bra[ψ])(Subscript[P, A])(Ket[ψ])= probability to find the electron on sublattice A, 
(Bra[ψ])(Subscript[P, B])(Ket[ψ])= probability to find the electron on sublattice B.

Chiral operator

or, in the momentum space,

S=(∑_k)((Ket[k,A])(Bra[k,A])-(Ket[k,B])(Bra[k,B])),

Represent in the {(Ket[k,A]),(Ket[k,B])} basis, S takes the matrix form

S≏(Underscript[⊕, k])(Subscript[S, k]), (Subscript[S, k])≏(({{1, 0}, {0, -1}}))=(σ^3).

Under chiral transformation

({{{, ({{S(Ket[i,A])=(Ket[i,A]),}, {S(Ket[i,B])=-(Ket[i,B]).}})}}) ({{{, ({{(Subscript[S, k])(Ket[k,A])=(Ket[k,A]),}, {(Subscript[S, k])(Ket[k,B])=-(Ket[k,B]).}})}})

B sublattice state gets a minus sign.

The chiral transformation forms a (Subscript[ℤ, 2]) group, as (S^2)=𝟙 (where 𝟙 stands for the identity operator).

For any state (Ket[ψ]) of an electron, (Bra[ψ])S(Ket[ψ])= sublattice polarization.
Such that
(Bra[ψ])S(Ket[ψ])>0: more likely to find the electron on sublattice A, 
(Bra[ψ])S(Ket[ψ])<0: more likely to find the electron on sublattice B.

Chiral symmetry: a unitary symmetry that anti-commutes with the Hamiltonian.

or, in the momentum space,

∀k:(Subscript[S, k]) (Subscript[H, k])(Subscript[S, k])=-(Subscript[H, k]).

Check Eq. () and Eq. ().

S H S=-v(∑_i)(S(Ket[i,B])(Bra[i,A])S+h.c.)-w(∑_i)(S(Ket[i+1,A])(Bra[i,B])S+h.c.)
=-v(∑_i)(-(Ket[i,B])(Bra[i,A])+h.c.)-w(∑_i)(-(Ket[i+1,A])(Bra[i,B])+h.c.)
=-H.

(Subscript[S, k])(Subscript[H, k])(Subscript[S, k])=(Subscript[f, k])(Subscript[S, k])(Ket[k,B])(Bra[k,A])(Subscript[S, k])+(Subscript[f, k]^( * ))(Subscript[S, k])(Ket[k,A])(Bra[k,B])(Subscript[S, k])
=-(Subscript[f, k])(Ket[k,B])(Bra[k,A])-(Subscript[f, k]^( * ))(Ket[k,A])(Bra[k,B])
=-(Subscript[H, k]).

Using the property (S^2)=𝟙, Eq. () is equivalent to S H=-H S.

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

H S(Ket[(Subscript[ψ, n])])=-S H(Ket[(Subscript[ψ, n])])=-S ((Subscript[E, n])(Ket[(Subscript[ψ, n])]))=-(Subscript[E, n])S(Ket[(Subscript[ψ, n])]),

If (Subscript[E, n])≠0, (Ket[(Subscript[ψ, n])]) and S(Ket[(Subscript[ψ, n])]) are two distinct eigenstate of H with opposite eigen energies (Subscript[E, n]) and -(Subscript[E, n]) => finite energy eigenstates always comes in ±(Subscript[E, n]) pairs.

For example, in the momentum space, (Subscript[H, k])(Ket[k,±])=(Subscript[ϵ, k,±])(Ket[k,±]), the eigen states indeed comes in pairs, related by the chiral transformation,

(Subscript[S, k])(Ket[k,±])=(Ket[k,∓]),

and the corresponding eigen energies are opposite to each other (Subscript[ϵ, k,±])=±|(Subscript[f, k])|.

However, even when the lattice translation symmetry is broken (e.g. on an open chain or in the presence of disorder) that quasi-momentum is not well-defined, the pairing of positive and negative energy states still holds through out the spectrum.

If (Subscript[E, n])=0, (Ket[(Subscript[ψ, n])]) and S(Ket[(Subscript[ψ, n])]) are both zero-energy eigenstates of H, but they could be the same state. So there is no requirement for the zero energy states to come in pairs.

(Bra[(Subscript[ψ, n])])H S(Ket[(Subscript[ψ, n])])=-(Bra[(Subscript[ψ, n])])S H(Ket[(Subscript[ψ, n])])
=>(Subscript[E, n])(Bra[(Subscript[ψ, n])])S(Ket[(Subscript[ψ, n])])=-(Subscript[E, n])(Bra[(Subscript[ψ, n])])S (Ket[(Subscript[ψ, n])])
=>2(Subscript[E, n])(Bra[(Subscript[ψ, n])])S(Ket[(Subscript[ψ, n])])=0.

If (Subscript[E, n])≠0, then (Bra[(Subscript[ψ, n])])S(Ket[(Subscript[ψ, n])])=0, i.e. sublattice polarization should vanish for any finite energy eigenstate (the electron must appear on sublattice A and B with equal probability).

This statement is manifest for momentum space eigenstates

(Ket[k,±])=(1/(Sqrt[2]))((Ket[k,A])±((Subscript[f, k])/(|(Subscript[f, k])|))(Ket[k,B])).

If (Subscript[E, n])=0, there is no restriction on (Bra[(Subscript[ψ, n])])S(Ket[(Subscript[ψ, n])]). In fact, one can choose (Ket[(Subscript[ψ, n])]) to be the common eigenstate of S and H (as they commute in the zero-energy sector). Given that S only has two eigenvalues ±1, the zero-energy eigenstates can all be classified by

(Bra[(Subscript[ψ, n])])S(Ket[(Subscript[ψ, n])])=±1,

i.e. the electron in zero-energy state is either fully polarized to sublattice A or to sublattice B.

The topological edge states are such examples.

Topological Phases of Band Insulators

(No symmetry)* Two band insulator Hamiltonians H and H' are in the same (invertible) topological phases, if they can be deformed into each other without closing the band gap.

* strictly speaking, band insulators already assumes the U(1) symmetry that corresponds to the electron number conservation. So “no symmetry” here means no additional symmetry apart form U(1).

(With symmetry) Two band insulator Hamiltonians H and H' are in the same symmetry protected topological (SPT) phases, if they can be deformed into each other without closing the band gap and without breaking the symmetry.

H≏(Underscript[⊕, k])(Subscript[H, k]),(Subscript[H, k])≏(Subscript[h, μ])(k)(σ^μ),

As a 2×2 Hermitian matrix, (Subscript[H, k]) can always be decomposed to real combinations of Pauli matrices (σ^μ)

(σ^0)=(({{1, 0}, {0, 1}})),(σ^1)=(({{0, 1}, {1, 0}})),(σ^2)=(({{0, -i}, {i, 0}})),(σ^3)=(({{1, 0}, {0, -1}})),

with the combination coefficients (Subscript[h, μ])(k)∈ℝ (μ=0,1,2,3).

For example, the SSH model:

(Subscript[h, 0])(k)=0,(Subscript[h, 1])(k)=Re (Subscript[f, k]), (Subscript[h, 2])(k)=Im (Subscript[f, k]), (Subscript[h, 3])(k)=0.

Deforming H is to deform the (Subscript[h, μ])(k) functions, subject to the periodic condition (Subscript[h, μ])(k+2π)=(Subscript[h, μ])(k) (i.e. k∈(S^1)).

(Subscript[ϵ, k,±])=(Subscript[h, 0])(k)±|h(k)|,

Gapping condition. To maintain the band gap is to keep the vector h(k) finite length everywhere in the Brillouin zone

∀k: \[LeftBracketingBar]h(k)\[RightBracketingBar]>0.

or h(k)∈(ℝ^3)\{0}.

Deformation without symmetry. All smooth deformations of the map h:(S^1)→(ℝ^3)\{0} are classified by the homotopy group

which is a trivial group, meaning that all maps h can be deformed into each other => Without additional symmetry (other than U(1)), all insulators in 1D are in the same phase (the trivial phase).

∀k:(Subscript[S, k])(Subscript[H, k])(Subscript[S, k])=-(Subscript[H, k]),

(Subscript[h, 0])(k)=(Subscript[h, 3])(k)=0.

Derive Eq. () from Eq. (), given Eq. ().

(Subscript[H, k])≏(Subscript[h, 0])(k)(σ^0)+(Subscript[h, 1])(k)(σ^1)+(Subscript[h, 2])(k)(σ^2)+(Subscript[h, 3])(k)(σ^3),

(Subscript[S, k])(Subscript[H, k])(Subscript[S, k])≏(σ^3)((Subscript[h, 0])(k)(σ^0)+(Subscript[h, 1])(k)(σ^1)+(Subscript[h, 2])(k)(σ^2)+(Subscript[h, 3])(k)(σ^3))(σ^3)
=(Subscript[h, 0])(k)(σ^3)(σ^0)(σ^3)+(Subscript[h, 1])(k)(σ^3)(σ^1)(σ^3)+(Subscript[h, 2])(k)(σ^3)(σ^2)(σ^3)+(Subscript[h, 3])(k)(σ^3)(σ^3)(σ^3)
=(Subscript[h, 0])(k)(σ^0)-(Subscript[h, 1])(k)(σ^1)-(Subscript[h, 2])(k)(σ^2)+(Subscript[h, 3])(k)(σ^3).

(Subscript[h, 0])(k)(σ^0)-(Subscript[h, 1])(k)(σ^1)-(Subscript[h, 2])(k)(σ^2)+(Subscript[h, 3])(k)(σ^3)=-(Subscript[h, 0])(k)(σ^0)-(Subscript[h, 1])(k)(σ^1)-(Subscript[h, 2])(k)(σ^2)-(Subscript[h, 3])(k)(σ^3),

Deformation under the chiral symmetry. All smooth deformation of the map h:(S^1)->(ℝ^2)\{0} are classified by the homotopy group

which is the integer addition group. Each class is labeled by an integer topological index ω, which corresponds to the winding number of h around the origin in the (Subscript[h, 1])-(Subscript[h, 2]) plane.

Protected by the chiral symmetry, it is not possible to smoothly deform 1D insulators of different topological indices => topological insulators (ω≠0) and the trivial insulator (ω=0) must be separated by gap-closing topological phase transitions.

However, if the chiral symmetry is broken, the symmetry protection is lifted, and it becomes possible to smoothly deforming among all 1D insulators.

H=-v(∑_i)((Ket[i,B])(Bra[i,A])+h.c.)-w(∑_i)((Ket[i+1,A])(Bra[i,B])+h.c.)-u(∑_i)((Ket[i,A])(Bra[i,A])-(Ket[i,B])(Bra[i,B])).

H≏(Underscript[⊕, k])(Subscript[H, k]),(Subscript[H, k])≏(({{u, (Subscript[f, k]^( * ))}, {(Subscript[f, k]), -u}}))=Re (Subscript[f, k])(σ^1)+Im (Subscript[f, k])(σ^2)+u (σ^3),

The energy dispersion becomes

(Subscript[ϵ, k,±])=±(Sqrt[((|(Subscript[f, k])|)^2)+(u^2)])=(Sqrt[(v^2)+(w^2)+2v w cos k+(u^2)]).

Manipulate[Plot[{-Sqrt[v^2+w^2+2v w Cos[k]+u^2],Sqrt[v^2+w^2+2v w Cos[k]+u^2]},{k,-π,π},PlotRange->{-3,3},PlotStyle->Black,PlotLabel->StringTemplate["\!\((v,w,u)=(`1`,`2`,`3`)\)"][v,w,u],FrameLabel->{"\!\( k\)","\!\(ϵ\_\(k,±\)\)"},FrameTicks->{{Automatic,Automatic},{Range[-π,π,π/2],Automatic}},ImageSize->180],{{v,1.5},0,2,ImageSize->Small},{{w,0.5},0,2,ImageSize->Small},{{u,0},0,2,ImageSize->Small},Paneled->False,LabelStyle->"InlineFormula"]

Phase diagram. With the chiral symmetry, u=0, the trivial and topological phases must be separated by phase transitions. However, when the chiral symmetry is broken by u≠0, the topological transitions can be avoided / circumvented.

Block[{shade},
shade[txt_]:={{White,Translate[txt,0.02#]&/@CirclePoints[16]},txt};Graphics[{{GrayLevel[0.95],Rectangle[{-2.5,-1},{2.5,1}]},{Arrowheads[0.05],AbsoluteThickness[1],Arrow[{#,0}&/@{-2.5,2.5}],Arrow[{0,#}&/@{-1,1}],Text["\!\(w/v\)",{2.5,0},{1,-1}],Text["\!\(u/v\)",{0,1},{-1.5,1}]},
Translate[{Dashed,Arrowheads[{-0.04,0.04}],Arrow[ReIm[0.3Exp[I #]]&/@Range[0,π,0.1π]]},{#,0}]&/@{-1,1},
{AbsoluteThickness[2],Green,Line[{#,0}&/@{-1,1}],Orange,Line[{#,0}&/@{-2.5,-1}],Line[{#,0}&/@{1,2.2}]},
{Red,AbsolutePointSize[6],Point[{#,0}&/@{-1,1}]},shade/@{Text["trivial",{0,0}],Text["topo.",{1.7,0}],Text["topo.",{-1.7,0}]}},ImageSize->220]]

Protected by the chiral symmetry, topological edge mode must exist on the boundary (the interface between trivial and topological states), otherwise one can proliferate the featureless boundary to avoid bulk transition (which contradict with the definition of SPT phases).

Block[{h=0.3,d=0.08},
Graphics[{{LightGray,Disk[{9,-h},0.5]},Table[If[If[i<9,EvenQ,OddQ][i],#,Table[Translate[#,d s {2h(-1)^i,1}],{s,{-1,1}}]]&@Line[{{i,h(-1)^i},{i+1,h(-1)^(i+1)}}],{i,15}],Table[Translate[{EdgeForm[Black],Lighter[If[EvenQ[i],Blue,Red],0.6],Disk[{0,0},0.18]},{i,h(-1)^i}],{i,16}],
Text["trivial",{4,-h},{0,1.5}],Text["topological",{14,-h},{0,1.5}],Text["boundary",{9,-h},{0,1.5}]},ImageSize->280]]

Interaction reduced classification. The band homotopy classification indicates that 1D band insulators with charge conservation and chiral symmetries are ℤ-classified. However, this is not the full story when interaction is taken into account. The stable classification is (Subscript[ℤ, 4]) for 1D interacting fermionic SPT states in symmetry class AIII, meaning that the ω=4 insulator can be trivialized by interaction. []

Lukasz Fidkowski, Alexei Kitaev. The effects of interactions on the topological classification of free fermion systems. arXiv:0904.2197.

2D Electronic Systems

Electron on Square Lattice

(Subscript[x, i])=(Subscript[n, i,1])(Subscript[a, 1])+(Subscript[n, i,2])(Subscript[a, 2]), ({{{, ({{(Subscript[a, 1])=(1,0)}, {(Subscript[a, 2])=(0,1)}})}}),

Block[{L=2},
Graphics[{{Arrowheads[0.08],Arrow[{{0,0},#}]&/@{{1,0},{0,1}},Text["\!\(a\_1\)",{1,0},{1,1}],Text["\!\(a\_2\)",{0,1},{1,1}]},{Rotate[Translate[Line[{0,#}&/@{-L-0.5,L+0.5}],{#,0}]&/@Range[-L,L],# π/2,{0,0}]&/@{0,1}},Translate[{White,EdgeForm[Black],Disk[{0,0},0.12]},#]&/@Tuples[Range[-L,L],2],Text["\!\(,t\)",{0.5,-1},{0,0.8}],Text["\!\(,i\)",{0,0},{2,1}],Text["\!\(,j\)",{1,0},{-2,1}]},ImageSize->140]]

Define the nearest neighboring bond vectors (the vector that point from one site to its nearest neighboring site)

({{{, ({{(Subscript[δ, 1])=(Subscript[a, 1])=(1,0)}, {(Subscript[δ, 2])=(Subscript[a, 2])=(0,1)}, {(Subscript[δ, 3])=-(Subscript[a, 1])=(-1,0)}, {(Subscript[δ, 4])=-(Subscript[a, 2])=(0,-1)}})}}).

i,j are site indices, which also label the tight binding basis (Ket[i]).

⟨i,j⟩ denotes i,j are adjacent on the lattice. (∑_(⟨i,j⟩)) sums over all such nearest neighbor links.

N=(Subscript[L, x])×(Subscript[L, y]) is the total number of unit cells (sites). (Subscript[L, x]),(Subscript[L, y]) are the (linear) size of the lattice along x and y directions.

(∑_k) sums over momentum k in the first Brillouin zone (((-π,π])^(×2)) on a grid with spacings 2π/(Subscript[L, x]) and 2π/(Subscript[L, y]) along x and y directions.

Derive Eq. ().

H=-t(∑_(⟨i j⟩))((Ket[i])(Bra[j])+h.c.)
=-(t/N)(∑_(⟨i j⟩))(∑_(k,k'))((e^(-i k · (Subscript[x, i])))(Ket[k])(Bra[k'])(e^(i k' · (Subscript[x, j])))+h.c.)
=-(t/N)(∑_(k,k'))(∑_(⟨i j⟩))((e^(-i (k-k') · (Subscript[x, i])))(e^(i k' · ((Subscript[x, j])-(Subscript[x, i]))))(Ket[k])(Bra[k'])+h.c.).

(∑_(⟨i j⟩))=(1/2)(∑_i)(∑_(j∈nn(i))),

H=-t(∑_(k,k'))((1/N)(∑_i)(e^(-i (k-k') · (Subscript[x, i])))(1/2)(∑_(j∈nn(i)))(e^(i k' · ((Subscript[x, j])-(Subscript[x, i]))))(Ket[k])(Bra[k'])+h.c.)
=-t(∑_(k,k'))((Subscript[δ, k,k'])(1/2)(∑_(j∈nn(i)))(e^(i k' · ((Subscript[x, j])-(Subscript[x, i]))))(Ket[k])(Bra[k'])+h.c.)
=-t(∑_k)((1/2)(∑_(j∈nn(i)))(e^(i k · ((Subscript[x, j])-(Subscript[x, i]))))(Ket[k])(Bra[k])+h.c.)
=-t(∑_k)((1/2)(∑_a)(e^(i k·(Subscript[δ, a])))(Ket[k])(Bra[k])+h.c.)
=-t(∑_k)(∑_a)(e^(i k·(Subscript[δ, a])))(Ket[k])(Bra[k]).

Row[{ContourPlot[-2(Cos[kx]+Cos[ky]),{kx,-π,π},{ky,-π,π},PlotRangePadding->0.8,Epilog->{{#1,Point[#3],Text[Style[#2,12],#3,#4]}&@@@{{White,"Γ",{0,0},{1,1}},{Black,"X",{π,0},{-1.5,0}},{Black,"Y",{0,π},{0,-0.7}},{Black,"M",{π,π},{-1.2,-0.7}}}},GridLines->{{-π,π},{-π,π}},FrameLabel->{"\!\(k\_x\)","\!\(k\_y\)"},FrameTicks->{{Range[-π,π,π/2],None},{Range[-π,π,π/2],None}},BaselinePosition->Scaled[0.65],ImageSize->180],ContourPlot[y,{x,0,1},{y,-4,4},AspectRatio->8,PlotRangePadding->None,ImageMargins->2,FrameTicks->{{None,Range[-4,4,2]},{None,None}},PlotLabel->"\!\(ϵ\_k/t\)",ImageSize->{Automatic,120}]}]

High symmetry points of the Brillouin zone: Γ(0,0), X(π,0), Y(0,π), M(π,π).

The band structure is often represented by plotting (Subscript[ϵ, k]) along straight line segments connecting high symmetry points.

Block[{Es},
Es[{kx_,ky_}]:={-2(Cos[kx]+Cos[ky])};
BZPlot[Es,{"Y","Γ","X","M","Γ"},BZRules->{"Γ"->{0,0},"X"->{π,0},"Y"->{0,π},"M"->{π,π}},PlotStyle->Black,GridLinesStyle->Directive[Opacity[0.2],AbsoluteThickness[1],Black],AxesStyle->Directive[Opacity[0.2],AbsoluteThickness[1],Black],FrameLabel->{None,"\!\(ϵ\_k/t\)"},AspectRatio->0.5,ImageSize->200]]

ν=(((∑_(k∈BZ))Θ((Subscript[ϵ, F])-(Subscript[ϵ, k])))/((∑_(k∈BZ))1)),

Manipulate[Block[{β=10.,ϵ},
ϵ[{kx_,ky_}]:=-2(Cos[kx]+Cos[ky]);
ContourPlot[ϵ[{kx,ky}],{kx,-π,π},{ky,-π,π},Contours->{μ},ContourStyle->None,ContourShading->{Gray,White},FrameLabel->{"\!\(k\_x\)","\!\(k\_y\)"},FrameTicks->{{Range[-π,π,π],None},{Range[-π,π,π],None}},PlotLabel->StringTemplate["\!\(ν=`1`\)"][Round[NIntegrate[(1-Tanh[β (ϵ[{kx,ky}]-μ)/2])/2,{kx,-π,π},{ky,-π,π}]/(2π)^2,0.02]],ImageSize->160]],
{{μ,0},-4,4,ImageSize->Small},Paneled->False,LabelStyle->"InlineFormula"]

Symmetry and the general shape of Fermi surface. In the free space, electrons should have a circular Fermi surface, due to the SO(2) rotation symmetry. However, on a lattice, the full SO(2) rotation symmetry is broken to its discrete subgroup, such as (Subscript[ℤ, 4]) for the square lattice. So the Fermi surface only have four-fold symmetry on a square lattice.

Small Fermi surface remains circular. Because electrons are also waves, when the wave length is much larger than the lattice spacing (which is the case near the Γ point), the electron will not resolve the lattice structure (as the structure is below the diffraction limit of the electron wave), and will not experience the symmetry breaking effect caused by the lattice, hence the Fermi surface remains circular near the Γ point.

Expand the dispersion relation (Subscript[ϵ, k]) in Eq. () around the Γ point (k→0) (here the lattice constant a=1 has been restored)

(Subscript[ϵ, k])=-4t +t (a^2) (k^2)+𝒪((k^4))
=-4t+(((ℏ^2) (k^2))/(2(Subscript[m, eff])))+𝒪((k^4))

with the effective mass (Subscript[m, eff])=(ℏ^2)/(2t (a^2)) of the electron set by the hopping parameter t, similar to Eq. ().

Large Fermi surface can be distorted. As the Fermi surface approaches the Brillouin zone boundary, the lattice scattering effect gets more prominent, which distorts the Fermi surface to non-circular shape (such as the square shape at half filling).

“Broken” Fermi surface. Beyond half filling, the Fermi surface can even touch the Brillouin zone boundary.

Note that whenever a Fermi surface touches the zone boundary, it must do so perpendicularly, because the normal component of the group velocity must be zero at the zone boundary.

The Fermi surface seems to be broken into segments, however, it is actually still continuous since the Brillouin zone is periodic.

Electron v.s. hole Fermi surface. When the filling fraction approaches one, the Fermi surface looks circular again (if we look at the empty region). It can be viewed as filling full band with holes. Holes are anti-particles of electrons. High filling of electrons = low filling of holes => Fermi surface becomes circular as the hole filling is small.

Electron on Honeycomb Lattice

Block[{L=1,a1,a2,Rs},
a1={Sqrt[3],0};
a2={Sqrt[3],3}/2;
Rs=Tuples[Range[-L-2,L+2],2].{a1,a2};
Graphics[{{Translate[Line[{{0,0},#}]&/@CirclePoints[{1,π/2},3],#+{Sqrt[3],1}]&/@Rs},{Arrowheads[0.06],{Arrow[{{0,0},#1}],Text[StringTemplate["\!\(\*StyleBox[\"a\",Bold]\_`1`\)"][#2],#1,#3]}&@@@{{a1,"1",{0,0.8}},{a2,"2",{-1.,-0.2}}},
Translate[{Arrow[{{0,0},#1}],Text[StringTemplate["\!\(δ\_`1`\)"][#2],#1,{-0.1+1.4Sign[#1[[1]]],0.3-0.8Sign[#1[[2]]]}]}&@@@Thread[{CirclePoints[{1,-π/2},3],Range[3]}],{-Sqrt[3]/2,1/2}]},
{EdgeForm[Black],Lighter[Blue,0.6],Disk[#+{Sqrt[3]/2,1/2},0.15]&/@Rs,Lighter[Red,0.6],Disk[#+{Sqrt[3],1},0.15]&/@Rs},
{Blue,Text["\!\(,B\)",{Sqrt[3]/2,1/2},{0,-0.8}],Red,Text["\!\(,A\)",{Sqrt[3]/2,-1/2},{0,1.2}]}},PlotRange->2L,PlotRangeClipping->True,PlotRangePadding->1,ImagePadding->{{Automatic, 5}, {Automatic, 5}},ImageSize->140]]

(Subscript[x, i])=(Subscript[x, n,s])=(Subscript[n, 1])(Subscript[a, 1])+(Subscript[n, 2])(Subscript[a, 2])+(Subscript[r, s]),

({{{, ({{(Subscript[a, 1])=((Sqrt[3]),0)}, {(Subscript[a, 2])=(((Sqrt[3])/2),(3/2))}})}}),({{{, ({{(Subscript[r, A])=((Sqrt[3]),1)}, {(Subscript[r, B])=(((Sqrt[3])/2),(1/2))}})}}).

The site index i↔(n,s) contains the unit cell and sublattice labels jointly.

Introduce the nearest neighboring displacement vectors for later convenience

({{{, ({{(Subscript[δ, 1])=(0,-1)}, {(Subscript[δ, 2])=(((Sqrt[3])/2),(1/2))}, {(Subscript[δ, 3])=(-((Sqrt[3])/2),(1/2))}})}}).

The tight binding basis (Ket[i]) is defined on each site, which may also be written as (Ket[n,s]) by the unit cell and sublattice labels.

(∑_(⟨i j⟩)) sums over nearest neighboring links ⟨i,j⟩ on the honeycomb lattice.

N is the number of unit cells.

The quasi-momentum is summed over the first Brillouin zone, which is a hexagon of side length ((4π)/(3(Sqrt[3]))).

Graphics[{GraphicsComplex[Join[4π/(3Sqrt[3])CirclePoints[6],ReIm/@{0,(2π/3)Exp[I π/6],(4π/3)Exp[-I π/6],(4π/3)I}],{{LightGray,EdgeForm[Black],Polygon[Range[6]]},Point@{2,3,7,8,9,10},Text@@@{{"Γ",7,{1,1}},{"K",2,{-1.5,0.2}},{"K'",3,{-1.5,-0.2}},{"M",8,{-1.5,-0.4}},{"Γ",9,{-1.5,0.2}},{"\!\(b\_1\)",9,{1.5,0.5}},{"Γ",10,{1,-0.4}},{"\!\(b\_2\)",10,{-1.5,1.5}}},
{Arrowheads[0.08],Arrow[{7,9}],Arrow[{7,10}]}}],
Text[Style["\!\(\(4π\)\/\(3\@3\)\)","InlineFormula"],{0,-2π/3},{0,1.2}]},ImageSize->120]

Recall that the first Brillouin zone is constructed as the Wigner-Seitz cell of the reciprocal lattice, spanned by the reciprocal basis

({{{, ({{(Subscript[b, 1])=(((2π)/(Sqrt[3])),-((2π)/3))}, {(Subscript[b, 2])=(0,((4π)/3))}})}}).

High symmetry points

Γ:(0,0)+(Subscript[m, 1])(Subscript[b, 1])+(Subscript[m, 2])(Subscript[b, 2]), 
K:(((4π)/(3(Sqrt[3]))),0)+(Subscript[m, 1])(Subscript[b, 1])+(Subscript[m, 2])(Subscript[b, 2]),
K':(-((4π)/(3(Sqrt[3]))),0)+(Subscript[m, 1])(Subscript[b, 1])+(Subscript[m, 2])(Subscript[b, 2]),
M:(0,((2π)/3))+(Subscript[m, 1])(Subscript[b, 1])+(Subscript[m, 2])(Subscript[b, 2]).

Derive Eq. ().

H=(∑_(i,j))(Ket[i])(Subscript[H, i,j])(Bra[j]),

(Subscript[H, i,j])=H((Subscript[x, i]),(Subscript[x, j]))=H((Subscript[x, (Subscript[n, i])])+(Subscript[r, (Subscript[s, i])]),(Subscript[x, (Subscript[n, j])])+(Subscript[r, (Subscript[s, j])])).

(Subscript[H, i,j])=H((Subscript[x, (Subscript[n, j])])-(Subscript[x, (Subscript[n, i])]),(Subscript[r, (Subscript[s, i])]),(Subscript[r, (Subscript[s, j])])).

H=(∑_(i,j))(Ket[i])(Subscript[H, i,j])(Bra[j])
=(∑_(i,j))(∑_(k,k'))(e^(-i k · (Subscript[x, i])))(Ket[k,(Subscript[s, i])])(Subscript[H, i,j])(Bra[k',(Subscript[s, j])])(e^(i k' · (Subscript[x, j])))
=(∑_(k,k'))(∑_(i,j))(Subscript[H, i,j])(e^(-i (k-k') · (Subscript[x, i])))(e^(i k' · ((Subscript[x, j])-(Subscript[x, i]))))(Ket[k,(Subscript[s, i])])(Bra[k',(Subscript[s, j])])
=(∑_(k,k'))(∑_(Subscript[n, i]))(∑_(Subscript[n, j]))(∑_((Subscript[s, i]),(Subscript[s, j])))H((Subscript[x, (Subscript[n, j])])-(Subscript[x, (Subscript[n, i])]),(Subscript[r, (Subscript[s, i])]),(Subscript[r, (Subscript[s, j])]))(e^(-i (k-k') · (Subscript[x, (Subscript[n, i])])))(e^(-i (k-k') · (Subscript[r, (Subscript[s, i])])))(e^(i k' · ((Subscript[x, j])-(Subscript[x, i]))))(Ket[k,(Subscript[s, i])])(Bra[k',(Subscript[s, j])])
=(∑_(k,k'))(∑_(Subscript[n, j]))(∑_((Subscript[s, i]),(Subscript[s, j])))H((Subscript[x, (Subscript[n, j])])-(Subscript[x, (Subscript[n, i])]),(Subscript[r, (Subscript[s, i])]),(Subscript[r, (Subscript[s, j])]))(e^(-i (k-k') · (Subscript[r, (Subscript[s, i])])))(e^(i k' · ((Subscript[x, j])-(Subscript[x, i]))))(Ket[k,(Subscript[s, i])])(Bra[k',(Subscript[s, j])])((∑_(Subscript[n, i]))(e^(-i (k-k') · (Subscript[x, (Subscript[n, i])]))))
=(∑_(k,k'))(∑_(Subscript[n, j]))(∑_((Subscript[s, i]),(Subscript[s, j])))H((Subscript[x, (Subscript[n, j])])-(Subscript[x, (Subscript[n, i])]),(Subscript[r, (Subscript[s, i])]),(Subscript[r, (Subscript[s, j])]))(e^(-i (k-k') · (Subscript[r, (Subscript[s, i])])))(e^(i k' · ((Subscript[x, j])-(Subscript[x, i]))))(Ket[k,(Subscript[s, i])])(Bra[k',(Subscript[s, j])])(Subscript[δ, k,k'])
=(∑_k)(∑_(Subscript[n, j]))(∑_((Subscript[s, i]),(Subscript[s, j])))H((Subscript[x, (Subscript[n, j])])-(Subscript[x, (Subscript[n, i])]),(Subscript[r, (Subscript[s, i])]),(Subscript[r, (Subscript[s, j])]))(e^(i k · ((Subscript[x, j])-(Subscript[x, i]))))(Ket[k,(Subscript[s, i])])(Bra[k,(Subscript[s, j])])

H=(∑_k)(∑_(Subscript[n, j]))(∑_((Subscript[s, i]),(Subscript[s, j])))(Subscript[H, i,j])(e^(i k · ((Subscript[x, j])-(Subscript[x, i]))))(Ket[k,(Subscript[s, i])])(Bra[k,(Subscript[s, j])]).

Represent in the {(Ket[k,A]),(Ket[k,B])} basis, H takes the matrix form

H≏(Underscript[⊕, k])(Subscript[H, k]),(Subscript[H, k])≏(({{0, (Subscript[f, k]^( * ))}, {(Subscript[f, k]), 0}})),

Block[{δ,f},
δ=<|1->{0,-1},2->{Sqrt[3]/2,1/2},3->{-Sqrt[3]/2,1/2}|>;
f[k_]:=-(Exp[I k.δ[1]]+Exp[I k.δ[2]]+Exp[I k.δ[3]]);
Row[{Show[Plot3D[{-Abs[f[{kx,ky}]],Abs[f[{kx,ky}]]},{kx,-4π/(3Sqrt[3]),4π/(3Sqrt[3])},{ky,-2π/3,2π/3},ColorFunction->"BlueRedSplit",Axes->None,BoxRatios->{1,1,1.5},MeshFunctions->{#3&},RegionFunction->Function[{kx,ky,e},Abs[kx]/(4π/(3Sqrt[3]))+Abs[ky]/(4π/3)<=1],ViewPoint->{0.4, -2, 0.2},ImageSize->{Automatic,200}],
Graphics3D[GraphicsComplex[Append[4π/(3Sqrt[3])#,0]&/@CirclePoints[6],{Line[{1,2,3,4,5,6,1}]}]]],
Column@Table[Graphics[{First@ContourPlot[ s Abs[f[{kx,ky}]],{kx,-4π/(3Sqrt[3]),4π/(3Sqrt[3])},{ky,-2π/3,2π/3},BoundaryStyle->Black,
ColorFunction->Function[e,ColorData["BlueRedSplit"][Rescale[e,{-3,3}]]],ColorFunctionScaling->False,RegionFunction->Function[{kx,ky,e},Abs[kx]/(4π/(3Sqrt[3]))+Abs[ky]/(4π/3)<=1]],
{#1,Point[#3],Text[Style[#2,12],#3,#4]}&@@@{{White,"Γ",{0,0},{1,1}},{Black,"K",{4π/(3Sqrt[3]),0},{-1.5,0}},{Black,"K'",{2π/(3Sqrt[3]),2π/3},{-1.2,-0.7}},{Black,"M",{π/Sqrt[3],π/3},{-1.2,-0.7}}}},ImageSize->100],{s,{1,-1}}],
ContourPlot[y,{x,0,1},{y,-3,3},Contours->11,AspectRatio->8,PlotRangePadding->None,ImageMargins->2,FrameTicks->{{None,Range[-3,3,1]},{None,None}},PlotLabel->"\!\(ϵ\_k/t\)",ImageSize->{Automatic,120}]}]]

Block[{δ,f,Es},
δ=<|1->{0,-1},2->{Sqrt[3]/2,1/2},3->{-Sqrt[3]/2,1/2}|>;
f[k_]:=-(Exp[I k.δ[1]]+Exp[I k.δ[2]]+Exp[I k.δ[3]]);
Es[k_]:={1,-1}Abs[f[k]];
BZPlot[Es,{"Γ","K","M","Γ","K'","M"},BZRules->{"Γ"->{0,0},"K"->{4π/(3Sqrt[3]),0},"K'"->{2π/(3Sqrt[3]),2π/3},"M"->{π/Sqrt[3],π/3}},PlotStyle->Black,GridLinesStyle->Directive[Opacity[0.2],AbsoluteThickness[1],Black],AxesStyle->Directive[Opacity[0.2],AbsoluteThickness[1],Black],FrameLabel->{None,"\!\(ϵ\_k/t\)"},AspectRatio->0.45,ImageSize->240]]

Two bands touch at the Brillouin zone corners (K and K' points).

(Subscript[H, K+κ])≏(Subscript[v, F])(({{0, (Subscript[κ, x])-i (Subscript[κ, y])}, {(Subscript[κ, x])+i (Subscript[κ, y]), 0}}))+𝒪((κ^2)),
=(Subscript[v, F])((Subscript[κ, x])(σ^1)+(Subscript[κ, y])(σ^2))+𝒪((κ^2)),

Block[{δ,f},
δ=<|1->{0,-1},2->{Sqrt[3]/2,1/2},3->{-Sqrt[3]/2,1/2}|>;
f[k_]:=-t(Exp[I k.δ[1]]+Exp[I k.δ[2]]+Exp[I k.δ[3]]);
Simplify@Series[f[{4π/(3Sqrt[3]),0}+κ{x,y}],{κ,0,1}]]

Expansion.

(Subscript[H, K'+κ])≏(Subscript[v, F])(({{0, -(Subscript[κ, x])-i (Subscript[κ, y])}, {-(Subscript[κ, x])+i (Subscript[κ, y]), 0}}))+𝒪((κ^2)),
=(Subscript[v, F])(-(Subscript[κ, x])(σ^1)+(Subscript[κ, y])(σ^2))+𝒪((κ^2)),

Block[{δ,f},
δ=<|1->{0,-1},2->{Sqrt[3]/2,1/2},3->{-Sqrt[3]/2,1/2}|>;
f[k_]:=-t(Exp[I k.δ[1]]+Exp[I k.δ[2]]+Exp[I k.δ[3]]);
Simplify@Series[f[{-4π/(3Sqrt[3]),0}+κ{x,y}],{κ,0,1}]]

Expansion.

These effective Hamiltonians describe 2D Dirac fermions of opposite chirality.

When the chemical potential is set to the zero energy (at charge neutrality), the electronic system is in the Dirac semi-metal phase.

Dirac cones (Dirac fermion dispersion). Near K and K' points, the energy is linear in the momentum deviation κ.

mimicking a relativistic particle with Fermi velocity

(Subscript[v, F])=((3t)/2).

Plot3D[{3/2Sqrt[kx^2+ky^2],-3/2Sqrt[kx^2+ky^2]},{kx,-1,1},{ky,-1,1},ColorFunction->"BlueRedSplit",Axes->None,BoxRatios->{1,1,1.5},MeshFunctions->{#3&,Arg[#1+I#2]&},RegionFunction->Function[{kx,ky,e},kx^2+ky^2<=1],ViewPoint->{2, -4, 3},ImageSize->{Automatic,140}]

The points K and K' are also called Dirac points, where the Dirac cones locate and the band gap close linearly. There are always even number of Dirac points in the Brillouin zone, known as the Fermion doubling.

Haldane Model

Unlike real hopping, imaginary hopping is directional, as

(SuperDagger[(i (Ket[i])(Bra[j]))])=-i (Ket[j])(Bra[i]).

Rule: the electron hopping from j to i (i.e. the (Ket[i])(Bra[j]) term) has a hopping coefficient +i if j to i is along the bond direction.

The bond directions are assigned in the following pattern

Block[{L=1,a1,a2,Rs},
a1={Sqrt[3],0};
a2={Sqrt[3],3}/2;
Rs=Tuples[Range[-L-2,L+2],2].{a1,a2};
Graphics[{{Arrowheads[{{0.07,0.65}}],Translate[{If[Norm[#]==0,Gray,LightGray],Rotate[Arrow[{# Sqrt[3]/2,1/2}&/@{-1,1}],# 2π/6,{0,0}]&/@Range[6]},#]&/@Rs},
Translate[{Arrowheads[0.06],Arrow[{{0,0},#1}],Text[StringTemplate["\!\(δ\_`1`\)"][#2],#1,{-0.1+1.4Sign[#1[[1]]],0.3-0.8Sign[#1[[2]]]}]}&@@@Thread[{CirclePoints[{1,-π/2},3],Range[3]}],{-Sqrt[3],2}],{Translate[Line[{{0,0},#}]&/@CirclePoints[{1,π/2},3],#+{Sqrt[3],1}]&/@Rs},{EdgeForm[Black],Lighter[Blue,0.6],Disk[#+{Sqrt[3]/2,1/2},0.15]&/@Rs,Lighter[Red,0.6],Disk[#+{Sqrt[3],1},0.15]&/@Rs}},PlotRange->2L,PlotRangeClipping->True,PlotRangePadding->1,ImagePadding->{{Automatic, 5}, {Automatic, 5}},ImageSize->140]]

Represent in the {(Ket[k,A]),(Ket[k,B])} basis, H takes the matrix form

H≏(Underscript[⊕, k])(Subscript[H, k]),(Subscript[H, k])≏(({{(Subscript[g, k]), (Subscript[f, k]^( * ))}, {(Subscript[f, k]), -(Subscript[g, k])}})),

The band dispersions are modified to

Manipulate[Block[{δ,f,g,Es},
δ=<|1->{0,-1},2->{Sqrt[3]/2,1/2},3->{-Sqrt[3]/2,1/2}|>;
f[k_]:=-(Exp[I k.δ[1]]+Exp[I k.δ[2]]+Exp[I k.δ[3]]);
g[k_]:=-2λ(Sin[k.(δ[1]-δ[2])]+Sin[k.(δ[2]-δ[3])]+Sin[k.(δ[3]-δ[1])]);
Es[k_]:={1,-1}Sqrt[Abs[f[k]]^2+g[k]^2];
BZPlot[Es,{"Γ","K","M","Γ","K'","M"},BZRules->{"Γ"->{0,0},"K"->{4π/(3Sqrt[3]),0},"K'"->{2π/(3Sqrt[3]),2π/3},"M"->{π/Sqrt[3],π/3}},PlotStyle->Black,PlotLabel->StringTemplate["\!\(λ/t=`1`\)"][Round[λ,0.01]],GridLinesStyle->Directive[Opacity[0.2],AbsoluteThickness[1],Black],AxesStyle->Directive[Opacity[0.2],AbsoluteThickness[1],Black],FrameLabel->{None,"\!\(ϵ\_k/t\)"},AspectRatio->0.45,ImageSize->220]],{{λ,0.1},-1./(3Sqrt[3]),1./(3Sqrt[3]),ImageSize->Small},Paneled->False,LabelStyle->"InlineFormula"]

The Haldane hopping term λ opens a gap between the two bands, gapping out Dirac fermions at K and K'.

(Subscript[H, K+κ])≏(({{m (Subscript[v, F]^2), (Subscript[v, F])((Subscript[κ, x])-i (Subscript[κ, y]))}, {(Subscript[v, F])((Subscript[κ, x])+i (Subscript[κ, y])), -m (Subscript[v, F]^2)}}))+𝒪((κ^2)),
=(Subscript[v, F])((Subscript[κ, x])(σ^1)+(Subscript[κ, y])(σ^2))+m (Subscript[v, F]^2)(σ^3)+𝒪((κ^2)),

Block[{δ,f,g,Es},
δ=<|1->{0,-1},2->{Sqrt[3]/2,1/2},3->{-Sqrt[3]/2,1/2}|>;
f[k_]:=-t(Exp[I k.δ[1]]+Exp[I k.δ[2]]+Exp[I k.δ[3]]);
g[k_]:=-2λ(Sin[k.(δ[1]-δ[2])]+Sin[k.(δ[2]-δ[3])]+Sin[k.(δ[3]-δ[1])]);
Simplify@Series[Through[{f,g}[{4π/(3Sqrt[3]),0}+κ{x,y}]],{κ,0,1}]]

Expansion.

(Subscript[H, K'+κ])≏(({{-m (Subscript[v, F]^2), (Subscript[v, F])(-(Subscript[κ, x])-i (Subscript[κ, y]))}, {(Subscript[v, F])(-(Subscript[κ, x])+i (Subscript[κ, y])), m (Subscript[v, F]^2)}}))+𝒪((κ^2)),
=(Subscript[v, F])(-(Subscript[κ, x])(σ^1)+(Subscript[κ, y])(σ^2))-m (Subscript[v, F]^2)(σ^3)+𝒪((κ^2)),

Block[{δ,f,g,Es},
δ=<|1->{0,-1},2->{Sqrt[3]/2,1/2},3->{-Sqrt[3]/2,1/2}|>;
f[k_]:=-t(Exp[I k.δ[1]]+Exp[I k.δ[2]]+Exp[I k.δ[3]]);
g[k_]:=-2λ(Sin[k.(δ[1]-δ[2])]+Sin[k.(δ[2]-δ[3])]+Sin[k.(δ[3]-δ[1])]);
Simplify@Series[Through[{f,g}[{-4π/(3Sqrt[3]),0}+κ{x,y}]],{κ,0,1}]]

Expansion.

(Subscript[v, F])=((3t)/2), m=((3(Sqrt[3])λ)/(Subscript[v, F]^2)).

Plot3D[{3/2Sqrt[kx^2+ky^2+0.15],-3/2Sqrt[kx^2+ky^2+0.15]},{kx,-1,1},{ky,-1,1},ColorFunction->"BlueRedSplit",Axes->None,BoxRatios->{1,1,1.5},MeshFunctions->{#3&,Arg[#1+I#2]&},RegionFunction->Function[{kx,ky,e},kx^2+ky^2<=1],ViewPoint->{2, -4, 3},ImageSize->{Automatic,140}]

When δ k is small, the energy expands to

(Subscript[ϵ, K/K'+κ,±])=±(m (Subscript[v, F]^2)+((κ^2)/(2m))+…),

which justifies the parameter m as the mass of the Dirac fermion (also the effective mass of electron/hole near the band edge).

Even if starting from non-relativistic electrons in graphene, relativistic behaviors can emerge at low-energy. This provides opportunities to simulate certain aspects of fundamental particle physics in condensed matter systems.

Chern Insulator and Integer Quantum Hall Effect

(Subscript[H, k])(Ket[k,±])=(Subscript[ϵ, k,±])(Ket[k,±]).

Berry connection of the vector bundle (Ket[k,±]), like a gauge potential in the momentum space,

At each momentum point k, the Berry connection has two components

(Subscript[𝒜, k,±]^x)=-i(Bra[k,±])(∂_(Subscript[k, x]))(Ket[k,±]),
(Subscript[𝒜, k,±]^y)=-i(Bra[k,±])(∂_(Subscript[k, y]))(Ket[k,±]).

Berry curvature: curl of Berry connection (taking the z component), like a magnetic field in the momentum space,

or in term of components

(Subscript[ℱ, k,±])=(∂_(Subscript[k, x]))(Subscript[𝒜, k,±]^y)-(∂_(Subscript[k, y]))(Subscript[𝒜, k,±]^x).

For 2×2 Hamiltonian of the form (Subscript[H, k])≏h(k)·σ, the Berry curvature is given by

(Subscript[ℱ, k,±])=±(1/(2((|h(k)|)^3)))h(k)·((∂_(Subscript[k, x]))h(k)×(∂_(Subscript[k, y]))h(k)).

Prove Eq. ().

Tr ρ(∂_μ)ρ(∂_ν)ρ
=Tr(Ket[ψ])(Bra[ψ])((Ket[(∂_μ)ψ])(Bra[ψ])+(Ket[ψ])(Bra[(∂_μ)ψ]))((Ket[(∂_ν)ψ])(Bra[ψ])+(Ket[ψ])(Bra[(∂_ν)ψ]))
=(BraKet[ψ, (∂_μ)ψ])(BraKet[ψ, (∂_ν)ψ])(BraKet[ψ, ψ])+(BraKet[ψ, (∂_μ)ψ])(BraKet[ψ, ψ])(BraKet[(∂_ν)ψ, ψ])+(BraKet[ψ, ψ])(BraKet[(∂_μ)ψ, (∂_ν)ψ])(BraKet[ψ, ψ])+(BraKet[ψ, ψ])(BraKet[(∂_μ)ψ, ψ])(BraKet[(∂_ν)ψ, ψ])
=(BraKet[(∂_μ)ψ, (∂_ν)ψ])+(BraKet[ψ, (∂_μ)ψ])(BraKet[ψ, (∂_ν)ψ])+(BraKet[ψ, (∂_μ)ψ])(BraKet[(∂_ν)ψ, ψ])+(BraKet[(∂_μ)ψ, ψ])(BraKet[(∂_ν)ψ, ψ]).

Tr ρ((∂_μ)ρ(∂_ν)ρ-(∂_ν)ρ(∂_μ)ρ)
=(BraKet[(∂_μ)ψ, (∂_ν)ψ])-(BraKet[(∂_ν)ψ, (∂_μ)ψ]).

(Subscript[𝒜, μ])=-i (BraKet[ψ, (∂_μ)ψ]),

(Subscript[ℱ, μ,ν])=(∂_μ)(Subscript[𝒜, ν])-(∂_ν)(Subscript[𝒜, μ])
=-i ((∂_μ)(BraKet[ψ, (∂_ν)ψ])-(∂_ν)(BraKet[ψ, (∂_μ)ψ]))
=-i ((BraKet[(∂_μ)ψ, (∂_ν)ψ])+(BraKet[ψ, (∂_μ)(∂_ν)ψ])-(BraKet[(∂_ν)ψ, (∂_μ)ψ])-(BraKet[ψ, (∂_ν)(∂_μ)ψ]))
=-i((BraKet[(∂_μ)ψ, (∂_ν)ψ])-(BraKet[(∂_ν)ψ, (∂_μ)ψ])).

(Subscript[ℱ, μ,ν])=-i Tr ρ ((∂_μ)ρ(∂_ν)ρ-(∂_ν)ρ(∂_μ)ρ).

(Subscript[ρ, k,±])=(Ket[k,±])(Bra[k,±])=(1/2)(𝟙±(OverHat[h])(k)·σ),

(Subscript[ℱ, k,±])=-i (Tr (Subscript[ρ, k,±])(∂_(Subscript[k, x]))(Subscript[ρ, k,±])(∂_(Subscript[k, y]))(Subscript[ρ, k,±])-(x↔y))
=-i ((1/8)Tr (𝟙±(OverHat[h])·σ)(∂_(Subscript[k, x]))(𝟙±(OverHat[h])·σ)(∂_(Subscript[k, y]))(𝟙±(OverHat[h])·σ)-(x↔y))
=-i ((1/8)Tr (𝟙±(OverHat[h])·σ)((∂_(Subscript[k, x]))(OverHat[h])·σ)((∂_(Subscript[k, y]))(OverHat[h])·σ)-(x↔y))
=-i ((1/8)Tr (𝟙±(OverHat[h])·σ)(((∂_(Subscript[k, x]))(OverHat[h])·(∂_(Subscript[k, y]))(OverHat[h]))𝟙+i ((∂_(Subscript[k, x]))(OverHat[h])×(∂_(Subscript[k, y]))(OverHat[h]))·σ)-(x↔y))
=-i ((1/8)(((∂_(Subscript[k, x]))(OverHat[h])·(∂_(Subscript[k, y]))(OverHat[h]))±i(OverHat[h])·((∂_(Subscript[k, x]))(OverHat[h])×(∂_(Subscript[k, y]))(OverHat[h])))Tr 𝟙 -(x↔y))
=-i ((1/4)(((∂_(Subscript[k, x]))(OverHat[h])·(∂_(Subscript[k, y]))(OverHat[h]))±i(OverHat[h])·((∂_(Subscript[k, x]))(OverHat[h])×(∂_(Subscript[k, y]))(OverHat[h])))-(x↔y))
=-i (1/2)(±i(OverHat[h])·((∂_(Subscript[k, x]))(OverHat[h])×(∂_(Subscript[k, y]))(OverHat[h])))
=± (1/2)(OverHat[h])·((∂_(Subscript[k, x]))(OverHat[h])×(∂_(Subscript[k, y]))(OverHat[h])).

∂(OverHat[h])=∂((h/(|h|)))=((∂h)/(|h|))+((h·∂h)/((|h|)^3))h,

(OverHat[h])·((∂_(Subscript[k, x]))(OverHat[h])×(∂_(Subscript[k, y]))(OverHat[h]))=(h/(|h|))·((((∂_(Subscript[k, x]))h)/(|h|))+((h·(∂_(Subscript[k, x]))h)/((|h|)^3))h)×((((∂_(Subscript[k, y]))h)/(|h|))+((h·(∂_(Subscript[k, y]))h)/((|h|)^3))h)
=(1/((|h|)^3))h·((∂_(Subscript[k, x]))h×(∂_(Subscript[k, y]))h)+((h·(∂_(Subscript[k, y]))h)/((|h|)^5))h·((∂_(Subscript[k, x]))h×h)+((h·(∂_(Subscript[k, x]))h)/((|h|)^5))h·(h×(∂_(Subscript[k, y]))h)+(((h·(∂_(Subscript[k, x]))h)(h·(∂_(Subscript[k, y]))h))/((|h|)^7))h·(h×h)
=(1/((|h|)^3))h·((∂_(Subscript[k, x]))h×(∂_(Subscript[k, y]))h).

(Subscript[ℱ, k,±])=± (1/(2((|h|)^3)))h·((∂_(Subscript[k, x]))h×(∂_(Subscript[k, y]))h).

Apply Eq. () to the Haldane model, one can show the distribution of Berry curvature in the Brillouin zone (for both upper and lower bands) as a function of λ

Block[{λ=0.1,h,F,plt,bg,bar},
h[{kx_,ky_}]:={-(Cos[ky]+2Cos[ky/2]Cos[Sqrt[3]kx/2]),-(-Sin[ky]+2Sin[ky/2]Cos[Sqrt[3]kx/2]),-2λ(Sin[Sqrt[3]kx]-2Sin[Sqrt[3]kx/2]Cos[3ky/2])};
F=Chop@Simplify[Det[Prepend[Transpose[Grad[h[{kx,ky}],{kx,ky}]],h[{kx,ky}]]]/(2Sqrt[h[{kx,ky}].h[{kx,ky}]]^3)];
plt=Association@Table[s->First@ContourPlot[s F,{kx,-4π/(3Sqrt[3]),4π/(3Sqrt[3])},{ky,-2π/3,2π/3},PlotRange->All,BoundaryStyle->Black,ContourStyle->Directive[AbsoluteThickness[1],Opacity[0.2]],
ColorFunction->Function[e,ColorData["BlueRedSplit"][Rescale[e,{-3,3}]]],ColorFunctionScaling->False,RegionFunction->Function[{kx,ky,e},Abs[kx]/(4π/(3Sqrt[3]))+Abs[ky]/(4π/3)<=1]],{s,{1,-1}}];
bg={#1,Point[#3],Text[Style[#2,12],#3,#4]}&@@@{{Black,"Γ",{0,0},{1,1}},{Black,"K",{4π/(3Sqrt[3]),0},{-1.5,0}},{Black,"K'",{2π/(3Sqrt[3]),2π/3},{-1.2,-0.7}},{Black,"M",{π/Sqrt[3],π/3},{-1.2,-0.7}}};
bar=ContourPlot[y,{x,0,1},{y,-3,3},Contours->11,ContourStyle->Directive[AbsoluteThickness[1],Opacity[0.2]],AspectRatio->8,PlotRangePadding->None,ImageMargins->2,FrameTicks->{{None,Range[-3,3,1]},{None,None}},PlotLabel->"\!\(ℱ\_k\)",ImageSize->{Automatic,120}];
Row[{Column[Prepend[#,Style[StringTemplate["\!\(λ/t=`1`\)"][-λ],"InlineFormula"]],Alignment->Center]&@Table[Graphics[{plt[-s],bg,
Text["\!\(ℱ\_\(k,"<>Switch[s,1,"+",-1,"-"]<>"\)\)",{0,2π/3},{0,-1}]},ImageSize->120],{s,{1,-1}}],
Column[Prepend[#,Style[StringTemplate["\!\(λ/t=`1`\)"][λ],"InlineFormula"]],Alignment->Center]&@Table[Graphics[{plt[s],bg,
Text["\!\(ℱ\_\(k,"<>Switch[s,1,"+",-1,"-"]<>"\)\)",{0,2π/3},{0,-1}]},ImageSize->120],{s,{1,-1}}],
bar}]]

What is the physical effect of magnetic field in the momentum space? Recall: magnetic field B in the real space → electron will experience a Lorentz force

((dk)/(dt))∝-v×B=-((dx)/(dt))×B.

By analogy, switching from real space to momentum space effectively exchanges coordinate and momentum,

((dx)/(dt))=-((dk)/(dt))×ℱ.

In the presence of an external electric field E

((ℏdk)/(dt))=-e E,

this will contribute to the current J by

J=-e ((dx)/(dt))=((e^2)/ℏ)ℱ×E.

So ((e^2)/ℏ)(Subscript[ℱ, k,α]) is the contribution to the Hall conductance from a single state (Ket[k,α]).

Chern number: integration of Berry curvature over the first Brillouin zone (in unit of 2π), like the total magnetic flux through the Brillouin zone,

The Chern number (Subscript[𝒞, ±]) is separately define for each band (labeled by ±). For two-band system, Chern numbers are opposite between upper and lower bands, (SubPlus[𝒞])=-(SubMinus[𝒞]).

It turns out that (Subscript[𝒞, ±])∈ℤ is always an integer. For Haldane model,

(Subscript[𝒞, ±])=± sgn λ.

As λ is tuned across λ=0, the Chern number exchanges between the upper and lower bands.

Plot[{-Sign[λ],Sign[λ]},{λ,-0.1,0.1},FrameLabel->{"\!\(λ\)","\!\(𝒞\)"},Epilog->{{Red,Text["\!\(𝒞\_+\)",{-0.05,-1},{0,-1}],Blue,Text["\!\(𝒞\_-\)",{-0.05,1},{0,1}]},{Red,Text["\!\(𝒞\_+\)",{0.05,1},{0,1}],Blue,Text["\!\(𝒞\_-\)",{0.05,-1},{0,-1}]}},ImageSize->200]

The Chern number is a topological number of the vector bundle (Ket[k,±]). The 2D insulators with non-trivial Chern numbers in occupied bands are called Chern insulators. The jump of Chern number at λ=0 signifies a topological transition between two distinct Chern insulators.

The Hall conductance of the Chern insulator (at the low-temperature limit) equals the total the Chern numbers (Subscript[𝒞, α]) of all occupied bands (bands below the Fermi level) times the quantum conductance (e^2)/h.

Argument: the Hall conductance of a band averages the contribution from all states in the band

(Subscript[σ, H])=((e^2)/ℏ)∫(((d^2)k)/(((2π))^2))(Subscript[ℱ, k,α])=((e^2)/h)(1/(2π))∫(d^2)k (Subscript[ℱ, k,α])=((e^2)/h)(Subscript[𝒞, α]).

The phenomenon that the Hall conductance quantize to integer multiplies of (e^2)/h is called the integer quantum Hall effect.

(Subscript[H, K+κ])≏(Subscript[κ, x])(σ^1)+(Subscript[κ, y])(σ^2)+m (σ^3),
(Subscript[H, K'+κ])≏-(Subscript[κ, x])(σ^1)+(Subscript[κ, y])(σ^2)-m (σ^3),

Using Eq. (), the Berry curvature can be evaluated

(Subscript[ℱ, K+κ,±])=(Subscript[ℱ, K'+κ,±])=±(m/(2((((κ^2)+(m^2)))^(3/2))))

Each massive Dirac fermion contribute ±1/2 to the band Chern number. For example, around K point

(Subscript[𝒞, K,±])=(1/(2π))∫(d^2)κ (Subscript[ℱ, K+κ,±])=±(1/2)(m/(|m|))=±(1/2)sgn m.

Same result will be obtained around K', (Subscript[𝒞, K,±])=(Subscript[𝒞, K',±])=±(1/2)sgn m, such that the Chern number follows

(Subscript[𝒞, ±])=(Subscript[𝒞, K,±])+(Subscript[𝒞, K',±])=± sgn m,

confirming the observation in Eq. ().

Classifying Topological Insulators in 2D

H≏(Underscript[⊕, k])(Subscript[H, k]),(Subscript[H, k])≏h(k)·σ.

The band gap is set by |h(k)|. To keep the gap open, must have

∀k:\[LeftBracketingBar]h(k)\[RightBracketingBar]>0,

meaning that h(k)∈(ℝ^3)\{0}. For 2D system, k∈(T^2) (2D torus).

All smooth deformations of the function h:(T^2)→(ℝ^3)\{0} are classified by the homotopy group

which is the integer addition group. Each class is labeled by an integer topological index ω, which corresponds to the winding number of h in the 2D Brillouin zone

ω=(1/(4π))(∫_BZ)(d^2)k (1/((|h(k)|)^3))h(k)·((∂_(Subscript[k, x]))h(k)×(∂_(Subscript[k, y]))h(k)).

By Eq. (), the winding number is the same as the Chern number. Thus topological insulators (with no additional symmetry other than U(1)) is classified by the band Chern number.

In the momentum space, the Hamiltonian still takes the form of

H≏(Underscript[⊕, k])(Subscript[H, k]),(Subscript[H, k])≏(({{(Subscript[g, k]), (Subscript[f, k]^( * ))}, {(Subscript[f, k]), -(Subscript[g, k])}})),

but the function (Subscript[g, k]) is updated to

(Subscript[g, k])=δ-2λ (sin k·((Subscript[δ, 1])-(Subscript[δ, 2]))+sin k·((Subscript[δ, 2])-(Subscript[δ, 3]))+sin k·((Subscript[δ, 3])-(Subscript[δ, 1]))).

Expanding around K and K',

(Subscript[H, K+κ])≏(Subscript[v, F])((Subscript[κ, x])(σ^1)+(Subscript[κ, y])(σ^2))+(Subscript[m, K]) (Subscript[v, F]^2)(σ^3)+𝒪((κ^2)),
(Subscript[H, K'+κ])≏(Subscript[v, F])(-(Subscript[κ, x])(σ^1)+(Subscript[κ, y])(σ^2))-(Subscript[m, K']) (Subscript[v, F]^2)(σ^3)+𝒪((κ^2)),

the Dirac masses are different

(Subscript[m, K])=(1/(Subscript[v, F]^2))(3(Sqrt[3])λ+δ),
(Subscript[m, K'])=(1/(Subscript[v, F]^2))(3(Sqrt[3])λ-δ).

The band Chern numbers

(Subscript[𝒞, ±])=(Subscript[𝒞, K,±])+(Subscript[𝒞, K',±])=± (1/2)(sgn (Subscript[m, K])+sgn (Subscript[m, K'])).

If λ=0 and δ≠0, the masses are opposite between K and K' valleys => the Chern number from both valleys cancels, i.e. (Subscript[𝒞, ±])=0 => the insulator is in the trivial phase.

Chiral Edge Mode

periodic boundary condition along x-direction (preserving translation symmetry along x-direction, (Subscript[k, x]) is still well-defined),

open boundary condition along y-direction (top and bottom edges breaks translation symmetry along y-direction, can not define (Subscript[k, y])),

Block[{Lx=3,Ly=4,a1,a2,row},
a1={Sqrt[3],0};
a2={Sqrt[3],3}/2;
row[dir_]:={{Line@Table[{Sqrt[3]x/2,-dir Cos[π x]/4},{x,-Lx-1/2,Lx+1/2,1/2}]},
{Arrowheads[{{0.07,0.65}}],LightGray,Table[Arrow[{{dir Sqrt[3](x-2)/2,-dir/4},{dir Sqrt[3]x/2,-dir/4}}],{x,-Lx+1,Lx+1,2}],
Table[Arrow[{{dir Sqrt[3](x+2)/2,dir/4},{dir Sqrt[3]x/2,dir/4}}],{x,-Lx-2,Lx,2}]},
{EdgeForm[Black],Table[{Lighter[If[dir(-1)^x>0,Blue,Red],0.6],Disk[{Sqrt[3]x/2,-dir Cos[π x]/4},0.15]},{x,-Lx,Lx}]}};
Graphics[{
Table[If[OddQ[x+y],{Translate[Line[{0,#}&/@{0,1}],{Sqrt[3]x/2,3y/2+1/4}]},{Arrowheads[{{0.07,0.65}}],LightGray,Translate[Rotate[Arrow[{{Sqrt[3]/2,1/2},{0,-1}}],2π/6 #,{0,0}]&/@{0,1,3,4},{Sqrt[3]x/2,3y/2+3/4}]}],{y,Ly-1},{x,-Lx-1,Lx+1}],
Table[Translate[row[(-1)^y],{0,3y/2}],{y,Ly}],
{Arrowheads[0.06],Translate[{Arrow[{{0,0},#1}],Text[StringTemplate["\!\(\*StyleBox[\"a\",Bold]\_`1`\)"][#2],#1,#3]}&@@@{{a1,"1",{0,0.8}},{a2,"2",{-1.,-0.2}}},{0,3+3/4}]}},PlotRange->{{-Lx-1/2,Lx+1/2}Sqrt[3]/2,All},
ImageSize->140]]

(Ket[(Subscript[k, x]),(Subscript[y, i]),s])=(∑_(Subscript[x, i]))(e^(-i (Subscript[k, x]) (Subscript[x, i])))(Ket[i]).

The Hamiltonian is diagonal in (Subscript[k, x]) (because (Subscript[k, x]) is a good quantum number)

H=(∑_(Subscript[k, x]))(Ket[(Subscript[k, x])])H((Subscript[k, x]))(Bra[(Subscript[k, x])]),

where

H((Subscript[k, x]))=(Subscript[H, t])((Subscript[k, x]))+(Subscript[H, λ])((Subscript[k, x]))+(Subscript[H, δ])((Subscript[k, x]))
(Subscript[H, t])((Subscript[k, x]))=-t(∑_y)((Ket[y,B])(Bra[y,A])+2cos(((Sqrt[3])(Subscript[k, x]))/2) (Ket[y+(Subscript[a, y]),A])(Bra[y,B])+h.c.)
(Subscript[H, λ])((Subscript[k, x]))=-2λ (∑_((Subscript[k, x]),y))(sin (Sqrt[3])(Subscript[k, x]) ((Ket[y,A])(Bra[y,A])-(Ket[y,B])(Bra[y,B]))-sin (((Sqrt[3])(Subscript[k, x]))/2) ((Ket[y+(Subscript[a, y]),A])(Bra[y,A])-(Ket[y+(Subscript[a, y]),B])(Bra[y,B])+h.c.))
(Subscript[H, δ])=δ(∑_((Subscript[k, x]),y))((Ket[y,A])(Bra[y,A])-(Ket[y,B])(Bra[y,B])),

with (Subscript[a, y])=3/2 being the inter-layer spacing.

Further diagonalize H((Subscript[k, x])) for every (Subscript[k, x])

H((Subscript[k, x]))(Ket[(Subscript[u, n])((Subscript[k, x]))])=(Subscript[ϵ, n])((Subscript[k, x]))(Ket[(Subscript[u, n])((Subscript[k, x]))])

Bulk modes (in black) and edge modes (in red)

Manipulate[Block[{L=20,t=1,H,Es},
H[kx_]:=KroneckerProduct[SparseArray[Band[{1,1}]->1,{L,L}],{{δ-2λ Sin[Sqrt[3] kx],-t},{-t,-δ+2λ Sin[Sqrt[3] kx]}}]+
(#+Transpose[#])&@KroneckerProduct[SparseArray[Band[{2,1}]->1,{L,L}],{{2λ Sin[Sqrt[3]kx/2],-2t Cos[Sqrt[3]kx/2]},{0,-2λ Sin[Sqrt[3]kx/2]}}];
Es[kx_]:=Join[#[[;;L]],Reverse@#[[L+1;;]]]&@Sort@Eigenvalues@N@H[kx];
ListLinePlot[Thread@Table[{kx,#}&/@Es[kx],{kx,π/Sqrt[3]Range[-1,1,1/32]}],PlotStyle->Append[ConstantArray[Black,L-1],Red],FrameLabel->{"\!\(k\_x\)","\!\(ϵ(k\_x)/t\)"},FrameTicks->{{Automatic,Automatic},{π/Sqrt[3]Range[-1,1],None}},PlotLabel->StringTemplate["\!\(λ/t=`1`,δ/t=`2`\)"][Round[λ,0.01],Round[δ,0.01]],ImageSize->220]
],{{λ,0},-0.1,0.1,ImageSize->Small},{{δ,0},-0.5,0.5,ImageSize->Small},Paneled->False]

Edge mode wave function (y-dependence). Edge modes are confined on the boundary. The wave functions decay exponentially into the bulk.

Block[{L=20,t=1,λ=0.1,δ=0.,k=0.,H,S,us},
H[kx_]:=KroneckerProduct[SparseArray[Band[{1,1}]->1,{L,L}],{{δ-2λ Sin[Sqrt[3] kx],-t},{-t,-δ+2λ Sin[Sqrt[3] kx]}}]+
(#+Transpose[#])&@KroneckerProduct[SparseArray[Band[{2,1}]->1,{L,L}],{{2λ Sin[Sqrt[3]kx/2],-2t Cos[Sqrt[3]kx/2]},{0,-2λ Sin[Sqrt[3]kx/2]}}];
S=KroneckerProduct[SparseArray[Band[{1,1}]->1,{L,L}],{{1,0},{0,-1}}];
us=Cases[Thread@Eigensystem@(H[k]+0.01S),{e_,u_}/;Abs[e]<Sqrt[9λ^2+δ^2]:>{u,u.S.u}];
Column[Join[{Style[StringTemplate["          \!\(k\_x=`1`\)"][k],"InlineFormula"]},Graphics[{If[EvenQ[#1],Blue,Red],Point[{#1,#2}],Opacity[0.5],Line[{{#1,#2},{#1,0}}]}&@@@Thread@{Range[2L],#1},PlotRange->{{0,2L+1},{-1,1}},Frame->True,FrameLabel->{None,"\!\(u\_\(i\)\%"<>If[#2>0,"btm","top"]<>"\)"},AspectRatio->0.4,ImageSize->220]&@@@us,{Style["          \!\(y\_i/a\_y\)","InlineFormula"]}],Alignment->Center]
]

Color every point ((Subscript[k, x]),(Subscript[ϵ, n])((Subscript[k, x]))) by the weight of the corresponding mode (Ket[(Subscript[u, n])((Subscript[k, x]))]) on a given site. Define the spectral function (Subscript[A, y,s]) on layer y sublattice s,

(Subscript[A, y,s])((Subscript[k, x]),ϵ)=(∑_n)((|(BraKet[y,s, (Subscript[u, n])((Subscript[k, x]))])|)^2)δ(ϵ-(Subscript[ϵ, n])((Subscript[k, x]))).

Bottom boundary (y=1,s=A)

Block[{L=20,t=1,λ=0.1,δ=0.,H,S,P,us},
H[kx_]:=KroneckerProduct[SparseArray[Band[{1,1}]->1,{L,L}],{{δ-2λ Sin[Sqrt[3] kx],-t},{-t,-δ+2λ Sin[Sqrt[3] kx]}}]+
(#+Transpose[#])&@KroneckerProduct[SparseArray[Band[{2,1}]->1,{L,L}],{{2λ Sin[Sqrt[3]kx/2],-2t Cos[Sqrt[3]kx/2]},{0,-2λ Sin[Sqrt[3]kx/2]}}];
S=KroneckerProduct[SparseArray[Band[{1,1}]->1,{L,L}],{{1,0},{0,-1}}];
P[i_,s_]:=SparseArray[{2(i-1)+s,2(i-1)+s}->1,{2L,2L}];
us[i_,s_]:=Table[Cases[Thread@Eigensystem@(H[kx]+0.01S),{e_,u_}:>{Opacity[u.P[i,s].u],Point[{kx,e}]}],{kx,π/Sqrt[3]Range[-1,1,0.02]}];
Graphics[us[1,1],Frame->True,FrameLabel->{"\!\(k\_x\)","\!\(ϵ(k\_x)/t\)"},FrameTicks->{{Automatic,Automatic},{π/Sqrt[3]Range[-1,1],None}},PlotLabel->Style[StringTemplate["\!\(λ/t=`1`,δ/t=`2`\)"][Round[λ,0.01],Round[δ,0.01]],14,Black],AspectRatio->1/GoldenRatio,ImageSize->220]]

Effective dispersion of boundary modes: ϵ((Subscript[k, x]))≃-v (Subscript[k, x]) => Left-moving modes.

Edge modes merge into bulk modes as the dispersion touches the bulk band.

Top boundary (y=(Subscript[L, y]),s=B)

Block[{L=20,t=1,λ=0.1,δ=0.,H,S,P,us},
H[kx_]:=KroneckerProduct[SparseArray[Band[{1,1}]->1,{L,L}],{{δ-2λ Sin[Sqrt[3] kx],-t},{-t,-δ+2λ Sin[Sqrt[3] kx]}}]+
(#+Transpose[#])&@KroneckerProduct[SparseArray[Band[{2,1}]->1,{L,L}],{{2λ Sin[Sqrt[3]kx/2],-2t Cos[Sqrt[3]kx/2]},{0,-2λ Sin[Sqrt[3]kx/2]}}];
S=KroneckerProduct[SparseArray[Band[{1,1}]->1,{L,L}],{{1,0},{0,-1}}];
P[i_,s_]:=SparseArray[{2(i-1)+s,2(i-1)+s}->1,{2L,2L}];
us[i_,s_]:=Table[Cases[Thread@Eigensystem@(H[kx]+0.01S),{e_,u_}:>{Opacity[u.P[i,s].u],Point[{kx,e}]}],{kx,π/Sqrt[3]Range[-1,1,0.02]}];
Graphics[us[L,2],Frame->True,FrameLabel->{"\!\(k\_x\)","\!\(ϵ(k\_x)/t\)"},FrameTicks->{{Automatic,Automatic},{π/Sqrt[3]Range[-1,1],None}},PlotLabel->Style[StringTemplate["\!\(λ/t=`1`,δ/t=`2`\)"][Round[λ,0.01],Round[δ,0.01]],14,Black],AspectRatio->1/GoldenRatio,ImageSize->220]]

Effective dispersion of boundary modes: ϵ((Subscript[k, x]))≃+v (Subscript[k, x]) => Right-moving modes.

Edge modes merge into bulk modes as the dispersion touches the bulk band.

Chiral edge mode: The edge mode circulates around the material in a chiral manner.

Block[{spec},
spec[kF_,s_,lbl_]:=Block[{ϵ},ϵ[k_]:=0.33Tanh[4 s k];{{LightGray,Rectangle[{-0.5,0.3},{0.5,0.5}],Line[{#,ϵ[#]}&/@Range[-0.46,0.46,0.01]]},{Black,Rectangle[{-0.5,-0.3},{0.5,-0.5}],Line[{#,ϵ[#]}&/@Range[-s 0.46,kF,s 0.01]]},{AbsoluteThickness[1],Dashed,Line[{#,ϵ[kF]}&/@{-0.5,0.5}]},Text["\!\(μ\_"<>lbl<>"\)",{0.5,ϵ[kF]},{-1,0}]}];Graphics[{{LightGray,EdgeForm[None],Rectangle[{-2,-1},{2,1}]},{Line[{#,1}&/@{-2,2}],Line[{#,-1}&/@{-2,2}]},{Red,Arrowheads[{{0.06,0.1},{0.06,0.5},{0.06,0.9}}],Arrow[{#,0.85}&/@{-2,2}],Arrow[{#,-0.85}&/@{2,-2}]},{Translate[spec[0,1,"top"],{2.8,1}],Translate[spec[0,-1,"btm"],{2.8,-1}]},
{Text["top edge",{0,1},{0,-1}],Text["bottom edge",{0,-1},{0,1}]}},ImageSize->220]]

Quantized Hall conductance. An electric field (Subscript[E, y]) creates electrostatic potential difference (Subscript[V, y]) between top and bottom boundaries => difference in electron chemical potentials (Subscript[μ, top]) and (Subscript[μ, btm])

-e (Subscript[E, y])(Subscript[L, y])=-e (Subscript[V, y])=(Subscript[μ, top])-(Subscript[μ, btm]),

Block[{spec},
spec[kF_,s_,lbl_]:=Block[{ϵ},ϵ[k_]:=0.33Tanh[4 s k];{{LightGray,Rectangle[{-0.5,0.3},{0.5,0.5}],Line[{#,ϵ[#]}&/@Range[-0.46,0.46,0.01]]},{Black,Rectangle[{-0.5,-0.3},{0.5,-0.5}],Line[{#,ϵ[#]}&/@Range[-s 0.46,kF,s 0.01]]},{AbsoluteThickness[1],Dashed,Line[{#,ϵ[kF]}&/@{-0.5,0.5}]},Text["\!\(μ\_"<>lbl<>"\)",{0.5,ϵ[kF]},{-1,0}]}];Graphics[{{LightGray,EdgeForm[None],Rectangle[{-2,-1},{2,1}]},{Line[{#,1}&/@{-2,2}],Line[{#,-1}&/@{-2,2}]},{Red,Arrowheads[{{0.06,0.1},{0.06,0.5},{0.06,0.9}}],AbsoluteThickness[3],Arrow[{#,0.85}&/@{-2,2}],AbsoluteThickness[1],Arrow[{#,-0.85}&/@{2,-2}]},{Translate[spec[0.1,1,"top"],{2.8,1}],Translate[spec[0.1,-1,"btm"],{2.8,-1}]},
{Text["top edge",{0,1},{0,-1}],Text["bottom edge",{0,-1},{0,1}]},{Arrow[{0,#}&/@{0.5,-0.5}],Text["\!\(E\_y\)",{0,-0.5},{-1.5,-0.5}],Red,Arrow[{#,0}&/@{-0.5,-1.5}],Text["\!\(j\_x\)",{-1,0},{0,-1}]}},ImageSize->220]]

more electrons occupy the top edge modes compared to bottom edge modes

(Subscript[N, top])-(Subscript[N, btm])=(((Subscript[k, x,top])-(Subscript[k, x,btm]))/((2π/(Subscript[L, x]))))=(((Subscript[μ, top])-(Subscript[μ, btm]))/(v ℏ (2π/(Subscript[L, x]))))=((-e (Subscript[E, y])(Subscript[L, x])(Subscript[L, y]))/(v h)),

more right-moving electrons than left-moving electrons contributes to a net current

(Subscript[I, x])=(Subscript[j, x]) (Subscript[L, y])=-e v (((Subscript[N, top])-(Subscript[N, btm]))/(Subscript[L, x]))=(((e^2) (Subscript[E, y])(Subscript[L, y]))/h).

Hall conductance (Subscript[σ, H]) is defined by

(Subscript[j, x])=(Subscript[σ, H]) (Subscript[E, y]),

Eq. () implies

(Subscript[σ, H])=((e^2)/h).

which is indeed quantized (integer) in unit of the quantum conductance (e^2)/h.

Bulk-boundary correspondence. The number of chiral edge modes (boundary property) ⇔ the quantized Hall conductance / the band Chern number (bulk property).

General Theory of Electrons in Crystals

Nearly Free Electron Model

(Ket[k]): plane wave state of electron (Subscript[ψ, k])(x)∼(e^(i k·x)), labeled by the wave vector k.

(Ket[k])=∫dx (Subscript[ψ, k])(x)(Ket[x])=(1/(V^(1/2)))∫dx (e^(i k·x))(Ket[x]).

V - volume of the system, to provide appropriate normalization for the wave function.

A key feature is that V(x) is periodic

V(x)=V(x+R),

where R is any lattice vector. Because atoms are arranged on a periodic lattice in the crystal.

Transform the potential energy to the momentum space

∫dx(Ket[x])V(x)(Bra[x])=(∑_(k,k'))(Ket[k'])(Subscript[V, k'-k])(Bra[k]),

with (Subscript[V, q]) being the Fourier transform of V(x)

(Subscript[V, q])=(1/V)∫dx V(x)(e^(-i q·x)).

Derive Eq. ().

∫dx(Ket[x])V(x)(Bra[x])=(∑_(k,k'))∫dx(Ket[k'])(BraKet[k', x])V(x)(BraKet[x, k])(Bra[k])
=(∑_(k,k'))(Ket[k'])((1/V)∫dx (e^(-i k'·x))V(x)(e^(i k·x)))(Bra[k])
=(∑_(k,k'))(Ket[k'])((1/V)∫dx V(x)(e^(-i (k'-k)·x)))(Bra[k]).

However, the periodicity of V(x) in Eq. () puts a strong constraint on its Fourier transform (Subscript[V, q]) in Eq. () that

(Subscript[V, q])=(e^(-i q·R))(Subscript[V, q]),

for any lattice vector R. For this equation to hold,

Either q satisfies the condition ∀R:(e^(-i q·R))=1, i.e. q=G∈ reciprocal lattice.

Or q is not on the reciprocal lattice, then (Subscript[V, q])=0 must vanish.

(Subscript[V, q])=(∑_G)(Subscript[V, G])δ(q-G).

G is summed over all reciprocal vectors on the reciprocal lattice,

(Subscript[V, G]) describes the scattering strength that an electron is scattered from momentum ℏ k to ℏ(k+G) on the lattice, and is given by

(Subscript[V, G])=(1/V)∫dx V(x)(e^(-i G·x)).

The (Ket[k]) state can only be mixed with (Ket[k+G]) states related by the lattice momentum. This is a manifestation of the Laue condition (the conservation of quasi-momentum).

Electron in 1D Periodic Potential

V(x)=(Subscript[V, 0])+2(Subscript[V, 1])cos(2π x)+2(Subscript[V, 2])cos(4π x).

Plot[-8Cos[2π x]-8Cos[4π x],{x,-2.5,2.5},PlotStyle->Black,Axes->False,AspectRatio->1/2,FrameLabel->{"\!\( x\)","\!\(V(x)\)"},ImageSize->200]

Lattice constant a=1 => Reciprocal lattice constant b=2π/a=2π => Reciprocal lattice vectors (Subscript[G, m])=m b=2π m (m∈ℤ - the Miller index, not the electron mass)

Scattering strength

(Subscript[V, G])=({{{, ({{(Subscript[V, 0]), G=0}, {(Subscript[V, 1]), G=±2π}, {(Subscript[V, 2]), G=±4π}, {0, otherwise}})}}),

Takes the set of basis states given a momentum k∈BZ in the first Brillouin zone

…,(Ket[k-4π]),(Ket[k-2π]),(Ket[k]),(Ket[k+2π]),(Ket[k+4π]),…

the Hamiltonian can be represented as a matrix

(Subscript[H, k])=(({{⋱, ⋱, ⋱,  ,  ,  ,  }, {⋱, (Subscript[ξ, k-4π]), (Subscript[V, 1]), (Subscript[V, 2]),  ,  ,  }, {⋱, (Subscript[V, 1]), (Subscript[ξ, k-2π]), (Subscript[V, 1]), (Subscript[V, 2]),  ,  }, { , (Subscript[V, 2]), (Subscript[V, 1]), (Subscript[ξ, k]), (Subscript[V, 1]), (Subscript[V, 2]),  }, { ,  , (Subscript[V, 2]), (Subscript[V, 1]), (Subscript[ξ, k+2π]), (Subscript[V, 1]), ⋱}, { ,  ,  , (Subscript[V, 2]), (Subscript[V, 1]), (Subscript[ξ, k+4π]), ⋱}, { ,  ,  ,  , ⋱, ⋱, ⋱}})),

where the diagonal energy is given by

(Subscript[ξ, k])=(((ℏ^2)(k^2))/(2m))+(Subscript[V, 0]).

Diagonalize the Hamiltonian to fine the eigen energies and eigenstates

(Subscript[H, k])(Ket[α,k])=(Subscript[ϵ, α,k])(Ket[α,k]).

α - band index, labels the nth eigen state of (Subscript[H, k]) at each k

k - quasi-momentum, there is no need to take k out side the first Brillouin zone as that do not leads to a new set of basis states.

(Subscript[ϵ, α,k]) - dispersion relation of the αth energy band (corresponds to the nth eigenvalue of (Subscript[H, k]))

Column@Table[Block[{M=5,m=1,V0=0,V1=-4V,V2=-4V,H},
H[k_]:=SparseArray[Join[{Band[{1,1}]->(k+2π Range[-M,M])^2/(2m)+V0,Band[{1,2}]->V1,Band[{2,1}]->V1,Band[{1,3}]->V2,Band[{3,1}]->V2}],{#,#}&@(2M+1)];
Row@{Rotate[Style[StringTemplate["\!\(V(x)`1`0\)"][<|0->"=",1->"≠"|>[V]],14,SingleLetterItalics->True,FontFamily->"LaTeX"],π/2],Graphics[{Line/@Thread@Table[{k,#}&/@Sort[Eigenvalues@H[k]][[;;3]],{k,-π,π,0.01π}]},PlotRange->{-5,40},PlotRangeClipping->True,Frame->True,FrameLabel->{"\!\( k\)","\!\(ϵ\_\(α,k\)\)"},FrameTicks->{{Automatic,Automatic},{Range[-3π,3π,π],None}},PlotLabel->Style["Folded",14,Black],AspectRatio->2,ImageSize->{Automatic,180}],Graphics[{{LightGray,Dashed,Line/@Thread@Table[{k,#}&/@Sort[Eigenvalues@H[k]][[;;3]],{k,-3π,3π,0.01π}]},Table[Line[SortBy[#,First]]&/@Thread@Table[MapIndexed[Block[{n},{n}=#2;{ k- 2π (-1)^n Floor[n/2]Sign[k],#1}]&,Sort[Eigenvalues@H[k]][[;;3]]],{k,s 0.001π,s π,s 0.01π}],{s,{-1,1}}]},GridLines->{DeleteCases[Range[-3π,3π,π],0],None},PlotRange->{-5,40},PlotRangeClipping->True,Frame->True,FrameLabel->{"\!\( k\)","\!\(ϵ(k)\)"},FrameTicks->{{Automatic,Automatic},{Range[-3π,3π,π],None}},PlotLabel->Style["Extended",14,Black],AspectRatio->2/3,ImageSize->{Automatic,180}]}],{V,0,1}]

Band gap opens at the Brillouin zone boundary.

The dispersion can be unfolded to extended Brillouin zones following (Subscript[ϵ, α,k])->ϵ(k-2π (((-))^α)⌊(α/2)⌋sgn(k)).

(Ket[α,k]) - eigenstate, as a linear combination of basis states

(Ket[α,k])=(∑_(m∈ℤ))(Subscript[(OverTilde[u]), α, k, 2π m])(Ket[k+2π m]),

in real space basis

(Subscript[ψ, α,k])(x)=(BraKet[x, α,k])=(∑_(m∈ℤ))(Subscript[(OverTilde[u]), α,k,2π m])(e^(i (k+2π m)x)),

Manipulate[Block[{M=5,m=1,V0=0,V1=-4,V2=-4,H,xs,ψ},
H[k_]:=SparseArray[Join[{Band[{1,1}]->(k+2π Range[-M,M])^2/(2m)+V0,Band[{1,2}]->V1,Band[{2,1}]->V1,Band[{1,3}]->V2,Band[{3,1}]->V2}],{#,#}&@(2M+1)];
xs=Range[-5.5,5.5,0.05];
ψ[n_,k_]:=Table[Exp[I(k+2π Range[-M,M]) x],{x,xs}].(#/Sign[Total[#]])&[SortBy[Thread@Eigensystem@N@H[k],First][[n,-1]]];
Row[{ListLinePlot[Thread[{xs,#}]&/@Thread@ReIm@ψ[n,k],PlotRange->{-2,2},BaselinePosition->Scaled@0.64,FrameLabel->{"\!\( x\)","\!\(ψ\_\(α,k\)(x)\)"},PlotLabel->StringTemplate["\!\(n=`1`, k=`2`π\)"][n,Round[k/π,0.01]],AspectRatio->0.4,ImageSize->220],Grid[{Graphics[{#1,Line[{{0,0},{1,0}}]},AspectRatio->0.1,ImageSize->20],Style[#2,FontFamily->"LaTeX"]}&@@@{{Blue,"Re"},{Red,"Im"}}]}]
],{n,Range[3]},{{k,π/5},-π,π,ImageSize->Small},Paneled->False]

Every eigenstate is a plane wave (e^(i k x)) modulated by a periodic function (Subscript[u, α,k])(x) (in response to the periodic potential). (Subscript[u, α,k])(x) is called the Bloch function, defined as

(Subscript[u, α,k])(x)=(∑_(m∈ℤ))(Subscript[(OverTilde[u]), α,k,2π m])(e^(i 2π m x)),

such that Eq. () explicitly takes the form of

(Subscript[ψ, α,k])(x)=(Subscript[u, α,k])(x)(e^(i k x)).

Manipulate[Block[{M=5,m=1,V0=0,V1=-4,V2=-4,H,xs,u},
H[k_]:=SparseArray[Join[{Band[{1,1}]->(k+2π Range[-M,M])^2/(2m)+V0,Band[{1,2}]->V1,Band[{2,1}]->V1,Band[{1,3}]->V2,Band[{3,1}]->V2}],{#,#}&@(2M+1)];
xs=Range[-5.5,5.5,0.05];
u[n_,k_]:=Table[Exp[I(2π Range[-M,M]) x],{x,xs}].SortBy[Thread@Eigensystem@N@H[k],First][[n,-1]];
Row[{ListLinePlot[Thread[{xs,#}]&/@Thread@ReIm[#/Sign[Total[#]]]&@u[n,k],PlotRange->{-2,2},BaselinePosition->Scaled@0.64,FrameLabel->{"\!\( x\)","\!\(u\_\(α,k\)(x)\)"},PlotLabel->StringTemplate["\!\(n=`1`, k=`2`π\)"][n,Round[k/π,0.01]],AspectRatio->0.4,ImageSize->220],Grid[{Graphics[{#1,Line[{{0,0},{1,0}}]},AspectRatio->0.1,ImageSize->20],Style[#2,FontFamily->"LaTeX"]}&@@@{{Blue,"Re"},{Red,"Im"}}]}]
],{n,Range[3]},{{k,π/5},-π,π,ImageSize->Small},Paneled->False]

The Bloch function is periodic in x. Because if x is shifted by a lattice vector x→x+n (n∈ℤ), by definition, we have

(Subscript[u, α,k])(x+n)=(∑_(m∈ℤ))(Subscript[(OverTilde[u]), α,k,2π m])(e^(i 2π m x))(e^(i 2π m n))
=(∑_(m∈ℤ))(Subscript[(OverTilde[u]), α,k,2π m])(e^(i 2π m x))=(Subscript[u, α,k])(x).

Bloch Theorem

Proof: in the Hamiltonian Eq. (), because the scattering term (∑_G)(Ket[k+G])(Subscript[V, G])(Bra[k]) only connects momentum eigenstates (Ket[k]) and (Ket[k+G]) that are related by the reciprocal lattice vector G. So the energy eigenstates (eigenstates of H) must be a linear combination of these momentum eigenstates, with some combination coefficients (Subscript[(OverTilde[u]), α,k,G])

(Ket[α,k])=(∑_G)(Subscript[(OverTilde[u]), α,k,G])(Ket[k+G]),

where α is the band index (labeling the αth eigenstate of H). In real space basis,

(BraKet[x, α,k])=(∑_G)(Subscript[(OverTilde[u]), α,k,G])(BraKet[x, k+G])
=(∑_G)(Subscript[(OverTilde[u]), α,k,G])(e^(i (k+G)·x))
=(Subscript[u, α,k])(x)(e^(i k·x)),

where the Bloch function (Subscript[u, α,k])(x) is introduced as the Fourier transform of the coefficients (Subscript[(OverTilde[u]), α,k,G]),

(Subscript[u, α,k])(x)=(∑_G)(Subscript[(OverTilde[u]), α,G])(e^(i G·x)),

which must be periodic in x, as

(Subscript[u, α,k])(x+R)=(∑_G)(Subscript[(OverTilde[u]), α,G])(e^(i G·x))(e^(i G·R))=(Subscript[u, α,k])(x),

given that (e^(i G·R))=1 by definition of the reciprocal lattice.

k is called the quasi-momentum because (Ket[α,k]) state is a superposition of different momentum eigenstates, that does not has a definite momentum, but it labels a equivalent class of momenta that are all related by the Laue condition.

({{(Ket[α,k]), =, (∑_G)(Subscript[(OverTilde[u]), α,k,G]), (Ket[k+G])}, {↑,  ,  , ↑}, {quasi-momentum,  ,  , true momentum}})

The electron still propagate in the form of a (modulated) plane wave in the crystal, with a quasi-momentum, similar to the plane wave state of a free electron.

Even though the potential that the electron feels from each atom is extremely strong, the electrons will still behave almost as if they do not see the atoms at all!

Equation of Motion for Bloch Electron

Classical mechanics (Newton’s 2nd law)

((dp)/(dt))=F=-e E.

Quantum mechanics (Schrödinger equation)

i ℏ (∂_t)ψ=H ψ,

formal solution:

ψ→(e^(-(i/ℏ) H t))ψ,

the electron matter wave will accumulate the phase in time in the rate that is set by the energy.

The potential energy of an electron in the electric field is

V(x)=e E·x.

A plane wave state (Ket[k]) will evolve in the electric field as

(e^(i k·x))→(e^(-(i/ℏ) V(x) t))(e^(i k·x))=(e^(i (k-(1/ℏ)e E t)·x)),

meaning that the momentum ℏ k of the electron will change in time as

ℏ k→ℏ k-e E t,

or described by

((d(ℏ k))/(dt))=-e E,

matching the classical mechanics result.

(Ket[α,k])=(∑_G)(Subscript[(OverTilde[u]), α,k,G])(Ket[k+G]),

(Ket[k])→(Ket[k - (e/ℏ) E t]).

(Ket[α,k])→(∑_G)(Subscript[(OverTilde[u]), α,k,G])(Ket[k - (e/ℏ) E t + G])
≃(∑_G)(Subscript[(OverTilde[u]), α,k - (e/ℏ) E t,G])(Ket[k - (e/ℏ) E t + G])
=(Ket[α,k - (e/ℏ) E t]),

Metal (Fermi Gas) Phase

Thermodynamic Properties

Fermi-Dirac Distribution

Pauli exclusion principle. A fermion mode can either be empty or occupied by one fermion, but not occupied by multiple fermions.

n=({{{, ({{0, empty}, {1, occupied}})}}).

Total energy of the system

E=ϵ n=({{{, ({{0, empty}, {ϵ, occupied}})}}).

Total fermion number of the system

N=n=({{{, ({{0, empty}, {1, occupied}})}}).

Probability to observe the system in the state with occupation number n is

(Subscript[p, n])=(1/Z)(e^(-β (E-μ N)))=(1/Z)(e^(-β (ϵ-μ)n)),

where Z is the partition function (quantum, grand canonical ensemble)

Z=(∑_(n=0,1))(e^(-β (ϵ-μ)n))=1+(e^(-β (ϵ-μ))).

β=1/((Subscript[k, B])T): (Subscript[k, B]) - Boltzmann constant, T - temperature.

μ - chemical potential.

More explicitly

(Subscript[p, 0])=(1/(1+(e^(-β (ϵ-μ))))), (Subscript[p, 1])=(1/((e^(β (ϵ-μ)))+1)).

Manipulate[Block[{β=1/T},
Plot[{1/(1+Exp[- β(ϵ-μ)]),1/(Exp[ β(ϵ-μ)]+1)},{ϵ,-2,2},Epilog->Block[{δ},δ= -2UnitStep[μ]+1;{Blue,Text["\!\(p\_0\)",{μ+δ,1/(1+Exp[-β δ])},{0,Sign[δ]}],Red,Text["\!\(p\_1\)",{μ+δ,1/(1+Exp[β δ])},{0,-Sign[δ]}]}],GridLines->{{μ},{0,1}},PlotRange->{0,1},PlotRangePadding->Scaled[.05],Axes->False,PlotLabel->StringTemplate["\!\(k\_B T=`1`,μ=`2`\)"][Round[T,0.01],Round[μ,0.01]],FrameLabel->{"\!\(ϵ\)","\!\(p\_n\)"},ImageSize->180]],{{T,0.2},0.01,1,ImageSize->Small},{{μ,0},-1,1,ImageSize->Small},Paneled->False]

Average occupation number

⟨n⟩=(∑_(n=0,1))n (Subscript[p, n])=(Subscript[p, 1])=(Subscript[n, F])(β(ϵ-μ)),

where (Subscript[n, F]) is the Fermi-Dirac distribution function

describing the average number of fermions occupying the mode of energy ϵ.

Average energy

⟨E⟩=(∑_(n=0,1))E (Subscript[p, n])
=ϵ (Subscript[n, F])(β(ϵ-μ)).

Average fermion (electron) number

⟨N⟩=(∑_(n=0,1))N (Subscript[p, n])
=(Subscript[n, F])(β(ϵ-μ)).

Note: Here ⟨N⟩=⟨n⟩ because the system has only one single mode.

Fermi Energy

(Subscript[m, eff]) - effective mass of electron in the crystal.

k - quasi-momentum. Each k labels a mode (Ket[k]) that can be occupied by the electron. More precisely, the electron also carries spin, that can be in one of the two degenerated spin states, so the electronic mode is jointly labeled by (Ket[k,σ]) (σ=↑,↓ ).

N=2(∑_k)(Subscript[n, F])(β((Subscript[ϵ, k])-μ))
=2(V/(((2π))^D))∫(d^D)k (Subscript[n, F])(β((Subscript[ϵ, k])-μ)).

Prefactor of 2 accounts for the two degenerated spin states (the ways that electrons occupy both spins states are identical).

D - dimension of space.

V - volume of the system. Such that the discretization unit in the momentum space is (((2π))^D)/V.

Fermi temperature (Subscript[T, F])

(Subscript[k, B])(Subscript[T, F])=(Subscript[ϵ, F]).

Fermi momentum (Subscript[k, F])

(((ℏ^2)(Subscript[k, F]^2))/(2(Subscript[m, eff])))=(Subscript[ϵ, F]).

Fermi velocity (Subscript[v, F])

(Subscript[v, F])=((ℏ (Subscript[k, F]))/(Subscript[m, eff])),

hence (Subscript[ϵ, F])=(1/2)(Subscript[m, eff])(Subscript[v, F]^2).

Number of electrons at T=0. At zero temperature, the Fermi-Dirac distribution function reduces to a step function (written as Θ(x)=1 for x≥0 and Θ(x)=0 for x<0), such that Eq. () becomes

N=2(V/(((2π))^D))∫(d^D)k Θ((Subscript[ϵ, F])-(Subscript[ϵ, k])),
=2(V/(((2π))^D))(∫^(k<(Subscript[k, F])))(d^D)k 
=2(V/(((2π))^D))(∫^(k<(Subscript[k, F])))(Subscript[A, D])(k^(D-1))dk
=2((V (Subscript[A, D])(Subscript[k, F]^D))/(D (((2π))^D))),

(Subscript[A, D])=2(π^(D/2))/Γ((D/2)) - area of a (D-1)-dimensional unit sphere (in the D-dimensional space). (Subscript[A, D])(Subscript[k, F]^D)/D is the volume of a ball of radius (Subscript[k, F]).

At T=0, the electrons simply fill a ball in the momentum space of radius (Subscript[k, F]). The surface is the ball is called the Fermi surface.

From Eq. (), we can solve for the Fermi momentum (Subscript[k, F])

(Subscript[k, F])=2π (((((D )/(2(Subscript[A, D])))(N/V)))^(1/D))

or more explicitly (let n=N/V)

(Subscript[k, F])=({{{, ({{π n/2, D=1}, {(((2π n))^(1/2)), D=2}, {(((3(π^2)n))^(1/3)), D=3}, {…, …}})}}).

Table[𝔇->PowerExpand[2π(𝔇/(2*2π^(𝔇/2)/Gamma[𝔇/2])n)^(1/𝔇)],{𝔇,3}]

Enumerate (Subscript[k, F]).

Correspondingly the Fermi energy

(Subscript[ϵ, F])=(((((2π))^2)(ℏ^2))/(2(Subscript[m, eff])))(((((D )/(2(Subscript[A, D])))(N/V)))^(2/D)).

or more explicitly (let n=N/V)

(Subscript[ϵ, F])=((ℏ^2)/(2(Subscript[m, eff])))({{{, ({{(((π n/2))^2), D=1}, {2π n, D=2}, {(((3(π^2)n))^(2/3)), D=3}, {…, …}})}}).

For typical metals, the Fermi energy is of the order (Subscript[ϵ, F])∼10eV, which corresponds to a Fermi temperature of the order (Subscript[T, F])∼(10^5)K, much higher than the room temperature (also much higher than the melting temperature of metals).

Therefore, most electrons in the metal are deep in the Fermi sea and do not participate in thermodynamic or transport process. Such many-body state of electrons is also called degenerate Fermi gas, which also happens in other fermion systems as long as T≪(Subscript[T, F]) (cold atomic gas (Subscript[T, F])∼(10^(-6))K, white dwarf star (Subscript[T, F])∼(10^9)K).

For typical metals, the Fermi velocity can be as large as (Subscript[v, F])~0.01c (c - the speed of light). The high velocity of electron originated from the Pauli exclusion principle — all low momentum (low velocity) states are filled, and the active electrons around the Fermi surface has to move with high velocity.

Density of States

E=2(∑_k)(Subscript[ϵ, k]) (Subscript[n, F])(β((Subscript[ϵ, k])-μ)),
N=2(∑_k)(Subscript[n, F])(β((Subscript[ϵ, k])-μ)).

(∑_k)=(V/(((2π))^D))∫(d^D)k =V∫d(Subscript[ϵ, k])g((Subscript[ϵ, k])),

Derive Eq. ().

(Subscript[ϵ, k])=(((ℏ^2)(k^2))/(2(Subscript[m, eff])))=>k=((Sqrt[2(Subscript[m, eff])ϵ])/ℏ),

(1/(((2π))^D))∫(d^D)k =(1/(((2π))^D))∫(Subscript[A, D])(k^(D-1))dk
=(((((2(Subscript[m, eff])))^(D/2))(Subscript[A, D]))/(2(((2π ℏ))^D)))∫(ϵ^(D/2-1))d ϵ.

g(ϵ)=(((((2(Subscript[m, eff])))^(D/2))(Subscript[A, D]))/(2(((2π ℏ))^D)))(ϵ^(D/2-1))
=((D n)/(4(Subscript[ϵ, F])))((((ϵ/(Subscript[ϵ, F]))))^(D/2-1)),

Column[Table[Plot[𝔇/4 x^(𝔇/2-1)UnitStep[x],{x,-1,3},Epilog->Text[StringTemplate["\!\(D=`1`\)"][𝔇],Scaled@{1,1},{1.2,1}],PlotStyle->Black,Filling->0,PlotRange->{0,1.7},PlotRangePadding->Scaled[.03],Exclusions->None,FrameLabel->{"\!\( ϵ/ϵ\_F\)","\!\(g(ϵ)ϵ\_F/n\)"},ImagePadding->{{40,2},{If[𝔇<3,2,42],2}},ImageSize->160],{𝔇,3}],Spacings->0]

E=2V∫dϵ g(ϵ)ϵ (Subscript[n, F])(β(ϵ-μ)),
N=2V∫dϵ g(ϵ)(Subscript[n, F])(β(ϵ-μ)).

Low-Temperature Expansion

I=(∫_(-∞)^∞)dϵ f(ϵ)(Subscript[n, F])(β(ϵ-μ)).

I=(∫_(-∞)^μ)dϵ f(ϵ)+(∑_(n=1)^∞)((2((((Subscript[k, B])T))^(2n)))/((2n-1)!))(∂_μ^(2n-1))f(μ) (Subscript[I, 2n-1]),

(Subscript[I, 2n-1])=(∫_0^∞)((x^(2n-1))/((e^x)+1))dx=(1-(2/(4^n)))Γ(2n)ζ(2n).

Derive Eq. () and Eq. ().

δ (Subscript[n, F])(β(ϵ-μ))=(Subscript[n, F])(β(ϵ-μ))-Θ(μ-ϵ).

Plot[1/(Exp[x]+1)-UnitStep[-x],{x,-5,5},PlotStyle->Black,FrameLabel->{"\!\( x\)","\!\(δ n\_F(x)\)"},ImageSize->180]

I=(∫_(-∞)^∞)dϵ f(ϵ)Θ(μ-ϵ)+(∫_(-∞)^∞)dϵ f(ϵ)δ (Subscript[n, F])(β(ϵ-μ))
=(∫_(-∞)^μ)dϵ f(ϵ)+(∫_(-∞)^∞)dϵ ((∑_(n=0)^∞)(((∂_μ^n)f(μ)(((ϵ-μ))^n))/(n!)))δ (Subscript[n, F])(β(ϵ-μ))
=(∫_(-∞)^μ)dϵ f(ϵ)+(∑_(n=0)^∞)(((∂_μ^n)f(μ))/(n!))((((Subscript[k, B])T))^(n+1))(∫_(-∞)^∞)dx (x^n)δ (Subscript[n, F])(x)

(∑_(n=1)^∞)(((∂_μ^(2n-1))f(μ))/((2n-1)!))((((Subscript[k, B])T))^(2n))(∫_(-∞)^∞)dx (x^(2n-1))δ (Subscript[n, F])(x),

(∫_(-∞)^∞)dϵ (x^(2n-1))δ (Subscript[n, F])(x)
=2(∫_0^∞)dϵ (x^(2n-1))δ (Subscript[n, F])(x)
=2(∫_0^∞)dϵ (x^(2n-1))(Subscript[n, F])(x)
=2(∫_0^∞)dϵ ((x^(2n-1))/((e^x)+1))
=2(Subscript[I, 2n-1]),

(∫_(-∞)^∞)dϵ f(ϵ)(Subscript[n, F])(β(ϵ-μ))=(∫_(-∞)^μ)dϵ f(ϵ)+((((π (Subscript[k, B])T))^2)/6)f'(μ)+…

E(μ,T)=((D )/(D+2))N (Subscript[ϵ, F])+2(μ-(Subscript[ϵ, F]))g((Subscript[ϵ, F]))(Subscript[ϵ, F])V+((((π (Subscript[k, B])T))^2)/3)(g(μ)+μ g'(μ))V+…,
N(μ,T)=N+2(μ-(Subscript[ϵ, F]))g((Subscript[ϵ, F]))V+((((π (Subscript[k, B])T))^2)/3)g'(μ)V+….

Derive Eq. ().

E=2V(∫_(-∞)^∞)dϵ g(ϵ)ϵ (Subscript[n, F])(β(ϵ-μ))
=2V((∫_(-∞)^μ)dϵ g(ϵ)ϵ+((((π (Subscript[k, B])T))^2)/6)(g(μ)+μ g'(μ))+…)
=2V((∫_(-∞)^(Subscript[ϵ, F]))dϵ g(ϵ)ϵ+(μ-(Subscript[ϵ, F]))g((Subscript[ϵ, F]))(Subscript[ϵ, F])+((((π (Subscript[k, B])T))^2)/6)(g(μ)+μ g'(μ))+…)
=((D )/(D+2))N (Subscript[ϵ, F])+2V(μ-(Subscript[ϵ, F]))g((Subscript[ϵ, F]))(Subscript[ϵ, F])+((((π (Subscript[k, B])T))^2)/3)V(g(μ)+μ g'(μ))+…,

(∫_(-∞)^(Subscript[ϵ, F]))dϵ g(ϵ)ϵ=((D n)/(4(Subscript[ϵ, F])))(∫_0^(Subscript[ϵ, F]))dϵ ((((ϵ/(Subscript[ϵ, F]))))^(D/2-1))ϵ
=((D n)/(4(Subscript[ϵ, F])))(Subscript[ϵ, F]^2)(∫_0^1)dx (x^(D/2))
=((2D n)/(4(D+2)))(Subscript[ϵ, F]).

N=2V(∫_(-∞)^∞)dϵ g(ϵ)(Subscript[n, F])(β(ϵ-μ))
=2V((∫_(-∞)^μ)dϵ g(ϵ)+((((π (Subscript[k, B])T))^2)/6)g'(μ)+…)
=2V((∫_(-∞)^(Subscript[ϵ, F]))dϵ g(ϵ)+(μ-(Subscript[ϵ, F]))g((Subscript[ϵ, F]))+((((π (Subscript[k, B])T))^2)/6)g'(μ)+…)
=N+2(μ-(Subscript[ϵ, F]))g((Subscript[ϵ, F]))V+((((π (Subscript[k, B])T))^2)/3)g'(μ)V+…,

(∫_(-∞)^(Subscript[ϵ, F]))dϵ g(ϵ)=((D n)/(4(Subscript[ϵ, F])))(∫_0^(Subscript[ϵ, F]))dϵ ((((ϵ/(Subscript[ϵ, F]))))^(D/2-1))
=((D n)/(4(Subscript[ϵ, F])))(Subscript[ϵ, F])(∫_0^1)dx (x^(D/2-1))
=(n/2).

Heat Capacity

Heat capacity: the rate of energy E change with respect to the temperature T

(Subscript[C, V])=((∂E)/(∂T))=(((π^2)(Subscript[k, B]^2)T)/3)(2g((Subscript[ϵ, F]))V).

For quadratic dispersion, the density of state g(ϵ) given in Eq. (),

(Subscript[C, V])=((π^2)/3)(((D N (Subscript[k, B]))/2))((T/(Subscript[T, F]))).

Specific heat: heat capacity per electron

For quadratic dispersion

(Subscript[c, V])=((π^2)/3)(((D (Subscript[k, B]))/2))((T/(Subscript[T, F]))).

D (Subscript[k, B])/2 is the classical result for specific heat of a gas, where every atom in the gas can exchange energy with the heat bath. However, in metals, on those electrons near the Fermi surface can absorb/release energy, so the specific heat of such degenerate Fermi gas is smaller than the classical gas by a factor of T/(Subscript[T, F]), which is a tiny ratio T/(Subscript[T, F])∼0.01 in typical metals.

At low-temperature, the specific heat of electrons in metals grows with temperature linearly (Subscript[c, V])∼T (regardless of the spatial dimension D of the system).

(Subscript[c, V])=γ T+α (T^D),

Charge Compressibility

Charge compressibility (quantum capacitance): the rate of charge  Q changes with respect to the electric potential U

(Subscript[C, Q])=((∂Q)/(∂U))=(e^2)((∂N)/(∂μ)).

e - electric charge of electron.

In the presence of electric potential U, the electron chemical potential shifts by

μ→μ-e U.

As a result the number of electron changes by

N→N+((∂N)/(∂μ))(-e U),

and the total charge Q=-e N changes by

Q→Q+(e^2)((∂N)/(∂μ))U.

Using the result in Eq. (),

N(μ,T)=N+2(μ-(Subscript[ϵ, F]))g((Subscript[ϵ, F]))V+…,

we have

((∂N)/(∂μ))=2g((Subscript[ϵ, F]))V.

So the charge compressibility is given by

(Subscript[C, Q])=(e^2)((∂N)/(∂μ))=(e^2)(2 g((Subscript[ϵ, F]))V).

For quadratic dispersion, the density of state g(ϵ) given in Eq. (),

(Subscript[C, Q])=(((e^2)D N)/(2(Subscript[ϵ, F])))

Specific charge compressibility

For quadratic dispersion,

(Subscript[c, Q])=(((e^2)D)/(2(Subscript[ϵ, F]))).

Spin Susceptibility

Spin susceptibility: the rate of magnetization M changes with respect to the external magnetic field B

χ=((∂M)/(∂B))

(Subscript[μ, B])=e ℏ/(2(Subscript[m, e])) - Bohr magneton, the magnetic moment carried by the electron.

In the presence of a magnetic field B, the up-spin and down-spin electron energy split

(Subscript[ϵ, k,↑ ])=(Subscript[ϵ, k])+(Subscript[μ, B])B,
(Subscript[ϵ, k,↓ ])=(Subscript[ϵ, k])-(Subscript[μ, B])B.

Magnetization (average magnetic moment per electron) responds to the magnetic field as

M=(Subscript[μ, B])(((Subscript[N, ↑])-(Subscript[N, ↓]))/N)
=((Subscript[μ, B])/N)(((∂(Subscript[N, ↑]))/(∂(Subscript[μ, ↑])))((Subscript[μ, B])B)-((∂(Subscript[N, ↓]))/(∂(Subscript[μ, ↓])))(-(Subscript[μ, B])B))
=((Subscript[μ, B]^2)/N)(((∂(Subscript[N, ↑]))/(∂(Subscript[μ, ↑])))+((∂(Subscript[N, ↓]))/(∂(Subscript[μ, ↓]))))B
=((Subscript[μ, B]^2)/N)((∂N)/(∂μ))B.

So the spin susceptibility can be calculated as

For quadratic dispersion, the density of state g(ϵ) given in Eq. (),

χ=(((Subscript[μ, B]^2)D)/(2(Subscript[ϵ, F]))).

Wilson Ratio

(Subscript[c, V])=(((π^2)(Subscript[k, B]^2)T)/3)((2g((Subscript[ϵ, F])))/n),
(Subscript[c, Q])=(e^2)((2g((Subscript[ϵ, F])))/n),
χ=(Subscript[μ, B]^2)((2g((Subscript[ϵ, F])))/n).

For non-interacting electrons, the Wilson ratio should be 1, which is expected to hold for typical metals.

Deviation of the Wilson ratio from unity signifies the interaction effect between electrons is important in the material.

Transport Properties

Boltzmann Transport

Equilibrium distribution. Without external perturbations, electrons relax to the equilibrium distribution

(Subscript[f, eq])(x,k)=(Subscript[n, F])(β((Subscript[ϵ, k])-μ))=(1/((e^(β((Subscript[ϵ, k])-μ)))+1)).

At the low-temperature limit (T→0),

(Subscript[f, eq])(x,k)=Θ((Subscript[ϵ, F])-(Subscript[ϵ, k])) 
=>((∂f)/(∂(Subscript[ϵ, k])))=-δ((Subscript[ϵ, F])-(Subscript[ϵ, k])).

Relaxation of near equilibrium distribution. If f deviated from (Subscript[f, eq]), it will relax to (Subscript[f, eq]) in a characteristic time scale τ, called the relaxation time,

f(t)=(Subscript[f, eq])+(f(0)-(Subscript[f, eq]))(e^(-t/τ)).

The relaxation dynamics can be described by the following differential equation

((df)/(dt))=-((f-(Subscript[f, eq]))/τ).

On the other hand, by the chain rule of differentiation

((df)/(dt))=((∂f)/(∂t))+((∂f)/(∂x))·((dx)/(dt))+((∂f)/(∂k))·((dk)/(dt)).

By semiclassical equation of motion

((dx)/(dt))=(Subscript[v, k])=(1/ℏ)(∂_k)(Subscript[ϵ, k]),
((dk)/(dt))=(1/ℏ)F,

(Subscript[v, k]) - group velocity of the electron at momentum k,

F - external force acting on the electron.

Eq. () can be written as

((df)/(dt))=(∂_t)f+(Subscript[v, k])·(∂_x)f+(1/ℏ)F·(∂_k)f.

(Subscript[v, k])·(∂_x)f+(1/ℏ)F·(∂_k)f=-((f-(Subscript[f, eq]))/τ).

Electrical Conductivity

F=-e E.

-(e/ℏ)E·(∂_k)f=-((f-(Subscript[f, eq]))/τ),

f=(Subscript[f, eq])+((e τ)/ℏ)E·(∂_k)f.

To the first order in E (linear response),

f≃(Subscript[f, eq])+((e τ)/ℏ)E·(∂_k)(Subscript[f, eq])
=(Subscript[f, eq])+e τ E·(Subscript[v, k])(∂_ϵ)(Subscript[f, eq]).

f is deviated from (Subscript[f, eq]). The main effect is to shift the Fermi sea in the momentum space by the amount of -(e τ/ℏ)E.

Block[{β=5,τ=0.2,μ=1,ϵ,f,df,g1,g2,g3},
ϵ[{kx_,ky_}]:=(kx^2+ky^2)/2;
f[k_]:=1/(Exp[β(ϵ[k]-μ)]+1);
df[k:{kx_,ky_}]:=-β τ  kx Exp[β(ϵ[k]-μ)]/(Exp[β(ϵ[k]-μ)]+1)^2;
{g1,g2,g3}=ContourPlot[Rescale[#1[{kx,ky}],{-1,1}],{kx,-2,2},{ky,-2,2},Epilog->{Translate[Block[{l=1},{AbsoluteThickness[1],Arrowheads[0.08],Arrow[{{0,0},#}]&/@{{l,0},{0,l}},Text["\!\(k\_x\)",{l,0},{0,-0.8}],Text["\!\(k\_y\)",{0,l},{-1.4,0.2}]}],{-1.8,-1.8}],Text[#2,Scaled@{0.98,1},{1,1}]},GridLines->{{0},{0}},PlotRangePadding->None,Contours->Range[0,1,0.1],ContourStyle->None,ColorFunctionScaling->False,FrameTicks->None,ImageSize->85]&@@@{{f,"\!\(f\_eq\)"},{df,"\!\(δ f\)"},{f[#]+df[#]&,"\!\( f\)"}};
Row[{g1,Style[" + ",18,FontFamily->"LaTeX"],g2,Style[" = ",18,FontFamily->"LaTeX"],g3}]
]

J=-2e∫(((d^D)k)/(((2π))^D))(Subscript[v, k])f
=-2e∫(((d^D)k)/(((2π))^D))(Subscript[v, k])(Subscript[f, eq])-2(e^2)τ∫(((d^D)k)/(((2π))^D))(Subscript[v, k]) E·(Subscript[v, k])(∂_ϵ)(Subscript[f, eq]).

The prefactor 2 comes from the spin degeneracy

The first term will vanish because there is no current when the system is at equilibrium.

Switch the integral from the momentum space to the energy domain

J=-2(e^2)τ∫d(Subscript[ϵ, k])g((Subscript[ϵ, k]))(Subscript[v, k]) E·(Subscript[v, k])(∂_ϵ)(Subscript[f, eq])
=-((2(e^2)τ (Subscript[v, F]^2))/3)E∫dϵ g(ϵ)(∂_ϵ)(Subscript[f, eq]).

The integral admits the following low-temperature expansion

∫dϵ g(ϵ)(∂_ϵ)(Subscript[f, eq])=-g(μ)-((π^2)/6)(((∂_μ^2)g(μ))/(β^2))+…,

Derive Eq. ().

(Subscript[f, eq])=(1/((e^(β (ϵ-μ)))+1))=>(∂_ϵ)(Subscript[f, eq])=-(β/4)(sech^2)((β(ϵ-μ))/2),

(∫_(-∞)^∞)dϵ g(ϵ)(∂_ϵ)(Subscript[f, eq])=-(β/4)(∫_(-∞)^∞)dϵ g(ϵ)(sech^2)((β(ϵ-μ))/2)
=-(β/4)(∫_(-∞)^∞)dϵ ((∑_(n=0)^∞)(((∂_μ^n)g(μ)(((ϵ-μ))^n))/(n!)))(sech^2)((β(ϵ-μ))/2)
=-(β/4)(∑_(n=0)^∞)(((∂_μ^n)g(μ))/(n!))(∫_(-∞)^∞)dϵ (((ϵ-μ))^n)(sech^2)((β(ϵ-μ))/2).

(∫_(-∞)^∞)dϵ g(ϵ)(∂_ϵ)(Subscript[f, eq])=-(β/4)(∑_(n=0)^∞)(((∂_μ^(2n))g(μ))/((2n)!))(∫_(-∞)^∞)dϵ (((ϵ-μ))^(2n))(sech^2)((β(ϵ-μ))/2)
=-(1/4)(∑_(n=0)^∞)(((∂_μ^(2n))g(μ))/((2n)!(β^(2n))))(∫_(-∞)^∞)dx (x^(2n))(sech^2)(x/2)

(∫_(-∞)^∞)dϵ g(ϵ)(∂_ϵ)(Subscript[f, eq])=-g(μ)-((π^2)/6)(((∂_μ^2)g(μ))/(β^2))-((7(π^4))/360)(((∂_μ^4)g(μ))/(β^4))+…

Table[2n->Integrate[x^(2n)Sech[x/2]^2,{x,-∞,∞}],{n,0,2}]

Integration.

Take the low-temperature limit and use μ=(Subscript[ϵ, F]) as T→0,

J=((2(e^2)τ (Subscript[v, F]^2))/3)g((Subscript[ϵ, F])) E .

J=σ E.

Thermal Conductivity

(∂_x)f=(∂_T)f (∂_x)T.

(Subscript[v, k])·(∂_x)T  (∂_T)f=-((f-(Subscript[f, eq]))/τ),

f=(Subscript[f, eq])-τ (Subscript[v, k])·(∂_x)T  (∂_T)f.

To the first order in (∂_x)T (linear response),

f≃(Subscript[f, eq])-τ (Subscript[v, k])·(∂_x)T  (∂_T)(Subscript[f, eq]).

Given (Subscript[f, eq])=((((e^((ϵ-μ)/T))+1))^(-1)), (∂_T)(Subscript[f, eq]) and (∂_ϵ)(Subscript[f, eq]) are related by

(∂_T)(Subscript[f, eq])=-((ϵ-μ)/T)(∂_ϵ)(Subscript[f, eq]),

therefore

f≃(Subscript[f, eq])+(((Subscript[ϵ, k])-μ)/T)τ (Subscript[v, k])·(∂_x)T  (∂_ϵ)(Subscript[f, eq]).

(Subscript[J, Q])=2∫(((d^D)k)/(((2π))^D))((Subscript[ϵ, k])-μ)(Subscript[v, k])f
=2∫(((d^D)k)/(((2π))^D))((Subscript[ϵ, k])-μ)(Subscript[v, k])(Subscript[f, eq])+2∫(((d^D)k)/(((2π))^D))((((((Subscript[ϵ, k])-μ))^2)τ)/T)(Subscript[v, k])  (Subscript[v, k])·(∂_x)T  (∂_ϵ)(Subscript[f, eq]).

The prefactor 2 comes from the spin degeneracy

The first term will vanish because there is no current when the system is at equilibrium.

Switch the integral from the momentum space to the energy domain

(Subscript[J, Q])=((2τ)/T)∫(((d^D)k)/(((2π))^D))((((Subscript[ϵ, k])-μ))^2)(Subscript[v, k])  (Subscript[v, k])·(∂_x)T  (∂_ϵ)(Subscript[f, eq])
=((2τ)/T)∫d(Subscript[ϵ, k])g((Subscript[ϵ, k]))((((Subscript[ϵ, k])-μ))^2)(Subscript[v, k])  (Subscript[v, k])·(∂_x)T  (∂_ϵ)(Subscript[f, eq])
=((2τ (Subscript[v, F]^2))/3)(((∂_x)T)/T)∫dϵ g(ϵ)(((ϵ-μ))^2) (∂_ϵ)(Subscript[f, eq]).

Use the low-temperature expansion in Eq. ()

∫dϵ g(ϵ)(((ϵ-μ))^2) (∂_ϵ)(Subscript[f, eq])=-((π^2)/3)((g(μ))/(β^2))+…,

therefore

(Subscript[J, Q])=-((2τ (Subscript[v, F]^2))/3)g((Subscript[ϵ, F]))(((π^2)(Subscript[k, B]^2)T)/3)(∂_x)T.

(Subscript[J, Q])=-κ (∂_x)T.

Lorenz Number

({{Material, L×(10^8)}, {Sodium (Na), 2.12}, {Magnesium (Mg), 2.14}, {Copper (Cu), 2.20}, {Lithium (Li), 2.22}, {Zinc (Zn), 2.33}, {Silver (Ag), 2.37}, {Gold (Au), 2.40}, {Cadmium (Cd), 2.43}, {Lead (Pb), 2.56}, {Platinum (Pt), 2.60}, {Iron (Fe), 2.61}, {Tungsten (W), 3.20}, {Bismuth (Bi), 3.53}})

The fact that many metals have roughly the same κ/(σ T) ratio is also known as the Wiedemann-Franz law.

Significance deviation of the Lorenz number of this value is an indication of the importance of electron interaction in the material.

Interacting Electrons

Interacting Two-Electron Systems

Warmup: Classical Interacting Particles

Graphics[{{AbsoluteThickness[1],Arrow[{{-1,0},{1,0}}],Text["\!\( x\)",{1,0},{0,1}]},Translate[{#3,Text["\!\(x\_"<>#1<>"\)",{0,0},{0,1}],Disk[{0,0},0.04]},#2]&@@@{{"1",{-0.5,0},Green},{"2",{0.2,0},Purple}}},ImageSize->180]

Suppose there is some uncertainty about the positions → should think about the probability distribution of the particles p((Subscript[x, 1])) and p((Subscript[x, 2]))

Block[{σ=0.2},Graphics[{{AbsoluteThickness[1],Arrow[{{-1,0},{1,0}}],Text["\!\( x\)",{1,0},{0,1}]},Translate[{#3,Text["\!\(x\_"<>#1<>"\)",{0,0},{0,1}],Table[{Opacity[1-0.9Abs[x]],Disk[{2σ x,0},0.04]},{x,SortBy[Range[-1,1,0.2],Abs]}],Line[Table[{σ x,0.2Exp[-x^2/2]},{x,-2.5,2.5,0.01}]]},#2]&@@@{{"1",{-0.5,0},Green},{"2",{0.2,0},Purple}}},ImageSize->180]]

Or more precisely, the joint probability distribution

p((Subscript[x, 1]),(Subscript[x, 2]))=p((Subscript[x, 1]))p((Subscript[x, 2])).

Block[{σ=0.4,xc={-0.5,0.2},p},
p[x_]:=Exp[-Total[(x-xc)^2]/(2σ^2)];
ContourPlot[Rescale[p[{x1,x2}],{-1,1}],{x1,-1,1},{x2,-1,1},ContourStyle->None,Epilog->{Dashed,Line[{{-1,-1},{1,1}}]},ColorFunction->"BlueRedSplit",ColorFunctionScaling->False,PlotRange->All,FrameLabel->{"\!\(x\_1\)","\!\(x\_2\)"},ImageSize->180]]

Suppose the interaction between the two particles is described by the interaction potential U((Subscript[x, 1])-(Subscript[x, 2])) (assuming translation symmetry and inversion symmetry), the average interaction energy is

(Subscript[E, U])=∫d(Subscript[x, 1])d(Subscript[x, 2])U((Subscript[x, 1])-(Subscript[x, 2]))p((Subscript[x, 1]),(Subscript[x, 2])).

(illustration of the integrand)

Block[{σ=0.4,xc={-0.5,0.2},p,U,gs},
p[x_]:=Exp[-Total[(x-xc)^2]/(2σ^2)];
U[x_]:=Block[{r},r=Abs[x[[1]]-x[[2]]];2Sech[3r/σ]-Sech[r/σ]];
gs=ContourPlot[Rescale[#1[{x1,x2}],#3],{x1,-1,1},{x2,-1,1},ContourStyle->None,Contours->Range[0,1,0.05],Epilog->{Dashed,Line[{{-1,-1},{1,1}}]},ColorFunction->#2,ColorFunctionScaling->False,PlotRange->All,PlotLabel->#4,FrameLabel->{"\!\(x\_1\)","\!\(x\_2\)"},FrameTicks->None,BaselinePosition->Scaled@0.5,ImageSize->100]&@@@{{U,"BlueRedSplit",{-1,1},"\!\(U(x\_i)\)"},{p,"BlueRedSplit",{-1,1},"\!\(p(x\_i)\)"},{p[#]U[#]&,"BlueRedSplit",{-0.8,0.8},"\!\(U(x\_i)p(x\_i)\)"}};
Row[{gs[[1]],Style["×",16,FontFamily->"LaTeX"],gs[[2]],Style["=",16,FontFamily->"LaTeX"],gs[[3]]}]]

Interacting Fermions

First, fermions are identical particles. There is no way to distinguish which electron is the first/second one. This implies that the probability distribution must be invariant under the exchange of the two fermions, i.e.

p((Subscript[x, 1]),(Subscript[x, 2]))=p((Subscript[x, 2]),(Subscript[x, 1])).

Secondly, fermions are quantum particles. The probability p((Subscript[x, 1]),(Subscript[x, 2])) should be modeled by the underlying probability amplitude (the wave function) ψ((Subscript[x, 1]),(Subscript[x, 2])) as

p((Subscript[x, 1]),(Subscript[x, 2]))=((|ψ((Subscript[x, 1]),(Subscript[x, 2]))|)^2).

((|ψ((Subscript[x, 1]),(Subscript[x, 2]))|)^2)=((|ψ((Subscript[x, 2]),(Subscript[x, 1]))|)^2).

ψ((Subscript[x, 1]),(Subscript[x, 2]))=({{{, ({{+ψ((Subscript[x, 2]),(Subscript[x, 1])), bosons,}, {-ψ((Subscript[x, 2]),(Subscript[x, 1])), fermions.}})}})

One can already see, if (Subscript[ψ, A])=(Subscript[ψ, B]) are the same state, the two-body wave function vanishes automatically => Pauli exclusion principle: fermions can not double occupy the same state.

Consider two localized orbitals

(Subscript[ψ, A])(x)∝exp(-((((x-(Subscript[x, A])))^2)/2)),
(Subscript[ψ, B])(x)∝exp(-((((x-(Subscript[x, B])))^2)/2)).

Block[{xA=-1,xB=1},
Plot[{Exp[-(x-xA)^2/2]/Sqrt[2π],Exp[-(x-xB)^2/2]/Sqrt[2π]},{x,-5,5},PlotStyle->{Green,Purple},Prolog->{{AbsoluteThickness[1],Dashed,Arrowheads[0.06],Green,Line[{xA,#}&/@{0,0.4}],Text[Style["\!\(x\_A\)",Background->White],{xA,0.19}],Arrow[{{xA,0.12},{xA,0.1}}],Purple,Line[{xB,#}&/@{0,0.4}],Text[Style["\!\(x\_B\)",Background->White],{xB,0.19}],Arrow[{{xB,0.12},{xB,0.1}}]}},Epilog->{{Green,Text["\!\(ψ\_A(x)\)",{-2.,0.3},{1,0}]},{Purple,Text["\!\(ψ\_B(x)\)",{2.,0.3},{-1,0}]}},FrameLabel->{"\!\( x\)","\!\(ψ(x)\)"},ImageSize->180]]

The two-body wave function can be constructed by antisymmetrization.

Block[{xA=-1,xB=1,ψA,ψB,gs},
ψA[x_]:=Exp[-(x-xA)^2/2];
ψB[x_]:=Exp[-(x-xB)^2/2];
gs=ContourPlot[Rescale[#1[x1,x2],#3],{x1,-3,3},{x2,-3,3},ContourStyle->None,Contours->Range[0,1,0.05],Epilog->{Dashed,Line[{{-5,-5},{5,5}}]},ColorFunction->#2,ColorFunctionScaling->False,PlotRange->All,PlotLabel->#4,FrameLabel->{"\!\(x\_1\)","\!\(x\_2\)"},FrameTicks->None,BaselinePosition->Scaled@0.5,ImageSize->100]&@@@{{ψA[#1]ψB[#2]&,"BlueRedSplit",{-1,1},"\!\(ψ\_A ψ\_B\)"},{ψB[#1]ψA[#2]&,"BlueRedSplit",{-1,1},"\!\(ψ\_B ψ\_A\)"},{ψA[#1]ψB[#2]-ψB[#1]ψA[#2]&,"BlueRedSplit",{-1,1},"\!\(ψ(x\_1,x\_2)\)"}};
Row[{gs[[1]],Style["-",16,FontFamily->"LaTeX"],gs[[2]],Style["=",16,FontFamily->"LaTeX"],gs[[3]]}]]

Probability distribution

p((Subscript[x, 1]),(Subscript[x, 2]))=((|ψ((Subscript[x, 1]),(Subscript[x, 2]))|)^2)
=(1/2)((|(Subscript[ψ, A])((Subscript[x, 1]))(Subscript[ψ, B])((Subscript[x, 2]))-(Subscript[ψ, B])((Subscript[x, 1]))(Subscript[ψ, A])((Subscript[x, 2]))|)^2)
=(1/2)(((|(Subscript[ψ, A])((Subscript[x, 1]))|)^2)((|(Subscript[ψ, B])((Subscript[x, 2]))|)^2)+((|(Subscript[ψ, B])((Subscript[x, 1]))|)^2)((|(Subscript[ψ, A])((Subscript[x, 2]))|)^2))-(1/2)((Subscript[ψ, A]^( * ))((Subscript[x, 1]))(Subscript[ψ, B]^( * ))((Subscript[x, 2]))(Subscript[ψ, B])((Subscript[x, 1]))(Subscript[ψ, A])((Subscript[x, 2]))+(Subscript[ψ, B]^( * ))((Subscript[x, 1]))(Subscript[ψ, A]^( * ))((Subscript[x, 2]))(Subscript[ψ, A])((Subscript[x, 1]))(Subscript[ψ, B])((Subscript[x, 2]))),

which contains two contributions

p((Subscript[x, 1]),(Subscript[x, 2]))=(Subscript[p, cl])((Subscript[x, 1]),(Subscript[x, 2]))-(Subscript[p, qu])((Subscript[x, 1]),(Subscript[x, 2])),

the classical part is just the symmetrized product distribution

(Subscript[p, cl])((Subscript[x, 1]),(Subscript[x, 2]))=(1/2)((Subscript[p, A])((Subscript[x, 1]))(Subscript[p, B])((Subscript[x, 2]))+(Subscript[p, B])((Subscript[x, 1]))(Subscript[p, A])((Subscript[x, 2]))),

where (Subscript[p, A/B])(x)=((|(Subscript[ψ, A/B])(x)|)^2).

the quantum part is a correction that originates from quantum interference between fermions

(Subscript[p, qu])((Subscript[x, 1]),(Subscript[x, 2]))=(1/2)((Subscript[ψ, A]^( * ))((Subscript[x, 1]))(Subscript[ψ, B]^( * ))((Subscript[x, 2]))(Subscript[ψ, B])((Subscript[x, 1]))(Subscript[ψ, A])((Subscript[x, 2]))+(Subscript[ψ, B]^( * ))((Subscript[x, 1]))(Subscript[ψ, A]^( * ))((Subscript[x, 2]))(Subscript[ψ, A])((Subscript[x, 1]))(Subscript[ψ, B])((Subscript[x, 2]))).

Block[{xA=-1,xB=1,ψA,ψB,gs},
ψA[x_]:=Exp[-(x-xA)^2/2];
ψB[x_]:=Exp[-(x-xB)^2/2];
gs=ContourPlot[Rescale[#1[x1,x2],#3],{x1,-3,3},{x2,-3,3},ContourStyle->None,Contours->Range[0,1,0.05],Epilog->{Dashed,Line[{{-5,-5},{5,5}}]},ColorFunction->#2,ColorFunctionScaling->False,PlotRange->All,PlotLabel->#4,FrameLabel->{"\!\(x\_1\)","\!\(x\_2\)"},FrameTicks->None,BaselinePosition->Scaled@0.5,ImageSize->120]&@@@{{(Abs[ψA[#1]]^2Abs[ψB[#2]]^2+Abs[ψB[#1]]^2Abs[ψA[#2]]^2)&,"BlueRedSplit",{-1,1},"\!\(p\_cl(x\_1,x\_2)\)"},{(ψA[#1]ψB[#2]ψB[#1]ψA[#2]+ψB[#1]ψA[#2]ψA[#1]ψB[#2])&,"BlueRedSplit",{-1,1},"\!\(p\_qu(x\_1,x\_2)\)"},{Abs[ψA[#1]ψB[#2]-ψB[#1]ψA[#2]]^2&,"BlueRedSplit",{-1,1},"\!\( p(x\_1,x\_2)\)"}};
Row[gs]]

The average interaction energy

(Subscript[E, U])=∫d(Subscript[x, 1])d(Subscript[x, 2])U((Subscript[x, 1])-(Subscript[x, 2]))p((Subscript[x, 1]),(Subscript[x, 2]))

will also contain two terms

the Hartree energy (direct energy) from the classical contribution

In terms of the Dirac notation

(Subscript[E, Hart])=(Bra[(Subscript[ψ, A])(Subscript[ψ, B])])U(Ket[(Subscript[ψ, A])(Subscript[ψ, B])]).

Block[{xA=-1,xB=1,ψA,ψB,p,U,gs},
ψA[x_]:=Exp[-(x-xA)^2/2];
ψB[x_]:=Exp[-(x-xB)^2/2];
p[x1_,x2_]:=(Abs[ψA[x1]]^2Abs[ψB[x2]]^2+Abs[ψB[x1]]^2Abs[ψA[x2]]^2);
U[x1_,x2_]:=Block[{r},r=Abs[x1-x2];Sech[r]];
gs=ContourPlot[Rescale[#1[x1,x2],#3],{x1,-3,3},{x2,-3,3},ContourStyle->None,Contours->Range[0,1,0.05],Epilog->{Dashed,Line[{{-5,-5},{5,5}}]},ColorFunction->#2,ColorFunctionScaling->False,PlotRange->All,PlotLabel->#4,FrameLabel->{"\!\(x\_1\)","\!\(x\_2\)"},FrameTicks->None,BaselinePosition->Scaled@0.5,ImageSize->100]&@@@{{U,"BlueRedSplit",{-1,1},"\!\( U\)"},{p,"BlueRedSplit",{-1,1},"\!\( p\_cl\)"},{p[##]U[##]&,"BlueRedSplit",{-0.4,0.4},"\!\(U p\_cl\)"}};
Row[{gs[[1]],Style["×",16,FontFamily->"LaTeX"],gs[[2]],Style["=",16,FontFamily->"LaTeX"],gs[[3]]}]]

the Fock energy (exchange energy) from the quantum correction

In terms of the Dirac notation

(Subscript[E, Fock])=(Bra[(Subscript[ψ, A])(Subscript[ψ, B])])U(Ket[(Subscript[ψ, B])(Subscript[ψ, A])]).

Block[{xA=-1,xB=1,ψA,ψB,p,U,gs},
ψA[x_]:=Exp[-(x-xA)^2/2];
ψB[x_]:=Exp[-(x-xB)^2/2];
p[x1_,x2_]:=(ψA[x1]ψB[x2]ψB[x1]ψA[x2]+ψB[x1]ψA[x2]ψA[x1]ψB[x2]);
U[x1_,x2_]:=Block[{r},r=Abs[x1-x2];Sech[r]];
gs=ContourPlot[Rescale[#1[x1,x2],#3],{x1,-3,3},{x2,-3,3},ContourStyle->None,Contours->Range[0,1,0.05],Epilog->{Dashed,Line[{{-5,-5},{5,5}}]},ColorFunction->#2,ColorFunctionScaling->False,PlotRange->All,PlotLabel->#4,FrameLabel->{"\!\(x\_1\)","\!\(x\_2\)"},FrameTicks->None,BaselinePosition->Scaled@0.5,ImageSize->100]&@@@{{U,"BlueRedSplit",{-1,1},"\!\( U\)"},{p,"BlueRedSplit",{-1,1},"\!\( p\_qu\)"},{p[##]U[##]&,"BlueRedSplit",{-0.4,0.4},"\!\(U p\_qu\)"}};
Row[{gs[[1]],Style["×",16,FontFamily->"LaTeX"],gs[[2]],Style["=",16,FontFamily->"LaTeX"],gs[[3]]}]]

Hartree-Fock Approximation

H=(∑_i)(Underscript[(Underscript[((Subscript[p, i]^2)/(2m))+V((Subscript[x, i])), ⏟]), (Subscript[H, 0])])+(∑_(i≠j))U((Subscript[x, i])-(Subscript[x, j])).

It assumes that the N-body ground state can be approximated by the state of N fermions separately occupying N single-particle orbitals.

By optimizing the orbital wave functions to minimize the energy, one can estimate the ground state energy as well as the approximate ground state wave function.

H=(Subscript[H, 0])((Subscript[p, 1]),(Subscript[x, 1]))+(Subscript[H, 0])((Subscript[p, 2]),(Subscript[x, 2]))+U((Subscript[x, 1])-(Subscript[x, 2]))

ψ((Subscript[x, 1]),(Subscript[x, 2]))=(1/(Sqrt[2]))((Subscript[ψ, A])((Subscript[x, 1]))(Subscript[ψ, B])((Subscript[x, 2]))-(Subscript[ψ, B])((Subscript[x, 1]))(Subscript[ψ, A])((Subscript[x, 2]))).

(Bra[ψ])H(Ket[ψ])=(Bra[(Subscript[ψ, A])])(Subscript[H, 0])(Ket[(Subscript[ψ, A])])+(Bra[(Subscript[ψ, B])])(Subscript[H, 0])(Ket[(Subscript[ψ, B])])+(Underscript[(Underscript[(Bra[(Subscript[ψ, A])(Subscript[ψ, B])])U(Ket[(Subscript[ψ, A])(Subscript[ψ, B])]), ⏟]), Hartree])(Underscript[-(Underscript[(Bra[(Subscript[ψ, A])(Subscript[ψ, B])])U(Ket[(Subscript[ψ, B])(Subscript[ψ, A])]), ⏟]), Fock]).

the A orbital energy is (similar for B orbital)

(Bra[(Subscript[ψ, A])])(Subscript[H, 0])(Ket[(Subscript[ψ, A])])=∫dx (Subscript[ψ, A]^( * ))(x)(-((ℏ^2)/(2m))(∂_x^2)+V(x))(Subscript[ψ, A])(x).

the Hartree energy

(Bra[(Subscript[ψ, A])(Subscript[ψ, B])])U(Ket[(Subscript[ψ, A])(Subscript[ψ, B])])=∫d(Subscript[x, 1])d(Subscript[x, 2])U((Subscript[x, 1])-(Subscript[x, 2]))((|(Subscript[ψ, A])((Subscript[x, 1]))|)^2)((|(Subscript[ψ, B])((Subscript[x, 2]))|)^2),

the Fock energy

(Bra[(Subscript[ψ, A])(Subscript[ψ, B])])U(Ket[(Subscript[ψ, B])(Subscript[ψ, A])])=∫d(Subscript[x, 1])d(Subscript[x, 2])U((Subscript[x, 1])-(Subscript[x, 2]))(Subscript[ψ, A]^( * ))((Subscript[x, 1]))(Subscript[ψ, B]^( * ))((Subscript[x, 2]))(Subscript[ψ, B])((Subscript[x, 1]))(Subscript[ψ, A])((Subscript[x, 2])).

The minimal value of (Bra[ψ])H(Ket[ψ]) provides an estimation of the ground state energy.

The optimal orbitals can be used to construct the many-body state ψ((Subscript[x, 1]),(Subscript[x, 2])), which is a variational approximation of the ground state.

Two Electrons in a Box

U(x)=g sech((x/ξ)),

Plot[Sech[x],{x,-5,5},FrameLabel->{"\!\(x/ξ\)","\!\(U(x)/g\)"},PlotStyle->Black,ImageSize->180]

Non-interacting limit (g=0): the electrons will occupy the lowest two orbitals in the square well. The single-electron density distribution can be calculated from

ρ(x)=∫((|ψ(x,x')|)^2)dx'=((|(Subscript[ψ, A])(x)|)^2)+((|(Subscript[ψ, B])(x)|)^2).

Block[{ξ=0.25,g=0.0001,M=4,H0,U,ψs0,energy,ortho,Egs,sol,ψgs},
H0=SparseArray[Band[{1,1}]->Range[M]^2,{M,M}];U=Chop@Table[4g NIntegrate[Sech[(x1-x2)/ξ]Sin[m1 π x1]Sin[m2 π x2]Sin[m3 π x1]Sin[m4 π x2],{x1,0,1},{x2,0,1},AccuracyGoal->3,PrecisionGoal->2],{m1,M},{m2,M},{m3,M},{m4,M}];
ψs0=Normal@SparseArray[{{1,1}->1,{2,2}->1},{2,M}];
energy[ψs_List]:=ψs[[1]].H0.ψs[[1]]+ψs[[2]].H0.ψs[[2]]+tTr[{U,ψs[[1]],ψs[[2]],ψs[[1]],ψs[[2]]},{{1,5},{2,6},{3,7},{4,8}}]-tTr[{U,ψs[[1]],ψs[[2]],ψs[[2]],ψs[[1]]},{{1,5},{2,6},{3,7},{4,8}}];
ortho[ψs_List]:=Total[Flatten[ψs.Transpose[ψs]-IdentityMatrix[2]]^2];
{Egs,sol}=FindMinimum[{energy[ψs],ortho[ψs]==0.},{ψs,ψs0}];
ψgs=ψs/.sol;
Column[{Plot[{ψgs[[1]].Sin[π Range[M] x],ψgs[[2]].Sin[π Range[M] x]},{x,0,1},PlotRange->{-1.5,1.5},PlotRangePadding->Scaled[.03],PlotLabel->StringTemplate["\!\(g=`1`,E\_GS=`2`\)"][Round[g,1.],Round[Egs,0.01]],FrameLabel->{"\!\( x\)","\!\(ψ(x)\)"},AspectRatio->1/2,ImagePadding->{{50, 5}, {2, 5}},ImageSize->200],
Plot[2(ψgs[[1]].Sin[π Range[M] x])^2+2(ψgs[[2]].Sin[π Range[M] x])^2,{x,0,1},PlotStyle->Black,PlotRange->{0,5},PlotRangePadding->Scaled[.03],FrameLabel->{"\!\( x\)","\!\(ρ(x)\)"},AspectRatio->1/3,ImagePadding->{{50, 5}, {Automatic, 5}},ImageSize->200]},Spacings->0]]

The electrons are repelled from the infinite walls => density vanishes approaching both walls.

The electrons do not want to stay with each other (Pauli exclusion) => a small dip in the middle part of the density distribution.

Strong-interacting limit (g→∞): the electrons will further repel each other => a deeper suppress of density in the middle, breaking the distribution into two peaks.

Block[{ξ=0.25,g=50.,M=5,H0,U,ψs0,energy,ortho,Egs,sol,ψgs},
H0=SparseArray[Band[{1,1}]->Range[M]^2,{M,M}];U=Chop@Table[4g NIntegrate[Sech[(x1-x2)/ξ]Sin[m1 π x1]Sin[m2 π x2]Sin[m3 π x1]Sin[m4 π x2],{x1,0,1},{x2,0,1},AccuracyGoal->3,PrecisionGoal->2],{m1,M},{m2,M},{m3,M},{m4,M}];
ψs0=Normal@SparseArray[{{1,1}->1,{1,2}->1,{2,1}->1,{2,2}->-1},{2,M}]/Sqrt[2];
energy[ψs_List]:=ψs[[1]].H0.ψs[[1]]+ψs[[2]].H0.ψs[[2]]+tTr[{U,ψs[[1]],ψs[[2]],ψs[[1]],ψs[[2]]},{{1,5},{2,6},{3,7},{4,8}}]-tTr[{U,ψs[[1]],ψs[[2]],ψs[[2]],ψs[[1]]},{{1,5},{2,6},{3,7},{4,8}}];
ortho[ψs_List]:=Total[Flatten[ψs.Transpose[ψs]-IdentityMatrix[2]]^2];
{Egs,sol}=FindMinimum[{energy[ψs],ortho[ψs]==0.},{ψs,ψs0}];
ψgs=ψs/.sol;
Column[{Plot[{ψgs[[1]].Sin[π Range[M] x],ψgs[[2]].Sin[π Range[M] x]},{x,0,1},PlotRange->{-1.6,1.6},PlotRangePadding->Scaled[.03],PlotLabel->StringTemplate["\!\(g=`1`,E\_GS=`2`\)"][g,Round[Egs,0.01]],FrameLabel->{"\!\( x\)","\!\(ψ(x)\)"},AspectRatio->1/2,ImagePadding->{{50, 5}, {2, 5}},ImageSize->200],
Plot[2(ψgs[[1]].Sin[π Range[M] x])^2+2(ψgs[[2]].Sin[π Range[M] x])^2,{x,0,1},PlotStyle->Black,PlotRange->{0,5},PlotRangePadding->Scaled[.03],FrameLabel->{"\!\( x\)","\!\(ρ(x)\)"},AspectRatio->1/3,ImagePadding->{{50, 5}, {Automatic, 5}},ImageSize->200]},Spacings->0]]

Two localized orbitals emerges. Each electron occupies one localized orbital, and does not want to travel to the other side of the box.

For more electrons (and in higher dimensions), this will lead to crystallization of electrons => Wigner crystal.

Quantum Origin of Magnetism

Exchange Interaction (Ferromagnetism)

Antisymmetric spin state → spin singlet state (s=0)

({{s, (Subscript[m, s]),  }, {0, 0, (1/(Sqrt[2]))((Ket[↑↓])-(Ket[↓↑])).}})

Spins are anti-aligned as ⟨(Subscript[S, 1])·(Subscript[S, 2])⟩=-3(ℏ^2)/4<0 => antiferromagnetic.

Symmetric spin states → spin triplet states (s=1)

({{s, (Subscript[m, s]),  }, {1, +1, (Ket[↑↑])}, {1, 0, (1/(Sqrt[2]))((Ket[↑↓])+(Ket[↓↑]))}, {1, -1, (Ket[↓↓])}})

Spins are aligned as ⟨(Subscript[S, 1])·(Subscript[S, 2])⟩=(ℏ^2)/4>0 => ferromagnetic.

Orbital symmetric, Spin antisymmetric (antiferromagnetic)

(1/(Sqrt[2]))((Ket[(Subscript[ψ, A])(Subscript[ψ, B])])+(Ket[(Subscript[ψ, B])(Subscript[ψ, A])]))(1/(Sqrt[2]))((Ket[↑↓])-(Ket[↓↑]))

→ Interaction energy

(Subscript[E, U])=(Bra[(Subscript[ψ, A])(Subscript[ψ, B])])U(Ket[(Subscript[ψ, A])(Subscript[ψ, B])])+(Bra[(Subscript[ψ, A])(Subscript[ψ, B])])U(Ket[(Subscript[ψ, B])(Subscript[ψ, A])]).

Orbital antisymmetric, Spin symmetric (ferromagnetic)

(1/(Sqrt[2]))((Ket[(Subscript[ψ, A])(Subscript[ψ, B])])-(Ket[(Subscript[ψ, B])(Subscript[ψ, A])]))(1/(Sqrt[2]))((Ket[↑↓])+(Ket[↓↑]))

→ Interaction energy

(Subscript[E, U])=(Bra[(Subscript[ψ, A])(Subscript[ψ, B])])U(Ket[(Subscript[ψ, A])(Subscript[ψ, B])])-(Bra[(Subscript[ψ, A])(Subscript[ψ, B])])U(Ket[(Subscript[ψ, B])(Subscript[ψ, A])]).

Spin exchange interaction: the ferromagnetic spin-spin interaction between two electrons occupying two degenerated orbitals.

This magnetic (spin) interaction originates from the quantum effect of the electric (charge) interaction (the exchange energy of Coulomb interaction). The fermionic nature of electron plays a key role in this mechanism.

[If the spin interaction was originated from the magnetic dipole interaction (electromagnetic force), the spins should tend to anti-align with each other, which makes it hard to explain the existence of ferromagnets.]

Hund’s rule: atomic orbitals in the same energy shell tend to be singly occupied by electrons of the same spin (to maximize the total spin).

This is a result of the spin exchange interaction between electrons on these orbitals.

Magnetic interaction in molecules and solids. The spin exchange interaction is less effective between electrons from different atoms, because the overlap (BraKet[(Subscript[ψ, A]), (Subscript[ψ, B])]) between orbital wave functions is small, and the exchange interaction is also small in general. There is a competing effect in this case: both electrons tends to occupy the bonding orbital (Ket[(Subscript[ψ, A])])+(Ket[(Subscript[ψ, B])]) to lower the kinetic energy, then the spin must be anti-aligned due to the Pauli exclusion principle → leading to the antiferromagnetic interaction between electrons in a covalent bond.

({{mechanism, type, origin}, {covalent bond, antiferro, kinetic energy}, {spin exchange, ferro, interaction energy}})

It is quit subtle to determine whether electron spins on neighboring atoms want to be aligned or anti-aligned. Generally either behavior is possible.

Super-Exchange Interaction (Antiferromagnetism)

Consider two sites (atoms) A and B adjacent to each other, with a total of two electrons.

For each electron, there are four possible states

(Ket[(Subscript[A, ↑])]), (Ket[(Subscript[A, ↓])]), (Ket[(Subscript[B, ↑])]), (Ket[(Subscript[B, ↓])]).

Two electrons together has six possible states ((Subscript[C, 4]^2)=6)

(Ket[(Subscript[A, ↑])(Subscript[B, ↑])])=-(Ket[(Subscript[B, ↑])(Subscript[A, ↑])]),
(Ket[(Subscript[A, ↑])(Subscript[B, ↓])])=-(Ket[(Subscript[B, ↓])(Subscript[A, ↑])]),
(Ket[(Subscript[A, ↓])(Subscript[B, ↑])])=-(Ket[(Subscript[B, ↑])(Subscript[A, ↓])]),
(Ket[(Subscript[A, ↓])(Subscript[B, ↓])])=-(Ket[(Subscript[B, ↓])(Subscript[A, ↓])]),
(Ket[(Subscript[A, ↑])(Subscript[A, ↓])])=-(Ket[(Subscript[A, ↓])(Subscript[A, ↑])]),
(Ket[(Subscript[B, ↑])(Subscript[B, ↓])])=-(Ket[(Subscript[B, ↓])(Subscript[B, ↑])]).

Hopping term (Subscript[H, hop]) turns a (Subscript[A, ↑]) orbital to a (Subscript[B, ↑]) orbital and vice versa, and similarly turns a (Subscript[A, ↓]) orbital to a (Subscript[B, ↓]) orbital and vice versa (electron hopping is independent of spin).

For example, starting from (Ket[(Subscript[A, ↑])(Subscript[B, ↓])]),

(Bra[(Subscript[A, ↑])(Subscript[A, ↓])])(Subscript[H, hop])(Ket[(Subscript[A, ↑])(Subscript[B, ↓])])=-t,
(Bra[(Subscript[B, ↑])(Subscript[B, ↓])])(Subscript[H, hop])(Ket[(Subscript[A, ↑])(Subscript[B, ↓])])=-t;

starting from (Ket[(Subscript[A, ↓])(Subscript[B, ↑])]),

(Bra[(Subscript[A, ↓])(Subscript[A, ↑])])(Subscript[H, hop])(Ket[(Subscript[A, ↓])(Subscript[B, ↑])])=-t,
(Bra[(Subscript[B, ↓])(Subscript[B, ↑])])(Subscript[H, hop])(Ket[(Subscript[A, ↓])(Subscript[B, ↑])])=-t,

which can also be written as

(Bra[(Subscript[A, ↑])(Subscript[A, ↓])])(Subscript[H, hop])(Ket[(Subscript[A, ↓])(Subscript[B, ↑])])=t,
(Bra[(Subscript[B, ↑])(Subscript[B, ↓])])(Subscript[H, hop])(Ket[(Subscript[A, ↓])(Subscript[B, ↑])])=t,

given the basis relations in Eq. ().

If the electron spins are aligned, no hopping can occur. As hopping processes like

(Ket[(Subscript[A, ↑])(Subscript[B, ↑])])→(Ket[(Subscript[B, ↑])(Subscript[B, ↑])])=0,
(Ket[(Subscript[A, ↑])(Subscript[B, ↑])])→(Ket[(Subscript[A, ↑])(Subscript[A, ↑])])=0,

are forbidden by the Pauli exclusion principle.

The hopping Hamiltonian can be represented as a matrix

(Subscript[H, hop])≏(({{0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, -t, -t}, {0, 0, 0, 0, t, t}, {0, 0, 0, 0, 0, 0}, {0, -t, t, 0, 0, 0}, {0, -t, t, 0, 0, 0}})).

Interaction term (Subscript[H, int]) gives energy penalty U>0 when two electrons are on the same site (such that they repel strongly)

(Bra[(Subscript[A, ↑])(Subscript[A, ↓])])(Subscript[H, int])(Ket[(Subscript[A, ↑])(Subscript[A, ↓])])=U,
(Bra[(Subscript[B, ↑])(Subscript[B, ↓])])(Subscript[H, int])(Ket[(Subscript[B, ↑])(Subscript[B, ↓])])=U,

which can be represented as a diagonal matrix for the last two states

(Subscript[H, int])≏(({{0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, U, 0}, {0, 0, 0, 0, 0, U}})).

H≏(({{0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, -t, -t}, {0, 0, 0, 0, t, t}, {0, 0, 0, 0, 0, 0}, {0, -t, t, 0, U, 0}, {0, -t, t, 0, 0, U}})).

The spin-triplet states are eigenstates of H of the eigen energy E=0

({{energy, state}, {0, (Ket[(Subscript[A, ↑])(Subscript[B, ↑])])}, {0, (1/(Sqrt[2]))((Ket[(Subscript[A, ↑])(Subscript[B, ↓])])+(Ket[(Subscript[A, ↓])(Subscript[B, ↑])]))}, {0, (Ket[(Subscript[A, ↓])(Subscript[B, ↓])])}}).

The following double occupied (anti-bonding) state is “dark” to hopping, which is also an eigenstate of H of energy E=U,

({{energy, state}, {U, (1/(Sqrt[2]))((Ket[(Subscript[A, ↑])(Subscript[A, ↓])])-(Ket[(Subscript[B, ↑])(Subscript[B, ↓])]))}}).

The spin-singlet state and the double occupied (bonding) states can turn into each other via hopping,

(Ket[s])=(1/(Sqrt[2]))((Ket[(Subscript[A, ↑])(Subscript[B, ↓])])-(Ket[(Subscript[A, ↓])(Subscript[B, ↑])])),
(Ket[d])=(1/(Sqrt[2]))((Ket[(Subscript[A, ↑])(Subscript[A, ↓])])+(Ket[(Subscript[B, ↑])(Subscript[B, ↓])])).

The Hamiltonian projected to the (Ket[s]),(Ket[d]) basis will be represented as

H≏(({{0, -2t}, {-2t, U}})).

Strong interaction limit U≫t, the ground state (Ket[s])+𝒪(t/U) is very close to the spin-singlet state (Ket[s]) with the eigen energy

E=-((4(t^2))/U),

which is also the energy splitting between singlet and triplet state. In a spin model H=J (Subscript[S, 1])·(Subscript[S, 2]), the energy splitting between the singlet (E=-3J/4) and triplet (E=J/4) states is simply J.

Therefore Eq. () implies that the electronic system can be effectively viewed as a spin system in the strong interaction regime with an antiferromagnetic spin coupling

which is also known as the super-exchange interaction.

Weak interaction limit U≪t, the ground state is

(1/(Sqrt[2]))((Ket[s])+(Ket[d]))
=(1/2)((Ket[(Subscript[A, ↑])(Subscript[A, ↓])])+(Ket[(Subscript[A, ↑])(Subscript[B, ↓])])-(Ket[(Subscript[A, ↓])(Subscript[B, ↑])])+(Ket[(Subscript[B, ↑])(Subscript[B, ↓])]))
=(1/2)((Ket[(Subscript[A, ↑])(Subscript[A, ↓])])+(Ket[(Subscript[A, ↑])(Subscript[B, ↓])])+(Ket[(Subscript[B, ↓])(Subscript[A, ↑])])+(Ket[(Subscript[B, ↑])(Subscript[B, ↓])]))
=(1/(Sqrt[2]))((Ket[(Subscript[A, ↑])])+(Ket[(Subscript[B, ↑])]))⊗(1/(Sqrt[2]))((Ket[(Subscript[A, ↓])])+(Ket[(Subscript[B, ↓])])),

with eigen energy -2t. This describes two electrons with opposite spins occupying the single-particle bonding state together. Each electron gains the energy t as the kinetic energy is reduced by hybridizing the two atomic orbitals.

In both limits, the ground state is a spin-singlet state (antiferromagnetic). As the interaction is turned on it simply suppress the component of the double occupied state (Ket[d]) in the ground state.

Instabilities of Metallic State

Itinerant Ferromagnetism

Because creating the spin imbalance (say (Subscript[N, ↑])>(Subscript[N, ↓])) corresponds to transferring some of the ↓ spin electrons from below the Fermi energy to ↑ spin electrons above the Fermi energy, which will always increase the energy.

Block[{w=0.5,FS},
FS[y_]:={{LightGray,Rectangle[{-w,0},{0,1+y}],Rectangle[{0,0},{w,1-y}]},
{Line[{w #,0}&/@{-1,1}],{Translate[Line[{0,#}&/@{0,1.5}],{w#,0}]&/@{-1,0,1}}},
{AbsoluteThickness[1],Dashed,Line[{w #,1}&/@{-1,1}],Text["\!\(ϵ\_F\)",{w,1},{-1.5,0}]},Translate[{Arrowheads[0.06],Arrow[{{0,0},{0,1.5}}],Text["\!\( ϵ\)",{0,1.5},{2,0}]},{-1.2w,0}],
{Text["↑",{-w/2,0.5}],Text["↓",{w/2,0.5}]}};
Graphics[{Translate[FS[0],{-1,0}],{Arrowheads[0.05],Arrow[{#,0.5}&/@{-0.2,0.2}]},Translate[{FS[0.2],Arrow@BSplineCurve[{{w/2,0.9},{w/4,1.5},{-w/4,1.5},{-w/2,1.1}}]},{1,0}]},ImageSize->200]]

This can also be seen from the positive spin susceptibility χ of electron gas calculated in Eq. (),

χ=((∂M)/(∂B))=(Subscript[μ, B]^2)((2g((Subscript[ϵ, F])))/n),

meaning that the total energy (Subscript[E, 0]) of a free electron system must grow with the magnetization M as

((Subscript[E, 0])/N)=(1/(2χ))(M^2)-M B,

such that the minimal energy solution is always given by M=χ B (matching the paramagnetic behavior characterized by χ).

Consider a short-range interaction between electrons modeled by an interaction potential

U(x-x')=(Subscript[U, 0]) δ(x-x').

Electron repels (assuming g>0) each other when they are in contact at a point. This mimics the behavior of (screened) Coulomb interaction.

For electrons of the same spin,

the Hartree energy

(Subscript[E, Hartree])=(∑_((Subscript[ϵ, k]),(Subscript[ϵ, k'])<(Subscript[μ, ↑])))(Bra[(Subscript[ψ, k,↑ ])(Subscript[ψ, k',↑ ])])U(Ket[(Subscript[ψ, k,↑ ])(Subscript[ψ, k',↑ ])])
=(Subscript[U, 0])(∑_((Subscript[ϵ, k]),(Subscript[ϵ, k'])<(Subscript[μ, ↑])))∫d(Subscript[x, 1])d(Subscript[x, 2])δ((Subscript[x, 1])-(Subscript[x, 2]))((|(Subscript[ψ, k])((Subscript[x, 1]))|)^2)((|(Subscript[ψ, k'])((Subscript[x, 2]))|)^2)
=(Subscript[U, 0])(∑_((Subscript[ϵ, k]),(Subscript[ϵ, k'])<(Subscript[μ, ↑])))∫dx((|(Subscript[ψ, k])(x)|)^2)((|(Subscript[ψ, k'])(x)|)^2)
=(Subscript[U, 0])(∑_((Subscript[ϵ, k]),(Subscript[ϵ, k'])<(Subscript[μ, ↑])))(1/V)
=(((Subscript[U, 0]) (Subscript[N, ↑]^2))/V),

the Fock energy is

(Subscript[E, Fock])=(∑_((Subscript[ϵ, k]),(Subscript[ϵ, k'])<(Subscript[μ, ↑])))(Bra[(Subscript[ψ, k,↑ ])(Subscript[ψ, k',↑ ])])U(Ket[(Subscript[ψ, k',↑ ])(Subscript[ψ, k,↑ ])])
=(Subscript[U, 0])(∑_((Subscript[ϵ, k]),(Subscript[ϵ, k'])<(Subscript[μ, ↑])))∫d(Subscript[x, 1])d(Subscript[x, 2])δ((Subscript[x, 1])-(Subscript[x, 2]))(Subscript[ψ, k]^( * ))((Subscript[x, 1]))(Subscript[ψ, k']^( * ))((Subscript[x, 2]))(Subscript[ψ, k'])((Subscript[x, 1]))(Subscript[ψ, k])((Subscript[x, 2]))
=(Subscript[U, 0])(∑_((Subscript[ϵ, k]),(Subscript[ϵ, k'])<(Subscript[μ, ↑])))∫dx((|(Subscript[ψ, k])(x)|)^2)((|(Subscript[ψ, k'])(x)|)^2)
=(Subscript[U, 0])(∑_((Subscript[ϵ, k]),(Subscript[ϵ, k'])<(Subscript[μ, ↑])))(1/V)
=(((Subscript[U, 0]) (Subscript[N, ↑]^2))/V),

the Hartree energy is precisely canceled by the Fock energy, and the total interaction energy vanishes for same-spin electrons under short-range interaction,

(Subscript[E, U])=(Subscript[E, Hartree])-(Subscript[E, Fock])=0.

This can be understood as a consequence of the Pauli exclusion principle that forbids the same-spin electrons to come to the same point to interact.

For electrons of opposite spins,

the Hartree energy is similar

(Subscript[E, Hartree])=(∑_((Subscript[ϵ, k])<(Subscript[μ, ↑])))(∑_((Subscript[ϵ, k'])<(Subscript[μ, ↓])))(Bra[(Subscript[ψ, k,↑ ])(Subscript[ψ, k',↓ ])])U(Ket[(Subscript[ψ, k,↑ ])(Subscript[ψ, k',↓ ])])
=(Subscript[U, 0])(∑_((Subscript[ϵ, k])<(Subscript[μ, ↑])))(∑_((Subscript[ϵ, k'])<(Subscript[μ, ↓])))∫d(Subscript[x, 1])d(Subscript[x, 2])δ((Subscript[x, 1])-(Subscript[x, 2]))((|(Subscript[ψ, k])((Subscript[x, 1]))|)^2)((|(Subscript[ψ, k'])((Subscript[x, 2]))|)^2)
=(Subscript[U, 0])(∑_((Subscript[ϵ, k])<(Subscript[μ, ↑])))(∑_((Subscript[ϵ, k'])<(Subscript[μ, ↓])))∫dx((|(Subscript[ψ, k])(x)|)^2)((|(Subscript[ψ, k'])(x)|)^2)
=(Subscript[U, 0])(∑_((Subscript[ϵ, k])<(Subscript[μ, ↑])))(∑_((Subscript[ϵ, k'])<(Subscript[μ, ↓])))(1/V)
=(((Subscript[U, 0]) (Subscript[N, ↑])(Subscript[N, ↓]))/V).

But there is no Fock energy contribution, because the orbital state can either be symmetric or antisymmetric (as the spin state can adjust accordingly), the two possibilities have opposite signed in front of their Fock energies, which cancel out.

So the total interaction energy for opposite-spin electrons is

(Subscript[E, U])=(Subscript[E, Hartree])=(((Subscript[U, 0]) (Subscript[N, ↑])(Subscript[N, ↓]))/V).

(Subscript[E, U])=(((Subscript[U, 0]) (Subscript[N, ↑])(Subscript[N, ↓]))/V)
=((Subscript[U, 0])/(4V))(((((Subscript[N, ↑])+(Subscript[N, ↓])))^2)-((((Subscript[N, ↑])-(Subscript[N, ↓])))^2))
=(((Subscript[U, 0]) (N^2))/(4V))(1-((((M/(Subscript[μ, B]))))^2)).

N=(Subscript[N, ↑])+(Subscript[N, ↓]): total number of electrons,

M=(Subscript[μ, B])(((Subscript[N, ↑])-(Subscript[N, ↓]))/N): magnetization, average magnetic moment per electron.

(((Subscript[E, 0])+(Subscript[E, U]))/N)=(1/(2χ))(M^2)-(((Subscript[U, 0]) n)/(4 (Subscript[μ, B]^2)))(M^2)+const
=(n/(4 (Subscript[μ, B]^2)))((1/(g((Subscript[ϵ, F]))))-(Subscript[U, 0]))(M^2)+const

When ((g((Subscript[ϵ, F])))^(-1))-(Subscript[U, 0])>0, the energy minimum is at zero magnetization M=0. => The electronic system is stable in the spin-balanced state => paramagnetic metal.

When ((g((Subscript[ϵ, F])))^(-1))-(Subscript[U, 0])<0, the system can lower its energy by keep increasing the magnetization M (until some higher-order effects set in to determine the saturation value of M). => The spin-balanced electronic system is unstable towards spontaneous magnetization => ferromagnetic metal.

Spin/Charge Density Wave

H=(Subscript[H, 0])+(Subscript[H, int]),
(Subscript[H, 0])=-t(∑_i)(∑_(σ=↑,↓ ))((Subscript[c, i+1,σ]^†)(Subscript[c, i,σ])+h.c.),
(Subscript[H, int])=U(∑_i)(Subscript[n, i,↑ ])(Subscript[n, i,↓ ]).

Block[{L=3,den,latt},
den[col_,δ_,txt_,shift_]:=Block[{wave,pts},
wave[x_]:=Sum[(1+δ(-1)^i)/2Sech[5(x-i)],{i,-L-2,L+2}];
pts=Table[{x,wave[x]},{x,-L-1/2,L+1/2,0.05}];
{{col,Opacity[0.2],Polygon[Join[pts,{#,0}&/@{L+1/2,-L-1/2}]]},{col,Opacity[0.5],Line[pts]},{col,Text[Style[txt,Bold],{(1-δ)/2,0.4},{0.8shift,0}]}}];
latt[phase_]:={{Line[{#,0}&/@{-L-1/2,L+1/2}]},{White,EdgeForm[Black],Disk[{#,0},0.12]&/@Range[-L,L]},Text[phase,{L+1/2,0},{-1.1,0}]};
Graphics[{{{den[Red,0,"↑",1],den[Blue,0,"↓",-1]},latt["metal (\!\(U=0\))"]},Translate[{{den[Red,1,"↑",0],den[Blue,-1,"↓",0]},latt["SDW (\!\(U>0\))"]},{0,-1.5}],Translate[{{den[Red,1,"↑",1],den[Blue,1,"↓",-1]},latt["CDW (\!\(U<0\))"]},{0,-3}]},ImageSize->260]]

At half filling, on average there is one electron per site (a site can at most host two electrons of ↑ and ↓ spins). Without interaction, the system is in the metallic phase.

When U>0 (repulsive interaction), ↑ and ↓ spin electrons do not want to stay on the same site to avoid repulsive interaction => different spins will spontaneously occupy alternative sites => spin density wave (SDW).

When U<0 (attractive interaction), ↑ and ↓ spin electrons wants to stay on the same site to gain attractive interaction => electrons will spontaneously doubly occupy half of the sites and empty the other half alternatively => charge density wave (CDW).

Both SDW/CDW order spontaneously breaks the lattice translation symmetry and doubles the primitive unit cell.

(Subscript[H, MF])=-(u/2)(∑_i)((Ket[i+1])(Bra[i])+h.c.)-v(∑_i)(((-1))^i)(Ket[i])(Bra[i]).

u - hopping strength of the mean-field model, set to u=1 in the end ((Subscript[H, MF]) can always be rescaled by an overall factor without affecting the variational ground state).

v - stagger potential, serves as the variational parameter, controls the ordering strength.

(Ket[k,A])=(L^(-1/2))(∑_(i∈A))(e^(i k (Subscript[x, i])))(Ket[i]),
(Ket[k,B])=(L^(-1/2))(∑_(i∈B))(e^(i k (Subscript[x, i])))(Ket[i]),

A sublattice: i is even, B sublattice: i is odd.

L - number of unit cell (not using N to avoid confusion with the total electron number).

Site coordinate (Subscript[x, i])=i (unit cell volume =2).

First Brillouin zone: k∈[-π/2,π/2).

(Subscript[H, MF])=-(∑_k)(u cos k(Ket[k,A])(Bra[k,B])+h.c.)-v(∑_k)((Ket[k,A])(Bra[k,A])-(Ket[k,B])(Bra[k,B])),

H≏(Underscript[⊕, k])(Subscript[H, k]), (Subscript[H, k])≏-(({{v, u cos k}, {u cos k, -v}}))

Band dispersion

(Subscript[ϵ, k,±])=±(Sqrt[(((u cos k))^2)+(v^2)])

The stagger potential v opens a gap, splitting spectrum into two bands. Half-filling the electronic system => fully occupies the lower band.

⟨(Subscript[H, MF])⟩=(∑_(k∈BZ))(Subscript[ϵ, k,-])=-(∑_(k∈BZ))(Sqrt[(((u cos k))^2)+(v^2)])
=-L(∫_(-π/2)^(π/2))((dk)/π)(Sqrt[(((u cos k))^2)+(v^2)]).

Define the expected bonding strength (for single spin species)

γ=2⟨(Subscript[c, i+1]^†)(Subscript[c, i])+h.c.⟩,

It can be evaluated from

(Subscript[⟨((∂(Subscript[H, MF]))/(∂u))⟩, u=1])=-L γ,

therefore

γ=(∫_(-π/2)^(π/2))((dk)/π)(((cos^2) k)/(Sqrt[(cos^2) k+(v^2)])).

Block[{γ},
γ[v_?NumericQ]:=NIntegrate[Cos[k]^2/Sqrt[Cos[k]^2+v^2],{k,-π/2,π/2}]/π;
Plot[γ[v],{v,-5,5},PlotStyle->Black,PlotRange->{0,2/π},PlotRangePadding->Scaled[.02],FrameLabel->{"\!\( v\)","\!\(γ\)"},ImageSize->180]]

Define the expected density imbalance (for single spin species)

δ=⟨(Subscript[n, A])⟩-⟨(Subscript[n, B])⟩.

It can be evaluated from

(Subscript[⟨((∂(Subscript[H, MF]))/(∂v))⟩, u=1])=-L δ,

therefore

δ=(∫_(-π/2)^(π/2))((dk)/π)(v/(Sqrt[(cos^2) k+(v^2)])).

Block[{δ},
δ[v_?NumericQ]:=NIntegrate[v/Sqrt[Cos[k]^2+v^2],{k,-π/2,π/2}]/π;
Plot[δ[v],{v,-5,5},PlotStyle->Black,PlotRange->{-1,1},PlotRangePadding->Scaled[.02],FrameLabel->{"\!\( v\)","\!\(δ\)"},ImageSize->180]]

Given that ⟨(Subscript[n, A])⟩+⟨(Subscript[n, B])⟩=1 (filling one band ⇔ one fermion per unit cell),

⟨(Subscript[n, A])⟩=((1+δ)/2),⟨(Subscript[n, B])⟩=((1-δ)/2).

The expectation energy for the SDW state is

((⟨H⟩)/L)=-2t γ+U (2⟨(Subscript[n, A])⟩⟨(Subscript[n, B])⟩)
=-2t γ+(U/2)(1-(δ^2)),

which has non-trivial minimum at v≠0 when U>0.

Block[{γ,δ,U=1},
γ[v_?NumericQ]:=NIntegrate[Cos[k]^2/Sqrt[Cos[k]^2+v^2],{k,-π/2,π/2}]/π;
δ[v_?NumericQ]:=NIntegrate[v/Sqrt[Cos[k]^2+v^2],{k,-π/2,π/2}]/π;
Plot[-2 γ[v]+U/2(1-δ[v]^2),{v,-0.5,0.5},PlotStyle->Black,PlotRange->All,PlotRangePadding->Scaled[.02],PlotLabel->StringTemplate["\!\(U/t=`1`\)"][U],FrameLabel->{"\!\( v\)","\!\(⟨H⟩/(t L)\)"},ImageSize->180]]

The energy minimum is achieved at (∂_v)⟨H⟩=0 (excluding the v=0 solution), which is given by the mean-field equation

Derive Eq. ().

-2t (∂_v)γ-U δ (∂_v)δ=0.

γ+v δ=(∫_(-π/2)^(π/2))((dk)/π)(((cos^2) k+(v^2))/(Sqrt[(cos^2) k+(v^2)]))
=(∫_(-π/2)^(π/2))((dk)/π)(Sqrt[(cos^2) k+(v^2)]),

(∂_v)(γ+v δ)=(∫_(-π/2)^(π/2))((dk)/π)(v/(Sqrt[(cos^2) k+(v^2)]))=δ,

(∂_v)γ=-v (∂_v)δ.

(2t v-U δ)(∂_v)δ=0
=>2t v-U δ=0.

1-(U/(2t))(δ/v)=1-(U/(2t))(∫_(-π/2)^(π/2))((dk)/π)(1/(Sqrt[(cos^2) k+(v^2)]))=0.

For small U/t (hence small v) the equation can be approximated by

1-(U/(2t))(2/π)ln((4/(|v|)))=0,

Assuming[v>0,Series[Integrate[1/Sqrt[Cos[k]^2+v^2]/π,{k,-π/2,π/2}],{v,0,1}]]

Integrate and Expand.

whose solution is

The SDW gap |v| opens with the repulsive interaction U non-perturbatively.

Plot[If[U>0,Exp[-1/U],0],{U,-1,5},PlotStyle->Black,FrameLabel->{"\!\(U/(t π)\)","\!\(\[LeftBracketingBar]v\[RightBracketingBar]/4\)"}]

Correspondingly, the SDW order parameter δ is also induced non-perturbatively

|δ|=(t/(2U))exp(-((t π)/U)).

The expectation energy for the CDW state is

((⟨H⟩)/L)=-2t γ+U (((⟨(Subscript[n, A])⟩)^2)+((⟨(Subscript[n, B])⟩)^2))
=-2t γ+(U/2)(1+(δ^2)),

which has non-trivial minimum at v≠0 when U<0.

Block[{γ,δ,U=-1},
γ[v_?NumericQ]:=NIntegrate[Cos[k]^2/Sqrt[Cos[k]^2+v^2],{k,-π/2,π/2}]/π;
δ[v_?NumericQ]:=NIntegrate[v/Sqrt[Cos[k]^2+v^2],{k,-π/2,π/2}]/π;
Plot[-2 γ[v]+U/2(1+δ[v]^2),{v,-0.5,0.5},PlotStyle->Black,PlotRange->All,PlotRangePadding->Scaled[.02],PlotLabel->StringTemplate["\!\(U/t=`1`\)"][U],FrameLabel->{"\!\( v\)","\!\(⟨H⟩/(t L)\)"},ImageSize->180]]

The energy minimum is achieved at (∂_v)⟨H⟩=0 (excluding the v=0 solution), which is given by the mean-field equation

For small U/t (hence small v) the equation can be approximated by

1-(((-U))/(2t))(2/π)ln((4/(|v|)))=0,

whose solution is

The CDW gap |v| opens with the attractive interaction (-U) non-perturbatively.

Plot[If[U<0,Exp[1/U],0],{U,-5,1},PlotStyle->Black,FrameLabel->{"\!\(U/(t π)\)","\!\(\[LeftBracketingBar]v\[RightBracketingBar]/4\)"}]

Correspondingly, the CDW order parameter δ is also induced non-perturbatively.

|δ|=(t/(2(-U)))exp(-((t π)/((-U)))).

For example, half-filled square lattice hopping model has a square-shaped Fermi surface, where parallel edges coincide when shifted by the nesting momentum Q.

Graphics[{{EdgeForm[Black],LightGray,Polygon[{{π,0},{0,π},{-π,0},{0,-π}}]},{Red,Arrowheads[0.08],Arrow[{{-π/2,-π/2},{π/2,π/2}}],Text["Q",{π/2,π/2},{0,1.5}],Arrow[{{-π/2,π/2},{π/2,-π/2}}],Text["Q",{π/2,-π/2},{0,-1.5}]}},Frame->True,FrameLabel->{"\!\(k\_x\)","\!\(k\_y\)"},FrameTicks->{{{-π,0,π},None},{{-π,0,π},None}},ImageSize->180]

The phenomenon of two finite segments of the Fermi surface are connected by the same moment shift is called Fermi surface nesting. In this case, the Fermi surface has an instability towards SDW/CDW ordering under interaction. The density wave vector is set by the nesting momentum.

Density wave order generally leads to gap opening on the Fermi surface. At the low temperature limit (T=0), the gap size Δ scales with interaction strength U as

Δ∼exp(-(1/(2U g((Subscript[ϵ, F]))))),

where g((Subscript[ϵ, F])) is the Fermi surface density of state.