A Chern number tells us whether something non-trivial is going on in the wavefunction and lets us distinguish between different topological phases.
Now let me clarify what I mean when I say "non-trivial".
Normally, when we are studying materials, we move from a spatial dependence for the wavefunction to a momentum dependence across a Brillouin Zone.
When we are in a non-trivial phase, we can't define a wavefunction across the entire Brillouin Zone at the same time. We can rewrite the wavefunction to cover the area that didn't work before, but then some other section isn't well-defined.
To be clear, the physics is well defined everywhere, and every way we write the wavefunction gives the same physics. The problem lies in our inability to write down a single "chart" for the whole Zone. This conundrum is similar to the problem with plotting a globe in 2 dimensions. We always have to make cuts, but the entire globe can be covered by "charts" that make up an "atlas".
A page of a Review of Modern Physics behemoth, "Classification of topological quantum matter with symmetries" [1], presented this model below. It inspired me to dig deeper and fill in the details. They attributed this model as describing the Quantum Anomolous Hall Effect, QAHE, but this form describes a wide variety topological phenomena.
\begin{equation} H(k) = R_0(k) \sigma_0 + \vec{R}(k) \cdot \vec{\sigma} \end{equation}\begin{equation} = R_0(k) \sigma_0 + R_1(k) \sigma_1 + R_2(k) \sigma_2 + R_3(k) \sigma_3 \end{equation}where $\sigma_0$ is the identity matrix and $\sigma_i$ are simply the Pauli matrices.
Combining the terms, $H(k)$ is a 2x2 matrix.
\begin{equation}=
\begin{pmatrix}
R_0 + R_3 & R_1 -i R_2 \\
R_1 + i R_2 & R_0 - R_3
\end{pmatrix}
\end{equation}
The wavefunction will then be 2-valued. This could denote two sublattices, or two different types of particles.
While we could use any values for $\vec{R}(k)$, we will use \begin{equation} \vec{R}(k) = \begin{pmatrix} -2 \sin k_x \\ -2 \sin k_y \\ \mu +2 \sum_{x,y} \cos k_i \\ \end{pmatrix} \end{equation} as it's a fairly simple form that gives us the physics we want and exhibits phase transitions depending on the value of $\mu$. I set $\mu$ as $1$ early on. Go through the notebook with that value, then go through the notebook with $\mu = -5, -1, 1, 5$.
In [1]:
using Plots
pyplot()
Out[1]:
In [2]:
μ=1
labels=["-π","-π/2","0","π/2","π"]
ticks=[-π,-π/2,0,π/2,π];
ks=range(-π,stop=π,length=314)
l=length(ks)
ka=Array{Array{Float64},2}(undef,l,l)
for ii in 1:l
x=ks[ii]
for jj in 1:l
ka[ii,jj]=[x,ks[jj]]
end
end
In [3]:
function R0(k::Array)
return 0
end
function R1(k::Array)
return -2*sin(k[1])
end
function R2(k::Array)
return -2*sin(k[2])
end
function R3(k::Array)
return μ+2*cos(k[1])+2*cos(k[2])
end
Out[3]:
In [4]:
function R(k::Array)
return sqrt(R1(k)^2+R2(k)^2+R3(k)^2)
end
Out[4]:
First, let's just take a look at the energy spectrum. For each $k$, we calculate the eigenvalues of the 2x2 matrix Hamiltonian.
To make the calculation simpler, I denote the function $R$ as \begin{equation} R=\sqrt{R_1^2+R_2^2+R_3^2} \end{equation} With that, energy can simply be written as \begin{equation} \lambda= R_0 \pm R \end{equation} $R_0$ moves the entire spectrum up and down but doesn't effect the gap, and won't affect the physics. That term won't even factor into the eigenvectors.
In [5]:
function λp(k::Array{Float64})
return R0(k)+R(k)
end
function λm(k::Array{Float64})
return R0(k)-R(k)
end
Out[5]:
Notation Note: I denote the Array evaluation of a function as the function name followed by 'a'.
In [6]:
λpa=zeros(Float64,l,l)
λma=zeros(Float64,l,l)
for ii in 1:l
for jj in 1:l
λpa[ii,jj]=λp(ka[ii,jj])
λma[ii,jj]=λm(ka[ii,jj])
end
end
In [7]:
surface(ks,ks,λpa)
surface!(ks,ks,λma)
plot!(xticks= (ticks,labels),
yticks=(ticks,labels),
xlabel="kx",
ylabel="ky",
zlabel="Energy",
title="Band Diagram")
Out[7]:
To find the eigenvectors, we find the nullspace of the following matrix, \begin{equation} \begin{vmatrix} R_0 +R_3 - R_0 \mp R & R_1 - i R_2 \\ R_1 + i R_2 & R_0 - R_3 -R_0 \mp R \\ \end{vmatrix} \end{equation}
Take the bottom row. Second column value multiples first column and first column value multiples second column. \begin{equation} \begin{pmatrix} \pm R + R_3 \\ R_1 + i R_2 \end{pmatrix} \end{equation}
And now we need to normalize the state \begin{equation} \frac{1}{\sqrt{2 R \left( R \pm R_3 \right) }} \begin{pmatrix} \pm R + R_3\\ R_1 + i R_2 \end{pmatrix} \end{equation} In low energy, just the minus state will be occupied (lower energy), but it has a singularity at $\vec{R} = ( 0 ,0, R) $.
We have two complex numbers in this vector, and thus four things to plot $r_1, \theta_1, r_2, \theta_2$, \begin{equation} |\psi \rangle = \begin{pmatrix} r_1 e^{i \theta_1} \\ r_2 e^{i \theta_2} \end{pmatrix} \end{equation}
If we chose to create our vector from the top row instead, \begin{equation} \begin{pmatrix} R_1 -i R_2 \
We can move where the singularity is, but we can't get rid of it. The problem-point doesn't show up in the physics, only in the wavefunction representation of it. We cannot well represent the state across the entire Brillouin Zone at the same time. This problem occurs because we have a topologically non-trivial state and a non-zero Chern number.
In [8]:
function up1(k::Array)
front=1/sqrt(2*R(k)*(R(k)+R3(k)))
return front*(R(k)+R3(k))
end
function up2(k::Array)
front=1/sqrt(2*R(k)*(R(k)+R3(k)))
return front*(R1(k)+im*R2(k))
end
up=Function[]
push!(up,up1)
push!(up,up2)
Out[8]:
In [9]:
function um1(k::Array)
front=1/sqrt(2*R(k)*(R(k)-R3(k)))
return front*(-R(k)+R3(k))
end
function um2(k::Array)
front=1/sqrt(2*R(k)*(R(k)-R3(k)))
return front*(R1(k)+im*R2(k))
end
um=Function[]
push!(um,um1)
push!(um,um2)
Out[9]:
In [10]:
upa=zeros(Complex{Float64},2,l,l)
uma=zeros(Complex{Float64},2,l,l)
for ii in 1:l
for jj in 1:l
upa[1,ii,jj]=up[1](ka[ii,jj])
upa[2,ii,jj]=up[2](ka[ii,jj])
uma[1,ii,jj]=um[1](ka[ii,jj])
uma[2,ii,jj]=um[2](ka[ii,jj])
end
end
In [11]:
# Plotting θ2
heatmap(ks,ks,angle.(uma[2,:,:])-angle.(uma[1,:,:]))
plot!(xticks= (ticks,labels),
yticks=(ticks,labels),
xlabel="kx",
ylabel="ky",
title="difference in phase, lower energy")
Out[11]:
In [12]:
surface(ks,ks,abs2.(uma[2,:,:]))
plot!(xticks= (ticks,labels),
yticks=(ticks,labels),
xlabel="kx", ylabel="ky",
title="magnitude,second component, lower energy")
Out[12]:
In [13]:
surface(ks,ks,abs2.(uma[1,:,:]))
plot!(xticks= (ticks,labels),
yticks=(ticks,labels),
xlabel="kx",
ylabel="ky",
title="magnitude,first component, lower energy")
Out[13]:
The first step in calculating the Chern number is evaluating the Berry Connection.
\begin{equation} \mathcal{A}^{i} = \langle u (k,r) | d_{i} u (k,r) \rangle \end{equation}Though $\mathcal{A}^i$ looks like a vector, it is not invariant under gauge tranformation. If a wavefunction transforms as \begin{equation} | u(k,r) \rangle \rightarrow e^{-i \phi} | u(k,r) \rangle \end{equation} then the connection transforms as \begin{equation} \mathcal{A}^i \rightarrow \mathcal{A}^i -i \partial_i \phi \end{equation}
Near the Dirac points in the Brillouin Zone, the wavefunction has quite a high curvature, which makes the numerical computation of the derivative prone to errors. I tried numerical derivation but was unable to arrive at a correct, stable answer. Therefore, I'm using the 'ForwardDiff' Package to take derivatives analytically.
In [14]:
using ForwardDiff
Before dealing with the physics, let's just look at the syntax for the package.
Here's just an example of taking the gradient of $x^2$.
In [15]:
ex=x->ForwardDiff.gradient(t->t[1]^2,x)
ex([1])
Out[15]:
Now let's apply that syntax to our wavefunctions.
The 'ForwardDiff' package seems to only work on purely real functions, so we have to take the derivative of the real and imaginary parts separately.
In [16]:
dum1=kt->ForwardDiff.gradient(um1,kt)
Rdum2=kt->ForwardDiff.gradient(t->real(um2(t)),kt)
Idum2=kt->ForwardDiff.gradient(t->imag(um2(t)),kt)
dum2(k)=Rdum2(k)+im*Idum2(k)
Out[16]:
With the derivatives, we can now calculate the connection.
In [17]:
Amkx(k)=conj(um1(k))*dum1(k)[1]+conj(um2(k))*dum2(k)[1]
Amky(k)=conj(um1(k))*dum1(k)[2]+conj(um2(k))*dum2(k)[2]
Out[17]:
In [18]:
Akxa=Array{Complex{Float64}}(undef,l,l)
Akya=Array{Complex{Float64}}(undef,l,l)
for ii in 1:l
for jj in 1:l
Akxa[ii,jj]=Amkx(ka[ii,jj])
Akya[ii,jj]=Amky(ka[ii,jj])
end
end
In [19]:
heatmap(ks,ks,log.(abs.(Akxa)))
plot!(xticks= (ticks,labels),
yticks=(ticks,labels),
xlabel="kx",ylabel="ky",
title="A_kx Magnitude Berry Connection")
Out[19]:
In [20]:
heatmap(ks,ks,log.(abs.(Akya)))
plot!(xticks= (ticks,labels),
yticks=(ticks,labels),
xlabel="kx",ylabel="ky",
title="A_ky Magnitude Berry Connection")
Out[20]:
In [21]:
heatmap(ks,ks,angle.(Akxa))
plot!(xticks= (ticks,labels),
yticks=(ticks,labels),
xlabel="kx",ylabel="ky",
title="A_kx Berry Connection angle")
Out[21]:
In [22]:
heatmap(ks,ks,angle.(Akya))
plot!(xticks= (ticks,labels),
yticks=(ticks,labels),
xlabel="kx",ylabel="ky",
title="A_ky Berry Connection angle")
Out[22]:
I previously had a section here about Berry phases, but I need to fix it. Until then, please believe the mass of physics literature that from the guage-dependent Berry connection, we can get a gauge independent quantity called the Berry Curvature. Hopefully I can clarify this mathematical magic when I understand it myself.
Of course it's called Berry. Every thing in this calculation is named after Berry. At least its not Euler. Done with side note, back to equations. Here's our version of Gauss's Theorem: \begin{equation} \oint \mathcal{A} \cdot ds = \iint F_{n}^{xy} d^2 k \end{equation}
And $F_n^{xy}$ will take on this form: \begin{equation} F_{n}^{xy} = \nabla \times \mathcal{A} = \partial_{k_x} \mathcal{A}^y_{n} - \partial_{k_y} \mathcal{A}^x_{n} \end{equation} \begin{equation} = \partial_{k_x} \langle u | \partial_{k_y} u \rangle - \partial_{k_y} \langle u | \partial_{k_x} u \rangle \end{equation}
TIME FOR MORE DERIVATIVES!
In [23]:
DRAmkx=kt->ForwardDiff.gradient(t->(real(Amkx(t))),kt )
DImAmkx=kt->ForwardDiff.gradient(t->(imag(Amkx(t))),kt )
DRAmky=kt->ForwardDiff.gradient(t->(real(Amky(t))),kt )
DImAmky=kt->ForwardDiff.gradient(t->(imag(Amky(t))),kt )
Out[23]:
In [24]:
F(k)=DRAmky(k)[1]+im*DImAmky(k)[1]-DRAmkx(k)[2]-im*DImAmkx(k)[2]
Out[24]:
In [25]:
Fa=Array{Complex{Float64}}(undef,l,l)
for ii in 1:l
for jj in 1:l
Fa[ii,jj]=F(ka[ii,jj])
end
end
In [26]:
heatmap(ks,ks,abs.(Fa))
plot!(xticks= (ticks,labels),
yticks=(ticks,labels),
xlabel="kx",ylabel="ky",
title="F Berry curvature norm")
Out[26]:
In [27]:
heatmap(ks,ks,angle.(Fa))
plot!(xticks= (ticks,labels),
yticks=(ticks,labels),
xlabel="kx",ylabel="ky",
title="F Berry Curvature phase")
Out[27]:
Last Step!
Easy calculation now in terms of code, but this number has a great deal of significance. I'm still trying to wrap my head around it. The Chern number seems to pop up in a variety of obscure mathematical stuff over this physicist's head, but hopefully none of that is necessary to grasp its incredible mind-blowing usefullness.
This single integer not only seperates out topological phases from topologically trivial phases, but seperates different topological phases from each other. And always evaluates to an integer.
\begin{equation} Ch = \frac{1}{2 \pi i} \iint_{BZ} F^{xy}_n d^2 k \end{equation}If you don't get approximately an integer, try a finer mesh.
In [28]:
sum(Fa)*(ks[2]-ks[1])^2/(2π*im)
Out[28]:
You made it!
Even if you didn't understand everything, or really understand much at all, pat yourself on the back. This is the very frontier of science.
[1] Ching Kai Chiu, Jeffrey C.Y. Teo, Andreas P. Schnyder, and Shinsei Ryu, Classification of topological quantum matter with symmetries, Reviews of Modern Physics 88 2016, no. 3, 1–70.
@article{RevelsLubinPapamarkou2016, title = {Forward-Mode Automatic Differentiation in Julia}, author = {{Revels}, J. and {Lubin}, M. and {Papamarkou}, T.}, journal = {arXiv:1607.07892 [cs.MS]}, year = {2016}, url = {https://arxiv.org/abs/1607.07892} }
In [ ]: