#!/usr/bin/env python
# coding: utf-8

# # The Su-Schrieffer–Heeger (SSH) model

# Saumya Biswas (saumyab@uoregon.edu)

# The celebrated SSH model is analyzed with QuTiP's lattice module below.

# ![title](images/SSH_E.png)

# The above figure shows a SSH model with 6 sites with periodic boundary condition. The same lattice with hardwall/aperiodic boundary condition would be the folloowing.
# ![title](images/SSH_Ea.png)

# In the secod quantized formalism, the periodic lattice Hamiltonian can be written as
# \begin{eqnarray}
# H_{per} = -t_{intra} (c_{-2}^{\dagger} c_{-1} + c_{-1}^{\dagger} c_{-2} ) -t_{intra} (c_{0}^{\dagger} c_{1} + c_{1}^{\dagger} c_{0} ) -t_{intra} (c_{2}^{\dagger} c_{3} + c_{3}^{\dagger} c_{2} ) \nonumber \\ 
#  -t_{inter} (c_{-1}^{\dagger} c_{0} + c_{0}^{\dagger} c_{-1} ) -t_{inter} (c_{1}^{\dagger} c_{2} + c_{2}^{\dagger} c_{1} ) -t_{inter} (c_{3}^{\dagger} c_{-2} + c_{-2}^{\dagger} c_{3} ) \nonumber 
# \end{eqnarray}
# The aperiodic lattice Hamiltonian can be obtained by discarding the very last term.
# \begin{eqnarray}
# H_{aper} = -t_{intra} (c_{-2}^{\dagger} c_{-1} + c_{-1}^{\dagger} c_{-2} ) -t_{intra} (c_{0}^{\dagger} c_{1} + c_{1}^{\dagger} c_{0} ) -t_{intra} (c_{2}^{\dagger} c_{3} + c_{3}^{\dagger} c_{2} ) \nonumber \\ 
#  -t_{inter} (c_{-1}^{\dagger} c_{0} + c_{0}^{\dagger} c_{-1} ) -t_{inter} (c_{1}^{\dagger} c_{2} + c_{2}^{\dagger} c_{1} ) \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \nonumber 
# \end{eqnarray}

# The representation in terms of unit cell blocks become obvious once we resolve the terms into unit cell operators.
# \begin{eqnarray}
# H_{per}= \begin{bmatrix}
#     c_{-2}^{\dagger} & c_{-1}^{\dagger}       
# \end{bmatrix} 
# \begin{bmatrix}
#       o     & -t_1 \\
#     -t_1    & 0 
# \end{bmatrix}
# \begin{bmatrix}
#     c_{-2}  \\
#     c_{-1}  
# \end{bmatrix} + 
# \begin{bmatrix}
#     c_{0}^{\dagger} & c_{1}^{\dagger}       
# \end{bmatrix} 
# \begin{bmatrix}
#       o     & -t_1 \\
#     -t_1    & 0 
# \end{bmatrix}
# \begin{bmatrix}
#     c_{0}  \\
#     c_{1}  
# \end{bmatrix} 
# \nonumber \\
# + \begin{bmatrix}
#     c_{2}^{\dagger} & c_{3}^{\dagger}       
# \end{bmatrix} 
# \begin{bmatrix}
#       o     & -t_1 \\
#     -t_1    & 0 
# \end{bmatrix}
# \begin{bmatrix}
#     c_{2}  \\
#     c_{3}  
# \end{bmatrix} \nonumber \\
# + \begin{bmatrix}
#     c_{-2}^{\dagger} & c_{-1}^{\dagger}       
# \end{bmatrix} 
# \begin{bmatrix}
#       o     & 0 \\
#     -t_2    & 0 
# \end{bmatrix}
# \begin{bmatrix}
#     c_{0}  \\
#     c_{1}  
# \end{bmatrix} + 
# \begin{bmatrix}
#     c_{0}^{\dagger} & c_{1}^{\dagger}       
# \end{bmatrix} 
# \begin{bmatrix}
#       o     & -t_2 \\
#       0     & 0 
# \end{bmatrix}
# \begin{bmatrix}
#     c_{-2}  \\
#     c_{-1}  
# \end{bmatrix} 
# \nonumber \\
# + \begin{bmatrix}
#     c_{0}^{\dagger} & c_{1}^{\dagger}       
# \end{bmatrix} 
# \begin{bmatrix}
#       o     & 0 \\
#     -t_2    & 0 
# \end{bmatrix}
# \begin{bmatrix}
#     c_{2}  \\
#     c_{3}  
# \end{bmatrix} + 
# \begin{bmatrix}
#     c_{2}^{\dagger} & c_{3}^{\dagger}       
# \end{bmatrix} 
# \begin{bmatrix}
#       o     & -t_2 \\
#       0     & 0 
# \end{bmatrix}
# \begin{bmatrix}
#     c_{0}  \\
#     c_{1}  
# \end{bmatrix} 
# \nonumber \\
# + \begin{bmatrix}
#     c_{2}^{\dagger} & c_{3}^{\dagger}       
# \end{bmatrix} 
# \begin{bmatrix}
#       o     & 0 \\
#     -t_2    & 0 
# \end{bmatrix}
# \begin{bmatrix}
#     c_{-2}  \\
#     c_{-1}  
# \end{bmatrix} + 
# \begin{bmatrix}
#     c_{-2}^{\dagger} & c_{-1}^{\dagger}       
# \end{bmatrix} 
# \begin{bmatrix}
#       o     & -t_2 \\
#       0     & 0 
# \end{bmatrix}
# \begin{bmatrix}
#     c_{2}  \\
#     c_{3}  
# \end{bmatrix} \nonumber
# \end{eqnarray}
# 

