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: enter image description here

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],

Any advice? I have indicated my confusion-points in bold. Thanks in advance!


1 Answer 1


I don't think I can use this Hamiltonian to recover the band structure. The band structure should use a tight-binding Hamiltonian with several unit cells. A colleague told me that this Hamiltonian is just for their discussion in Figure 4.

  • 1
    $\begingroup$ Nice to see you figured this one out, and thanks for making the answer a community wiki rather than writing an ordinary self-answer for this! $\endgroup$ Jul 16, 2021 at 2:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.