# Help with definitions in k-space twisted bilayer graphene model

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 # 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!