# Hence, $H_{TB}$ can be succinctly put in the form:
# \begin{eqnarray}
# H_{per/aper} = \sum_i \psi_i^{\dagger} D \psi_i + \sum_{i} \left( \psi_i^{\dagger} T \psi_{i+1} +  \psi_{i+1}^{\dagger} T^{\dagger} \psi_i  \right) \label{eq:TB_block} 
# \end{eqnarray}
# where $D = \begin{bmatrix}
#     0 & -t_{intra}  \\
#     -t_{intra} & 0  
# \end{bmatrix}$ and $T = \begin{bmatrix}
#     0 & 0  \\
#     -t_{inter} & 0  
# \end{bmatrix}$. The $\psi_i = \begin{bmatrix}
#     c^0_{i}  \\
#     c^1_{i}  
# \end{bmatrix}$ is associated with $\vec{R}=\vec{R}_i$ and $\psi_{i+1}$ is associated with $\vec{R}=\vec{R}_i + \hat{x}$, with $\hat{x}$ being the unit vector along the x direction.
# 
# The equation above can be put into the alternate form (also changing the summation variables to m,n to distinguish them from the imaginary i):
# 
# \begin{eqnarray}
# H_{per/aper} = \sum_{m,n} \psi_m^{\dagger} (D \delta_{m,n} + T \delta_{m,n-1} + T^{\dagger} \delta_{m,n+1} ) \psi_n \label{eq:TB_block_1} 
# \end{eqnarray}
# 
# 

# So far, we have see that a SSH lattice can be put into the unit cell representations and the lattice can be completely defied with two matrices D and T. In declaring an instance of Qutip.lattice.Lattice1d we only need to input the two matrices. 

# In declaring an instance of Lattice1d class, the two arguments cell_Hamiltonian
# and inter_hop are set to D and T respectively. The matrix structure of D and T can be
# obtained from the aide function cell_structures().

# In[1]:


from qutip import *
import matplotlib.pyplot as plt
import numpy as np


# In[2]:


val_s = ['site0','site1']
(H_cell_form,T_inter_cell_form,H_cell,T_inter_cell) = cell_structures( val_s)


# In[3]:


H_cell_form


# In[4]:


T_inter_cell_form


# Guided by cell_H_form and inter_cell_T_form, we can set values to cell_H and inter_cell_T which were initialized to all zero elements by cell_structures().

# In[5]:


t_intra = -0.5
t_inter = -0.5
H_cell[0,1] = t_intra
H_cell[1,0] = t_intra
T_inter_cell[1,0] = t_inter
H_cell = Qobj(H_cell)
T_inter_cell = Qobj(T_inter_cell)


# Using cell_structures() is completely optional. The user could equally well have defined cell_H and inter_cell_T directly.

# In[6]:


H_cell = Qobj( np.array( [[ 0, t_intra ],[t_intra,0]] )  )
T_inter_cell =  Qobj( np.array( [[ 0, 0 ],[t_inter,0]] )  ) 


