%run ../code/init_mooc_nb.ipy
randn = np.random.randn
from matplotlib import gridspec
Task 1: Simulatneous spinful time-reversal and particle-hole symmetries¶
I want to start by copying here the methods for creating Hamiltonians respecting these symmetries separately. Since I do this kind of calculation for the first time, I prefer to go step-by-step. So, we have:
def make_random_symplectic_ham(N):
if N % 2:
raise ValueError('Matrix dimension should be a multiple of 2')
sy = np.kron(np.eye(N // 2), np.array([[0, -1j], [1j, 0]]))
h = randn(N, N) + 1j * randn(N, N)
h += h.T.conj()
Th = np.dot(sy, np.dot(h.conj(), sy))
return (h + Th) / 4
np.random.seed(123)
H0 = make_random_symplectic_ham(N=4)
pprint_matrix(H0)
which is just a piece of code from the lecture, which creates a Hamiltonian with spinful TRS.
Next, the particle-hole symmetric Hamiltonian, in our case it is a BdG Hamiltonian. We will use here the recommended basis, i.e. $H = -H^{\ast}$, for completing the task:
def make_BdG_ham(N):
#This is antisymmetric basis
H = 1j * randn(2*N, 2*N)
H += H.T.conj()
return H/2
np.random.seed(123)
H0=make_BdG_ham(N=2)
pprint_matrix(H0)
which, obviously, indeed satisfies the condition described above.
Okay, now we can simply combine these symmetries and see what happens:
def make_BdG_ham_TRS(N):
#create first a BdG Hamiltonian
H = make_BdG_ham(N) #Note: it has already dimension 2N
#then we simply follow the TRS algorithm
sy = np.kron(np.eye(N), np.array([[0, -1j], [1j, 0]]))
H += H.T.conj()
TH = np.dot(sy, np.dot(H.conj(), sy))
return (H + TH) / 4
np.random.seed(123)
H0=make_BdG_ham_TRS(N=2)
pprint_matrix(H0)
Alright, we can first of all see that the two symmetries together imply zero diagonal blocks, otherwise the structure looks similar to the BdG Hamiltonian above.
We move on and want to see the spectrum now. I am a lazy person, so I will just reuse the code from the lecture again:
def plot_spectrum(alphas, spectrum, ylim=[-3.0, 4.0]):
fig, ax = plt.subplots()
ax.plot(alphas, spectrum, c='k', linewidth = 1.2)
ax.axhline(0., ls='--', linewidth=2., c='b')
ax.set_ylabel('$E$')
ax.set_xlabel(r'$\alpha$')
ax.set_xticks([0, 0.5, 1.0])
ax.set_ylim(ylim)
Now we are ready to do some plotting. We start from the spectrum:
modes = 10
np.random.seed(5)
H0 = make_BdG_ham_TRS(modes)
H1 = make_BdG_ham_TRS(modes)
ens = []
pfaffian = []
alphas = np.linspace(0, 1, 100)
for a in alphas:
H = (1 - a) * H0 + a * H1
pfaffian.append(np.real(np.sign(pf.pfaffian(1j*H))))
energies = np.linalg.eigvalsh(H)
ens.append(energies)
plot_spectrum(alphas, ens, [-3.0, 3.0])
Just to compare with the case without the TRS, we redo the same thing with BdG:
np.random.seed(5)
h0 = make_BdG_ham(modes)
h1 = make_BdG_ham(modes)
ens_noTRS = []
pfaffian_noTRS = []
for a in alphas:
h = (1 - a) * h0 + a * h1
pfaffian_noTRS.append(np.real(np.sign(pf.pfaffian(1j*h))))
energies = np.linalg.eigvalsh(h)
ens_noTRS.append(energies)
plot_spectrum(alphas, ens_noTRS, [-3.0, 3.0])
Well, the biggest difference is that there avoided crossings instead of crossings. Crossing in the latter case indicate fermion parity switches, when a ground state has to host an even/odd number of quasiprticles. In this case we do not have this feature. Why is that? Well, we learned from the TRS alone that the levels come in pairs (Kramers pairs) and thus cross the zero-energy level also in pairs. By continuing this line of discussion I would suppose that the same happens with quasiparticles in the BdG, however adding or removing quasiparticles in pairs does not alter the fermion parity and we therefore obtain equivalent states before and after the "crossing point". This is easy to understand if we remind that BdG does not conseve the number of particles and having $N$ or $N\pm 2$ is essentially the same ground state. So, at the crossing points the two states, i.e. with N and $N\pm 2$ quasiparticles, become degenerate and repell each other.
Let's see what happens to the Pfaffian invariant in this case:
def plot_spec_pf(alphas, ens, pfaffian):
fig = plt.figure(figsize=(4.5, 4.0))
gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1])
ax0 = plt.subplot(gs[0])
ax0.plot(alphas, ens, c='k', linewidth = 1.2)
ax0.axhline(0., ls='--', linewidth=2., c='b')
ax0.set_xticks([])
ax0.set_ylim(-1.5, 1.5)
ax0.set_yticks([ -1, 0, 1])
ax1 = plt.subplot(gs[1])
ax1.plot(alphas, pfaffian, color='k', linewidth=1.5)
ax1.fill_between(alphas, pfaffian, -1, color='r', alpha=0.5)
ax1.grid(True)
ax1.set_ylim(-1.2, 1.2)
ax1.set_yticks([ -1, 1])
plt.tight_layout()
ax0.set_ylabel('$E$')
ax1.set_ylabel('$Q_{BdG}$')
ax1.set_xlabel(r'$ \alpha $')
ax1.set_xticks([0, 0.5, 1.0]);
plot_spec_pf(alphas,ens, pfaffian)
plot_spec_pf(alphas,ens_noTRS, pfaffian_noTRS) #just for comparison
We can clearly see the difference: $Q_{BdG}$ is no longer giving any useful information, so no topological phase transition happens in the system.
P.S. I must admit it was a fun to do this kind of game when all the tools are already there. Thank you, the organizers, for the very nice piece of work that you do.
Task 2: Su-Schrieffer-Heeger (SSH) model¶
Alright, it was quite a fun to do the first task. Let's go to the second one. As suggested in the task, we first have to adapt the Kitaev chain Hamiltonian to the SSH Hamiltonian. So, we have:
sigma0 = np.array([[1., 0.], [0., 1.]])
sigmax = np.array([[0., 1.], [1., 0.]])
sigmay = np.array([[0., -1j], [1j, 0.]])
sigmaz = np.array([[1., 0.], [0., -1.]])
def make_ssh_chain(L=None, t1=1, t2=1):
lat = kwant.lattice.chain()
if L is None:
sys=kwant.Builder(kwant.TranslationalSymmetry((-1,)))
L = 1
infinite=True
else:
sys = kwant.Builder()
infinite=False
# transformation to antisymmetric basis
U = np.matrix([[1.0, 1.0], [1.j, -1.j]])/np.sqrt(2)
# onsite terms (every site has two orbitals connected with t1)
for x in xrange(L):
sys[lat(x)] = - t1 * U.dot(sigmax.dot(U.H))
# hopping terms (t2 connects the 2nd orbital of site n to the 1st orbital of site n+1)
if infinite:
sys[lat(1), lat(0)] = - t2 * U.dot((sigmax - 1j * sigmay).dot(U.H))/2
else:
sys[kwant.HoppingKind((1,), lat)] = - t2 * U.dot((sigmax - 1j * sigmay).dot(U.H))/2
return sys
Good! Now, as we have our system set up we can do some simulations, similarly to the lecture on Kitaev chain.
We are going to reuse the code given in the lecture:
def plot_wf(sys, wf1, wf2, ax=None, style='b-'):
"""Function to plot a wave function for the SSH model
#At zero energy there is a degeneracy. To make the graphs consistent we plot the sum of both wavefunctions absolute value squared.
#At non-zero energy use the wave function and its chiral symetric partner.
Parameters
----------
sys : The system
wf: vector of the wave function
"""
L = sys.graph.num_nodes
xs = []
for i in xrange(sys.graph.num_nodes):
xs.append(sys.sites[i].pos[0])
xs = np.array(xs)
indx = np.argsort(xs)
if ax is None:
ax = plt.subplot(111)
ax.set_xlim(-1, L)
ax.set_xlabel("$x$")
ax.set_ylabel("$|\Psi_A|^2+|\Psi_B|^2$")
temp = np.linalg.norm(wf1.reshape(-1, 2), axis=1)**2 + np.linalg.norm(wf2.reshape(-1, 2), axis=1)**2
ax.plot(xs[indx], temp[indx], style)
def plot_ssh(L, t2, fig):
#At mu=0 the first excited state is not well defined due to the massive degeneracy.
#That is why we add a small offset to mu.
sys = make_ssh_chain(L, 1, t2+1e-4)
sys = sys.finalized()
# - wave function of lowest and first excited mode(left)
# - energy diagram as function of delta (right)
# Wave functions
ham = sys.hamiltonian_submatrix()
ev, evec = np.linalg.eigh(ham)
# for plotting: the indices of the vector corrsponds to randomly
# ordered coordinates, resort those
ax = fig.add_subplot(1, 2, 2)
plot_wf(sys, evec[:, L], evec[:, L-1], ax=ax)
plot_wf(sys, evec[:, L+1] , evec[:, L-2], ax=ax, style='r--')
plt.yticks([])
plt.xticks(range(0, L, 10))
def plot_spectrum(L, t2, fig):
# Plot the eigenspectrum of the finite system
t2s = np.linspace(0.0, 4.0)
spectrum = []
for _t2 in t2s:
sys = make_ssh_chain(L, 1, _t2)
sys = sys.finalized()
ham = sys.hamiltonian_submatrix()
evs = np.linalg.eigh(ham)[0]
spectrum.append(evs)
ymax = 3.
ax = fig.add_subplot(1, 2, 1)
ax.set_ylim(-ymax, ymax)
ax.set_yticks(np.array(np.arange(-ymax, ymax+.1)))
ax.set_xlim(0.0, 4.0)
ax.set_xticks(np.arange(5))
ax.set_xlabel("$t_2/t_1$")
ax.set_ylabel("$E/t_1$")
ax.plot(t2s, spectrum, "black")
ax.plot(np.array([t2, t2]), np.array([-ymax, ymax]), 'b--')
def plot_both(L, t2):
fig = plt.figure(figsize=(9, 3.5))
plot_spectrum(L, t2, fig)
plot_ssh(L, t2, fig)
return fig
# StaticInteract precomputes the data for all the desired value before showing the plot,
# so it works well for shared notebooks.
# For 'live' calculations, you can use interact instead.
StaticInteract((lambda t2: plot_both(25, t2)),
t2=RangeWidget(0, 4.0, 0.2, name='t2', default=0, show_range=True))
# The same thing using interact looks like this:
#interact((lambda mu: plot_both(25, 1.0, 1.0, mu)),
# mu=(0, 4.0, 0.2))
We can clearly see that similarly to the case considered in the lecture, we have a topologically non-trivial phase with edge modes when $t_2>t_1$. What do we have so special in this model? Well, we have two sublattices, just like in the case of graphene. And just like in graphene 'B' sites couple only to 'A' sites, though asymmetrically to the right and to the left. So, I guess that the sublattice symmetry must be responsible for the topological properties.
On the other hand, it is useful to have a look at an SSH chain as a simple linear chain of sites, where each site has two orbitals. This has a direct similarity with the Majorana basis of the Kitaev chain model. Now, the hopping $t_1$ couples the orbitals on the same site, while the hopping $t_2$ couples the $2nd$ orbital of the site $n$ with the $1st$ orbital of site $n+1$. When $t_2=0$, we have a chain of disconnected two-orbital atoms. Energy levels are degenerate and are given by $E=\pm|t_1|$. No we start increasing $t_2$. The gap in the spectrum persists, but decreases until we reach $t_2=t_1$, where the system undergoes a topological phase transition. To understand the latter we can consider another extreme case $t_1=0$, or $t_2/t_1\rightarrow\infty$ . Then the orbitals on the same site are decoupled, but the sites themselves are coupled via $t_2$. Clearly, one can draw such a system and notice that we get two 'dangling' orbitals at the edges disconnected from the rest of the chain, which are both at $E=0$, and they are Majorana modes. Because of the sublattice symmetry, we cannot move them from $E=0$ because the spectrum has to be symmetric with respect to this energy. They can repell each other only if they overlap, but for this to happen, the energy gap $E_g\propto|t_2-t_1|$ has to close. That's why, whenever we sweep $t_2$ over the value of $t_1$ we see a phase transition from the trivial case $\left[t_1\neq 0, t_2=0\right]$ to the topologically non-trivial case $\left[t_1\neq 0, t_2>t_1\right]$.
By the way, to make connection to the intuitive picture of the gap in the Dirac spectrum going from 'positive' to 'negative' creating a Majorana mode at the domain wall we can see the following. In our case a trivial part of the space is the vacuum around the chain. Then, the domain wall for an energy gap will be at the edges of the chain when it is in a topologically non-trivial phase. Imagine that we define an energy gap, by convention, as a difference between the lower unoccupied state (at $E>0$) and the upper occupied state (at $E < 0$) of the chain, when it is in the trivial phase, $\left[t_1\neq 0, t_2=0\right]$. These two states evolve as we change the hopping $t_2$ and finally 'exchange places', meaning that the state wich was above before now shifts below the other one. Then, one could think of this situation as if an energy gap $E_g\propto(t_1-t_2)$ goes from positive value to a negative one. Then, as soon as we have a negative gap, there will be a Majorana state at the chain edges.
Let's now test whether we can break it by breaking the sublattice symmetry, which I argues was a responsible for this behavior. It is easy to break it by simply adding an extra hopping between the same-type-orbitals. We have:
def make_ssh_chain_broken_SLS(L=None, t1=1.0, t2=1.0, t_br=1.0):
lat = kwant.lattice.chain()
if L is None:
sys=kwant.Builder(kwant.TranslationalSymmetry((-1,)))
L = 1
infinite=True
else:
sys = kwant.Builder()
infinite=False
# transformation to antisymmetric basis
U = np.matrix([[1.0, 1.0], [1.j, -1.j]])/np.sqrt(2)
# onsite terms (every site has two orbitals connected with t1)
for x in xrange(L):
sys[lat(x)] = - t1 * U.dot(sigmax.dot(U.H))
# hopping terms (t2 connects the 2nd orbital of site n to the 1st orbital of site n+1)
if infinite:
sys[lat(1), lat(0)] = - t2 * U.dot((sigmax - 1j * sigmay).dot(U.H))/2
else:
#sys[kwant.HoppingKind((1,), lat)] = - t2 * U.dot((sigmax - 1j * sigmay).dot(U.H))/2
sys[kwant.HoppingKind((1,), lat)] = U.dot((-0.5*t2*(sigmax-1j*sigmay)-t_br*sigma0).dot(U.H))
return sys
def plot_ssh_broken_SLS(L, t2, t_br, fig):
#At mu=0 the first excited state is not well defined due to the massive degeneracy.
#That is why we add a small offset to mu.
sys = make_ssh_chain_broken_SLS(L, 1, t2+1e-4, t_br)
sys = sys.finalized()
# - wave function of lowest and first excited mode(left)
# - energy diagram as function of delta (right)
# Wave functions
ham = sys.hamiltonian_submatrix()
ev, evec = np.linalg.eigh(ham)
# for plotting: the indices of the vector corrsponds to randomly
# ordered coordinates, resort those
ax = fig.add_subplot(1, 2, 2)
plot_wf(sys, evec[:, L], evec[:, L-1], ax=ax)
plot_wf(sys, evec[:, L+1] , evec[:, L-2], ax=ax, style='r--')
plt.yticks([])
plt.xticks(range(0, L, 10))
def plot_spectrum_broken_SLS(L, t2, t_br, fig):
# Plot the eigenspectrum of the finite system
t2s = np.linspace(0.0, 4.0)
spectrum = []
for _t2 in t2s:
sys = make_ssh_chain_broken_SLS(L, 1, _t2, t_br)
sys = sys.finalized()
ham = sys.hamiltonian_submatrix()
evs = np.linalg.eigh(ham)[0]
spectrum.append(evs)
ymax = 3.
ax = fig.add_subplot(1, 2, 1)
ax.set_ylim(-ymax, ymax)
ax.set_yticks(np.array(np.arange(-ymax, ymax+.1)))
ax.set_xlim(0.0, 4.0)
ax.set_xticks(np.arange(5))
ax.set_xlabel("$t_2/t_1$")
ax.set_ylabel("$E/t_1$")
ax.plot(t2s, spectrum, "black")
ax.plot(np.array([t2, t2]), np.array([-ymax, ymax]), 'b--')
def plot_both_broken_SLS(L, t2, t_br):
fig = plt.figure(figsize=(9, 3.5))
plot_spectrum_broken_SLS(L, t2, t_br, fig)
plot_ssh_broken_SLS(L, t2, t_br, fig)
return fig
# StaticInteract precomputes the data for all the desired value before showing the plot,
# so it works well for shared notebooks.
# For 'live' calculations, you can use interact instead.
t_br = 1.0
StaticInteract((lambda t2: plot_both_broken_SLS(25, t2, t_br)),
t2=RangeWidget(0, 4.0, 0.2, name='t2', default=0, show_range=True))
# The same thing using interact looks like this:
#interact((lambda mu: plot_both(25, 1.0, 1.0, mu)),
# mu=(0, 4.0, 0.2))
Finally, let us plot the same thing, but as function of $t_{br}$ -- hopping term, which breaks sublattice symmetry. We will change $t_{br}$ starting from $t_{br}=0$ in a topological phase, $t_1=1, t_2=2t_1$ up to $t_{br}=4t_1$.
def plot_spectrum_broken_SLS2(L, t2, t_br, fig):
# Plot the eigenspectrum of the finite system
t_brs = np.linspace(0.0, 4.0)
spectrum = []
for _t_br in t_brs:
sys = make_ssh_chain_broken_SLS(L, 1, t2, _t_br)
sys = sys.finalized()
ham = sys.hamiltonian_submatrix()
evs = np.linalg.eigh(ham)[0]
spectrum.append(evs)
ymax = 3.
ax = fig.add_subplot(1, 2, 1)
ax.set_ylim(-ymax, ymax)
ax.set_yticks(np.array(np.arange(-ymax, ymax+.1)))
ax.set_xlim(0.0, 4.0)
ax.set_xticks(np.arange(5))
ax.set_xlabel("$t_2/t_1$")
ax.set_ylabel("$E/t_1$")
ax.plot(t_brs, spectrum, "black")
ax.plot(np.array([t_br, t_br]), np.array([-ymax, ymax]), 'b--')
def plot_both_broken_SLS2(L, t2, t_br):
fig = plt.figure(figsize=(9, 3.5))
plot_spectrum_broken_SLS2(L, t2, t_br, fig)
plot_ssh_broken_SLS(L, t2, t_br, fig)
return fig
From both simulations we can see that it is possible to still have edge states in the chain, but these states are not topologically protected. Their position is a function of $t_2$ and $t_{br}$.
# StaticInteract precomputes the data for all the desired value before showing the plot,
# so it works well for shared notebooks.
# For 'live' calculations, you can use interact instead.
t2 = 2.0
StaticInteract((lambda t_br: plot_both_broken_SLS2(25, t2, t_br)),
t_br=RangeWidget(0, 4.0, 0.2, name='t_br', default=0, show_range=True))
# The same thing using interact looks like this:
#interact((lambda mu: plot_both(25, 1.0, 1.0, mu)),
# mu=(0, 4.0, 0.2))
#I removed the result of this simulation because there is a size limit on Sage Cloud