I am trying to numerically do calculations using Eq. 8 of MacDonald's simple model for twisted Bilayer graphene. I only want to calculate the Berry phase. However, I don't think I have my definitions right (as my band diagram looks very off). Would someone mind clarifying where I went wrong with my definitions and code below? Thanks.
$$H_k=\begin{bmatrix} h_k(\theta/2) & T_b & T_{tr} & T_{tl}\\ T_b^\dagger & h_{k_{b}}(-\theta/2) & 0 & 0 \\ T_{tr}^\dagger & 0 & h_{k_{tr}}(-\theta/2) & 0 \\ T_{tl}^\dagger & 0 & 0 & h_{k_{tl}}(-\theta/2) \end{bmatrix},$$
where $k_j = k + q_j$. Additionally: $$ T_1=\begin{bmatrix} 1 & 1 \\ 1 & 1 \end{bmatrix}, \text{ } T_2=e^{-i G^{(2)'}\cdot d}\begin{bmatrix} e^{-i \phi} & 1 \\ e^{i \phi} & e^{-i \phi} \end{bmatrix}, \text{ and } T_3=e^{-i G^{(3)'}\cdot d}\begin{bmatrix} e^{i \phi} & 1 \\ e^{-i \phi} & e^{i \phi} \end{bmatrix} \text{ (for } \phi = 2\pi/3\text{)} $$
and $$ h_k(\theta)=-v k \begin{bmatrix} 0 & e^{i(\theta_k - \theta)} \\ e^{-i(\theta_k - \theta)} & e^{-i \phi} \end{bmatrix}. $$
In Figure 1, they give the hopping directions of the other Dirac points as: \begin{align} j &= 1 (\equiv b?) \rightarrow (0,-1) \equiv (q_{x_1},q_{y_1}) \\ j &= 2 (\equiv tr?) \rightarrow (\sqrt{3}/2,1/2) \equiv (q_{x_2},q_{y_2}) \\ j &= 3 (\equiv tl?) \rightarrow (-\sqrt{3}/2,1/2) \equiv (q_{x_3},q_{y_3}). \end{align}
Below eq. 7, they state that the spectrum is independent of the vector $d$ for $\theta \neq 0$, and so I set $d = (0,0)$. Is this valid?
Therefore, I set the exponential factor involving $G\cdot d$ in $T_1$ and $T_2$ to $=1$. With this interpretation, $T_1\equiv T_b,T_2\equiv T_{tr}$ and $T_3\equiv T_{tl}$ seem pretty clear. However, is this correct?
Next, I set $v=1$ (the Dirac velocity). However, I am not sure whether this is justifiable because the authors involve some 'renormalized velocity' in eq. 11, involving a parameter $\alpha$. If I don't care about the analysis in Fig. 4, am I okay with setting $v=1$?
Finally, my biggest confusion is with the definition of $h_i(\theta)$. Since I have my parameter space as $(k_x,k_y)$, am I right in setting $k=\sqrt{k_x^2+k_y^2}$? Since I am taking $\theta$ = twist angle, I leave that as a constant - correct? But, how do I define $\theta_k$? It is defined as the 'momentum orientation relative to the x axis'. My guess was that: \begin{align} h_k(\theta/2) &\leftarrow \arctan2(ky+q_{y_0},kx+q_{x_0}) - \theta/2 \\ h_{k_b}(-\theta/2) &\leftarrow \arctan2(ky+q_{y_1},kx+q_{x_1}) + \theta/2 \\ h_{k_{t_r}}(-\theta/2) &\leftarrow \arctan2(ky+q_{y_2},kx+q_{x_2}) + \theta/2 \\ h_{k_{t_l}}(-\theta/2) &\leftarrow \arctan2(ky+q_{y_3},kx+q_{x_3}) + \theta/2, \end{align}
where above, I indicated only the argument fed into the function; and used $(q_{x_0},q_{y_0}) = (0,0)$. I assumed this because the authors did not seem to mention this. Also, I assumed the 'orientation' meant 'angle', and that it should be calculated with respect to the directional vectors above the horizontal line in this post.
For $-1\leq k_x \leq 1$ and $k_y=0$, I get the following 'band diagram' that is clearly incorrect:
My Python code for defining the Hamiltonian is as follows:
t = params[0] # twist angle
v = 1 # dirac velocity
k = np.sqrt(kx**2+ky**2)
q1x = 0; q1y = 0
q2x = 0; q2y = -1
q3x = np.sqrt(3)/2; q3y = 1/2
q4x = -np.sqrt(3)/2; q4y = 1/2
tk1 = np.arctan2(ky+q1y,kx+q1x)-t/2
tk2 = np.arctan2(ky+q2y,kx+q2x)+t/2
tk3 = np.arctan2(ky+q3y,kx+q3x)+t/2
tk4 = np.arctan2(ky+q4y,kx+q4x)+t/2
hk = -v*k*np.matrix([[0,np.exp(1j*tk1)],[np.exp(-1j*tk1),0]])
hkb = -v*k*np.matrix([[0,np.exp(1j*tk2)],[np.exp(-1j*tk2),0]])
hktr = -v*k*np.matrix([[0,np.exp(1j*tk3)],[np.exp(-1j*tk3),0]])
hktl = -v*k*np.matrix([[0,np.exp(1j*tk4)],[np.exp(-1j*tk4),0]])
phi = 2*np.pi/3
G2 = 0 # below eq 7 say G.d = 0
G3 = 0 # below eq 7 say G.d = 0
Tb = np.matrix([[1,1],[1,1]])
Ttr = np.exp(-1j*G2)*np.matrix([[np.exp(-1j*phi),1],[np.exp(1j*phi),np.exp(-1j*phi)]])
Ttl = np.exp(-1j*G3)*np.matrix([[np.exp(1j*phi),1],[np.exp(-1j*phi),np.exp(1j*phi)]])
zer = np.zeros((2,2))
H = np.block([[hk,Tb,Ttr,Ttl],
[Tb.H,hkb,zer,zer],
[Ttr.H,zer,hktr,zer],
[Ttl.H,zer,zer,hktl]])
Any advice? I have indicated my confusion-points in bold. Thanks in advance!