# For our SSH lattice with 3 unit cells, 2 sites in each unit cell, and [1] degree of freedom per each site, we can initiate an instance of the Lattice1d class at this stage. 

# In[7]:


boundary_condition = "periodic"
cells = 3
cell_sites = 2
site_dof = [1]
SSH_lattice = Lattice1d(num_cell=cells, boundary = boundary_condition, cell_num_site = cell_sites, cell_site_dof = site_dof, Hamiltonian_of_cell = H_cell, inter_hop = T_inter_cell )


# The model can be visualized with the display functions.

# In[8]:


H = SSH_lattice.display_unit_cell(label_on = True)
T = SSH_lattice.display_lattice()


# The Hamiltonian of the lattice can be obtained with the method Hamiltonian()

# In[9]:


SSH_Haml = SSH_lattice.Hamiltonian()
SSH_Haml


# ### Sublattice Projectors and Chiral Symmetry of the SSH model:
# We explain now that, the Hamiltonian of the SSH model is bipartite. We can set the sites into two partitions. Each partition consists of the set of every other sites starting from site 0(or 1). Since, the Hamiltonian  does not enable any transition between sites of the same sublattice, only between them, the Hamiltonian is said to be bipartite.

# \begin{eqnarray}
# \hat{P}_{A} = \sum\limits_{m=1}^{N} |m,A\rangle \langle n, A|,\ \ \ \ \ \ \hat{P}_{B} = \sum\limits_{m=1}^{N} |m,B\rangle \langle n, B|, \ \ \ \ \ \
# \hat{\Sigma}_z = \hat{P}_{A} - \hat{P}_{B}
# \end{eqnarray}

# The projectors $\hat{P}_A$(or $\hat{P}_B$) project out a ket vector into its projection into sublattice A(or B). Because, the Hamiltonian is bipartite,
# \begin{eqnarray}
# \hat{P}_{A} H \hat{P}_{A} = \hat{P}_{B} H \hat{P}_{B} = 0
# \end{eqnarray}
# 

# The SSH Hamiltonian is said to have chiral symmetry, since
# \begin{eqnarray}
# \hat{\Sigma}_z H \hat{\Sigma}_z = -H
# \end{eqnarray}
# which is a generalization of the symmetry concept of a Hamiltonian, since, ordinarily, a Hamiltonian is said to have a symmetry if an unitary operator leaves it invariant.
# \begin{eqnarray}
# \hat{U} H \hat{U}^{\dagger} = H
# \end{eqnarray}
# Now, we form the operators $\hat{P}_A$,$\hat{P}_B$ and $\hat{\Sigma}_z$, and verify the properties of the SSH Hamiltonian.

# In[10]:


chiral_op = SSH_lattice.distribute_operator(sigmaz())        
anti_commutator_chi_H = chiral_op * SSH_Haml + SSH_Haml * chiral_op
is_null = (np.abs(anti_commutator_chi_H) < 1E-10).all()
print(is_null)


# Hence, it is verified that $\hat{\Sigma}_z$ and H indeed anticommute and the SSH Hamiltonian has chiral symmetry.

# The dispersion relationship for the lattice can be obtained with plot_dispersion() method.

# In[11]:


SSH_lattice.plot_dispersion()


# plot_dispersion() plots the 3(since number of unit cells in 3) points in k-space 
# (the first and last one are the same) over the dispersion relation of an infinite
# crystal. 

# In[12]:


[V,[S0,S1,S2,S3,S4,S5]]=SSH_Haml.eigenstates()


# In[13]:


V


# First, the eigen-values are the same as the ones obtained from the dispersion calculation. Second, they are symmetric about the value 0.
# The second is a consequence of the chiral symmetry of the Hamiltonian, as we explain now.
# \begin{eqnarray}
# \hat{\bf{1}} = \hat{P}_A + \hat{P}_B, \ \ \ \ \, \hat{P}_A = \frac{1}{2}(\hat{\bf{1}}+\hat{\Sigma}_z), \ \ \ \ \, \hat{P}_B = \frac{1}{2}(\hat{\bf{1}}-\hat{\Sigma}_z) \nonumber \\
# H |\psi_n\rangle = E_n | \psi_n \rangle \implies H (\hat{\Sigma}_z | \psi_n \rangle) = - \hat{\Sigma}_z H | \psi_n \rangle = -E_n (\hat{\Sigma}_z|\psi_n\rangle)
# \end{eqnarray}
# So, if $|\psi_n\rangle$ is an eigenstate with eigenenergy $E_n$, $\hat{\Sigma}_z | \psi_n \rangle$ is also an eigenstate with energy $-E_n$ and the eigen-spectrum is symmetric about 0.
# 

