diff --git a/.DS_Store b/.DS_Store
new file mode 100644
index 0000000000000000000000000000000000000000..5008ddfcf53c02e82d7eee2e57c38e5672ef89f6
Binary files /dev/null and b/.DS_Store differ
diff --git a/__pycache__/tSnS.cpython-38.pyc b/__pycache__/tSnS.cpython-38.pyc
index 45462ff99d850960856b55c316c47e6591ea6ea7..18e5191016c906076dfecd7be4d8892488f30e6a 100644
Binary files a/__pycache__/tSnS.cpython-38.pyc and b/__pycache__/tSnS.cpython-38.pyc differ
diff --git a/spin_texture.py b/spin_texture.py
new file mode 100644
index 0000000000000000000000000000000000000000..a9f56354730c342089e930cdfc99902878775987
--- /dev/null
+++ b/spin_texture.py
@@ -0,0 +1,59 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+from momentum import Momentum
+from tSnS import tSnS
+
+
+def spin_texture(Nk=30, bandset=2):
+
+ """
+ Compute Spin texture in BZ
+ """
+
+ #Pauli matrices and orbital identity
+ sigma_x_tau_0 = np.array([[0,0,1,0],[0,0,0,1],[1,0,0,0],[0,1,0,0]])
+ sigma_y_tau_0 = 1j*np.array([[0,0,-1,0],[0,0,0,-1],[1,0,0,0],[0,1,0,0]])
+ sigma_z_tau_0 = np.array([[1,0,0,0],[0,1,0,0],[0,0,-1,0],[0,0,0,-1]])
+
+
+ #Create Mesh and compute spin polarization
+ mesh = BZ.Bravais_mesh_simple(Nk)
+ ekf = np.zeros((Nk**2, model.N))
+ ukf = np.zeros((Nk**2, model.N, model.N), dtype=complex)
+ spin_pol = np.zeros((Nk**2, 3, model.N), dtype=complex)
+
+ for i in range(len(mesh)):
+ ekf, ukf = np.linalg.eigh(model.BH_Hamiltonian(mesh[i], bandset=bandset))
+ spin_pol[i,0,:] = (np.diagonal(ukf.T.conjugate()@sigma_x_tau_0@ukf)).real
+ spin_pol[i,1,:] = (np.diagonal(ukf.T.conjugate()@sigma_y_tau_0@ukf)).real
+ spin_pol[i,2,:] = (np.diagonal(ukf.T.conjugate()@sigma_z_tau_0@ukf)).real
+
+
+ fig, axes = plt.subplots(3, model.N, sharex=True, sharey=True, figsize=(4,3))
+ for i in range(3):
+ for j in range(model.N):
+ axes[i,j].scatter(mesh[:,0], mesh[:,1], c=spin_pol[:,i,j], cmap=plt.cm.RdBu, vmin=-1, vmax=1, s=1)
+ axes[i,j].axis('scaled')
+ axes[i,j].set_xticks([])
+ axes[i,j].set_yticks([])
+
+ axes[0,0].set_ylabel(r'$\langle S_x \rangle$')
+ axes[1,0].set_ylabel(r'$\langle S_y \rangle$')
+ axes[2,0].set_ylabel(r'$\langle S_z \rangle$')
+
+ axes[0,0].set_xlabel(r'$n_b = 1$', labelpad=-70)
+ axes[0,1].set_xlabel(r'$n_b = 2$', labelpad=-70)
+ axes[0,2].set_xlabel(r'$n_b = 3$', labelpad=-70)
+ axes[0,3].set_xlabel(r'$n_b = 4$', labelpad=-70)
+
+
+ plt.subplots_adjust(hspace=0, wspace=0)
+ plt.show()
+
+
+if __name__ == '__main__':
+
+ model = tSnS()
+ BZ = Momentum(model.G1, model.G2)
+ spin_texture()
diff --git a/tSnS.py b/tSnS.py
index 7bfe6bee2a70bedac23f4f0e829702eee6fe48b9..9f1ed2fc46da2862e84ba55b12f9c28498b0936c 100644
--- a/tSnS.py
+++ b/tSnS.py
@@ -61,25 +61,15 @@ class tSnS:
if bandset == 1:
- popt = [-0.28147216773927675, 0.003966634424731179,
- 0.0029073602294481787, -4.4983518626626615e-07,
- -0.004052637332926278, 0.00016072142349528914,
- -3.7721284264216035e-05, -0.0009141524777857822,
- -0.0011772772860306072, -0.0014951352789262786,
- 0.0009438170971867384, -9.088937950755711e-05,
- -0.000692733427968305, 0.0002266622565665427,
- 0.0005191476715863555]
+ popt = [-0.28159930344528794, 0.003968154118035716, 0.0030637767514870547, -0.00010069737918143397,
+ -0.004016436852695664, -0.0008890332864148689, -0.001246870994836976, 0.000848707931621465,
+ -0.0012602009377081112, -0.0007056692333075815, -0.0005572460653001037, -0.00045337702969447147, np.pi/30]
else:
- popt = [-0.32047590155917033, 0.0021542300052420408,
- 0.001469771721045236, 9.574936436961549e-05,
- 0.005574114270507838, -0.00026463100214217314,
- 0.0004064637000665132, 0.0004858014275756446,
- 0.0006724769177247459, -0.00029328415914492434,
- -0.0002932514194152098, 1.0106788563322311e-05,
- 0.0007016914741439542, 0.0021456042488447995,
- 0.0018535546929481545]
+ popt = [-0.32044985418524075, 0.002055488150036208, 0.0018580964148542333, 9.451241186629052e-05,
+ 0.0055709631952728005, 0.0005463154794551496, 0.0006627259014835454, -0.000292263072270747,
+ -0.00029205370932945825, 0.0006839685785136236, -0.0021245180394823893, -0.0017931410640271215, np.pi/8]
@@ -87,11 +77,11 @@ class tSnS:
t_AA_0 = popt[0]
t_AB_0, t_AB_1 = popt[1], popt[2]
t_AA_x, t_AA_y = popt[3], popt[4]
- t_AB_mx, t_AB_px = popt[5], popt[6]
- t_AB_my, t_AB_py = popt[7], popt[8]
- t_AA_m, t_AA_p = popt[9], popt[10]
- t_AA_2x, t_AA_2y = popt[11], popt[12]
- tR_AB_0, tR_AB_1 = popt[13], popt[14]
+ t_AB_my, t_AB_py = popt[5], popt[6]
+ t_AA_m, t_AA_p = popt[7], popt[8]
+ t_AA_2y = popt[9]
+ tR_AB_0, tR_AB_1 = popt[10], popt[11]
+ phi_AB = popt[12]
##################################################
@@ -102,15 +92,6 @@ class tSnS:
ham[:,0,:,0] += t_AA_0*pauli[0]
ham[:,1,:,1] += t_AA_0*pauli[0]
- #Nearest-Neighbor AB/BA coupling (1. nn)
- hop = t_AB_0*(1.+self.phi_lm(k,0,-1))*pauli[0]
- ham[:,0,:,1] += hop
- ham[:,1,:,0] += np.conjugate(hop).T
-
- hop = t_AB_1*(self.phi_lm(k,1,0)+self.phi_lm(k,1,-1))*pauli[0]
- ham[:,0,:,1] += hop
- ham[:,1,:,0] += np.conjugate(hop).T
-
#Nearest-Neighbor AA/BB coupling (2. nn)
hop = t_AA_x*(self.phi_lm(k,1,0))*pauli[0]
ham[:,0,:,0] += hop
@@ -150,12 +131,6 @@ class tSnS:
ham[:,1,:,1] += np.conjugate(hop)*pauli[0]
#Second shell (5. nn)
- hop = t_AA_2x*self.phi_lm(k,2,0)*pauli[0]
- ham[:,0,:,0] += hop
- ham[:,0,:,0] += np.conjugate(hop).T
- ham[:,1,:,1] += hop
- ham[:,1,:,1] += np.conjugate(hop).T
-
hop = t_AA_2y*self.phi_lm(k,0,2)*pauli[0]
ham[:,0,:,0] += hop
ham[:,0,:,0] += np.conjugate(hop).T
@@ -166,18 +141,10 @@ class tSnS:
#AB coupling
################
#Edge states of second shell (4. nn)
- hop = t_AB_mx*(self.phi_lm(k,2,0)+self.phi_lm(k,2,-1))*pauli[0]
- ham[:,0,:,1] += hop
- ham[:,1,:,0] += np.conjugate(hop).T
-
hop = t_AB_my*(self.phi_lm(k,1,1)+self.phi_lm(k,1,-2))*pauli[0]
ham[:,0,:,1] += hop
ham[:,1,:,0] += np.conjugate(hop).T
- hop = t_AB_px*(self.phi_lm(k,-1,0)+self.phi_lm(k,-1,-1))*pauli[0]
- ham[:,0,:,1] += hop
- ham[:,1,:,0] += np.conjugate(hop).T
-
hop = t_AB_py*(self.phi_lm(k,0,1)+self.phi_lm(k,0,-2))*pauli[0]
ham[:,0,:,1] += hop
ham[:,1,:,0] += np.conjugate(hop).T
@@ -203,6 +170,23 @@ class tSnS:
ham[:,0,:,1] += hop
ham[:,1,:,0] += np.conjugate(hop).T
+
+ ######################################
+ #SOC Kane-Mele Coupling
+ ######################################
+ #Nearest-neighbors (1. nn)
+ km = np.array([[np.exp(1j*phi_AB),0],[0,np.exp(1j*phi_AB).conjugate()]])
+ hop = t_AB_0*(1.+self.phi_lm(k,0,-1))*km
+ ham[:,0,:,1] += hop
+ ham[:,1,:,0] += np.conjugate(hop).T
+
+ km = np.array([[np.exp(-1j*phi_AB),0],[0,np.exp(-1j*phi_AB).conjugate()]])
+ hop = t_AB_1*(self.phi_lm(k,1,0)+self.phi_lm(k,1,-1))*km
+ ham[:,0,:,1] += hop
+ ham[:,1,:,0] += np.conjugate(hop).T
+
ham = ham.reshape((4,4))
+ assert np.allclose(ham, ham.conjugate().T, atol=1e-12)
+
return ham