# Here, S0 is the eigenvector with eigenvalue -1 and S5 is the eigenvector with eigenvalue +1. So, we can verify if S5 is the same eigenvector(withinn a phase factor) as ($\hat{\Sigma}_z*$S0).

# In[14]:


print(S0)
print(chiral_op*S0)
print(S5)


# Clearly, S5 is the same eigenvector as ($\hat{\Sigma}_z*$S0).

# Since, $\hat{\Sigma}_z | \psi_n \rangle$ and $| \psi_n \rangle$ are eigenvectors of a Hermitian oerator,H with distinct eigenvalues, they must be orthogonal.
# \begin{eqnarray}
# E_n \ne 0 \implies 0 = \langle \phi_n| \hat{\Sigma}_z | \phi_n \rangle = \langle \phi_n \hat{P}_A |\phi_n \rangle - \langle \phi_n| \hat{P}_B |\phi_n \rangle
# \end{eqnarray}
# i.e., an eigenstate with nonzero eigenvalue has equal support on both sublattices. We, now check this for S5.
# 

# In[15]:


is_null = (np.abs(S0.dag()*S5) < 1E-10).all()
print(is_null) # Are S0 and S5 orthogonal?


# In[16]:


dimH = chiral_op.dims
identity_H = Qobj(np.identity(6), dims=dimH)
identity_H
PA = 0.5*( identity_H + chiral_op)
PB = 0.5*( identity_H - chiral_op)
support_S5_A = S5.dag() * PA * S5
support_S5_B = S5.dag() * PB * S5
print(support_S5_A == support_S5_B) # Does S5 have equal support on A and B?


# We discuss the implications of $E_n = 0$ for an eigenvector later in the context of edge states later.

# Unsurprisingly, diagonalizing the Hamiltonian gives the same spectrum of eigen-values
# as the one obtained from the plot_dispersion() function. We shall soon illustrate, translational symmetry is a very useful computational hack. Here, we see how they can produce the eigen-values and eigen-vectors of the Hamiltonian ($6\times6$, in our example) from diagonalizing a $2\times2$ matrix. The reduction in size by a factor of 3 comes from the fact that the lattice of 3 cells repeats itself infinitely on both ends.

# ### Using Translational Symmetry:
# Any periodic lattice Hamiltonian can be diagonalized easier exploiting translational symmetry, this feature is not specific to SSH model alone.

# \begin{eqnarray}
# |\psi_n(k) \rangle = |k \rangle \otimes | u_{n}(k) \rangle   \nonumber \\
# | u_{n}(k) \rangle = a_n(k)|a\rangle  + b_n(k)|b\rangle \nonumber \\
# \end{eqnarray}
# The vectors $| u_{nk}(r) \rangle \in H_{internal}$ are the eigenstates of the bulk momentum space Hamiltonian $H(k)$ defined as
# \begin{eqnarray}
# \langle k | H_{bulk} | k \rangle = \sum\limits_{\alpha,\beta \in {A,B}} \langle k, \alpha | H_{bulk} | k, \beta \rangle | \alpha \rangle \langle \beta |  \nonumber \\
# H(k)|u_{n}(k) \rangle = E(k)| u_{n}(k) \rangle
# \end{eqnarray}
# In a lattice with N unitcells, since $|\psi_n(k) \rangle$ is required to be periodic over N cells, it needs to be invariant with a translation by N cells, and the valid/good quantum number for k are $0,1,2,...,\frac{2\pi}{N}$.
# 

# The dispersion, $E(k)$ can be obtained from get_dispersion() method. A list of $H(k)$ at the valid quantum numbers can be produced from the bulk_Hamiltonian_array() method.

# In[17]:


(knxA, qH_ks) = SSH_lattice.bulk_Hamiltonians()
qH_ks


# The array of $|u_n(k) \rangle$(in terms of its expansion in {a(k),b(k)}) at the good quantum numbers, k can be produced with the method array_of_unk().

# In[18]:


(knxA, vec_kns) = SSH_lattice.cell_periodic_parts()
vec_kns


# In both cases, knxA is simply an array containing the valid values of k, in units of $2\pi/a$, a being the length of the unit cell i.e. 1. a is always 1 in all methods().

# In[19]:


knxA


# bloch_wave_functions() yields an ordered array for the eigenvalues and eigenvectors(which are bloch wave functions) for the Hailtonian of the lattice.

# In[20]:


eigen_states = SSH_lattice.bloch_wave_functions()
eigen_states


# Knowing the cell periodic part of a bloch wavefunction, $u_n(k)$ suffices to calculate it. The translational symmetry enables calculation of the eigenstates of a $6\times6$ matrix through diagonalizing a $2\times2$ matrix.

# ### Topology: Winding Number:
# Due to the chiral symmetry, the bulk momentum space Hamiltonian,$H(k)$ can be written in terms of $\sigma_x$ and $\sigma_y$ components alone.
# \begin{eqnarray}
# H(k)= h_x(k)\sigma_x + h_y(k)\sigma_y
# \end{eqnarray}
# For this specific model, where $\bf{h}(k)$ moves about in a 2d plane($h_x-h_y$), winding number is a topological invariant characterizing the topology of the model. It enumerates the number of times, $\bf{h}(k)$ traverses around the origin in a positive sense as k is varied from 0 to $2\pi$.
# The method winding_number() evaluates the following integral to determine it as well as plots the trajectory of $\bf{h}(k)$ in the $h_x-h_y$ plane.
# \begin{eqnarray}
# \nu = \frac{1}{2\pi i}\int\limits_{-\pi}^{\pi} dk \frac{d}{dk}Log(h(k)) \ \ \ \text{,where} \ \ \ \ \ \ \ \ h(k) = h_x(k)- ih_y(k)
# \end{eqnarray}
# 
# 
# 
# 
# 

# ### The trivial Insulator:
# The latice become a gapped system when $t_{inter}$ does not equal $t_{intra}$. For the case of $|t_{inter}| < |t_{intra}|$, the lattice is a topologically trivial insulator.

# In[22]:


t_intra = -0.5
t_inter = -0.35
H_cell = Qobj( np.array( [[ 0, t_intra ],[t_intra,0]] )  )
T_inter_cell =  Qobj( np.array( [[ 0, 0 ],[t_inter,0]] )  ) 
SSH_lattice_TrI = Lattice1d(num_cell=100, boundary = "periodic", cell_num_site = 2, cell_site_dof = [1], Hamiltonian_of_cell = H_cell, inter_hop = T_inter_cell )


# In[23]:


SSH_lattice_TrI.plot_dispersion()


# In[24]:


SSH_H_t = SSH_lattice_TrI.Hamiltonian()
D = SSH_H_t.eigenenergies()


# In[25]:


plt.plot(D,'ro')
plt.xlabel('index of eigen values')
plt.ylabel('eigen values')
plt.show()
plt.close()


# Unsurprisingly, we see a gap in the dispersion relationship and also in the spectrum of
# eigen values of the Hamiltonian.

# However, a calculation of the winding number yields 0, since it has trivial topology.

# In[26]:


SSH_lattice_TrI.winding_number()


# ### The Topologically Nontrivial Insulator
# The other limit ($|t_{inter}| < |t_{intra}| $) is the interesting one, although the same
# calculations on a lattice with periodic boundary condition reveal nothing interesting.

# In[27]:


t_intra = -0.5
t_inter = -0.65
H_cell = Qobj( np.array( [[ 0, t_intra ],[t_intra,0]] )  )
T_inter_cell =  Qobj( np.array( [[ 0, 0 ],[t_inter,0]] )  ) 
pSSH_lattice_nTrI = Lattice1d(num_cell=100, boundary = "periodic", cell_num_site = 2, cell_site_dof = [1], Hamiltonian_of_cell = H_cell, inter_hop = T_inter_cell )


# In[28]:


pSSH_lattice_nTrI.plot_dispersion()


# In[29]:


pSSH_H_nt = pSSH_lattice_nTrI.Hamiltonian()
nD = pSSH_H_nt.eigenenergies()


# In[30]:


plt.plot(nD,'ro')
plt.ylabel('eigen values')
plt.show()
plt.close()


# ### Topologically Nontrivial Insulator: Hardwall Boundary Condition
# To reveal the topological nontriviality, we have to put the insulator next to a 
# boundary, i.e. use a hardwall boundary condition and the spectrum of eigenvalues of
# the Hamiltonian shows the conspicuous edge states which plotted as a eigenvector
# are dense towards the edge considerably. However, the dispersion looks completely the
# same.

# In[ ]:


t_intra = -0.5
t_inter = -0.65
cell_H = Qobj( np.array( [[ 0, t_intra ],[t_intra,0]] )  )
inter_cell_T =  Qobj( np.array( [[ 0, 0 ],[t_inter,0]] )  ) 
apSSH_lattice_nTrI = Lattice1d(num_cell=100, boundary = "aperiodic", cell_num_site = 2, cell_site_dof = [1], cell_Hamiltonian = cell_H, inter_hop = inter_cell_T )


# In[ ]:


apSSH_lattice_nTrI.winding_number()


# In[ ]:


apSSH_H_nt = apSSH_lattice_nTrI.Hamiltonian()
[n_D,Vx] = apSSH_H_nt.eigenstates()


# In[ ]:


plt.plot(n_D,'ro')
plt.ylabel('eigen values')
plt.show()
plt.close()


# We observe the much celebrated mid-gap edge states, who appear at mid-gaps. The nontriviality of the topological insulator is revealed with their placement next to a boundary. The edge states are the consequences of the bulk-bundary correspondence. Plotting the absolute values of the eigen-functions reveal that they are localized at the edges. We plot the absolute values of the eigenfunctions corresponding to the two eigenvalues found at mid-gaps.

# In[ ]:


xA = [i for i in range(200)]
Es0 = np.abs(Vx[99])
Es1 = np.abs(Vx[100])
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4), sharey=True)
ax1.plot(xA, Es0, label="State0")
ax2.plot(xA, Es1, label="State1")
ax1.legend()
ax2.legend()
plt.show()
fig.suptitle('Mid-gap Edge states')
plt.close()


# It is quite obvious, the edge states are localized towards the edges.
# The eigenstates are not unique, they are only unique within a phase factor. The routine used in QuTiP only yields eigenstates that are one of the possibilities within the phase factor.

# Let us check out the two edge state eigenvalues as well as the three preceding and three succeeding ones.

# In[ ]:


n_D[96:104]


# It is obvious that, n_D[99] and n_D[100] are values 0 within numerical precision. There is an interesting property of the eigenstates of a chiral Hamiltonian with eigenvalue 0.

# \begin{eqnarray}
# H|\psi_n \rangle = 0 \implies H(\hat{P}_{A/B}| \psi_n\rangle) = \frac{1}{2}H(\hat{\bf{1}} \pm \Sigma_z )| \psi_n \rangle = \frac{1}{2}(\hat{\bf{1}} \mp \Sigma_z )H| \psi_n \rangle = 0
# \end{eqnarray}

# Therefore, $\hat{P}_{A/B} |\psi_n\rangle$ is an eigenstate with eigenvalue 0 as well. As a result, the edge-states can be chosen to have support in only one sublattice.

# In[ ]:


chiral_op_nTrI = apSSH_lattice_nTrI.distribute_operator(sigmaz())
dimH_nTrI = chiral_op_nTrI.dims
identity_H_nTrI = Qobj(np.identity(200), dims=dimH_nTrI)
identity_H_nTrI
PA_200 = 0.5*( identity_H_nTrI + chiral_op_nTrI)
PB_200 = 0.5*( identity_H_nTrI - chiral_op_nTrI)

xA = [i for i in range(200)]
Es0 = np.abs(PA_200*Vx[99])
Es1 = np.abs(PB_200*Vx[99])
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4), sharey=True)
ax1.plot(xA, Es0, label="State0")
ax2.plot(xA, Es1, label="State1")
ax1.legend()
ax2.legend()
plt.show()
fig.suptitle('Mid-gap Edge states')
plt.close()


# $\hat{P}_{A} |\psi_n\rangle$ and $\hat{P}_{B} |\psi_n\rangle$ are orthogonal eigenstates with eigenvalue 0, and we can choose the two edge states to be localized at two edges of the lattice. Just to be sure, we also check if $H\hat{P}_{A/B} |\psi_n\rangle $ is indeed 0.

# In[ ]:


is_null = (np.abs(apSSH_H_nt*PA_200*Vx[99]) < 1E-10).all()
print(is_null)


# In[ ]:


qutip.about()


# In[ ]:


qutip.cite()