I am trying to finalize my band structure plot for twisted bilayer graphene. I have been having some problems with the plot itself.

To troubleshoot my code, I first just looked at the diagonal part of my Hamiltonian matrix which consists of block diagonal single-layer graphene Dirac Hamiltonians.

The symmetry path I chose for the Moiré BZ in its zigzag direction was $\Gamma \rightarrow K' \rightarrow K \rightarrow \Gamma \rightarrow M $. I choose around 100 wave-vector values between each symmetry point, however when I plot the bands, I get the correct dispersion along the $K'$ to $K$ path, while the others just look weird. See below...

I think this has to do with the actual scaling of the symmetry path. I literally took the values of the corners in momentum space, and that is how I constructed the path, but I am not sure if there's another convention on plotting it. Have you had this same issue, or know of a specific way to plot the bands against the symmetry paths?

Sorry if this question is too simple...

Edit: The Hamiltonian I am using for this system is one with intralayer Dirac terms and the interlayer terms "turned off". I am following the Bistritzer and Macdonald model.

$H = \hbar v_F \sigma^{\theta, \mu}_\zeta \cdot (q - \zeta K^\mu_\zeta + G_i)$

Here the $\sigma$ is the rotated Pauli operator: $e^{-i\theta/2}\sigma e^{i\theta/2}$ and the specific rotated Dirac points are the following. $K$ are the rotated Dirac points and hence the corners of the Moire BZ and the $G_i$ can be a sum of the Moire reciprocal lattice vectors.

For example: The new Moire BZ corners which I will call $K^1$ and $K^2$ are

$K^1 = \frac{4 \pi }{3 \sqrt{3} a}(\cos(\theta/2), \sin(\theta/2))$ $K^2 = \frac{4 \pi }{3 \sqrt{3} a}(\cos(\theta/2), -\sin(\theta/2))$

When I was choosing my path I used the $K^2$ vector as reference. Using what @Franz had suggested I went from $\Gamma-K^2-M-\Gamma$ and I coded the path in the following way:

KX = list(np.linspace(0, K_2[0]*np.sqrt(3)/2, num))+ list(np.linspace(K_2[0]*np.sqrt(3)/2,K_2[0]*np.sqrt(3)/2, 50)) + list(np.linspace(K_2[0]*np.sqrt(3)/2,0,86))
KY = list(np.linspace(0, -K_2[1]/2, num)) + list(np.linspace(-K_2[1]/2,K_1[1]/2, 50)) + list(np.linspace(K_1[1]/2,0, 86))

From here I feed these KX and KY to the Hamiltonian function I wrote, and then will obtain a (len(KX), 4*n) array for the respective eigenvalues.

  • $\begingroup$ Try calculating the eigenvalues for a single point: $\Gamma$. This should remove one of the terms in the parentheses from your $H$, since $K=0$. Do you get a sensible result then? Or do you still get a band gap $>10000$? $\endgroup$
    – Ruslan
    Feb 17, 2021 at 20:45
  • $\begingroup$ BTW, what are the units in the vertical axis of your plot? $\endgroup$
    – Ruslan
    Feb 17, 2021 at 20:48
  • $\begingroup$ @Ruslan I tried just now and the eigenvalues at the gamma point range from -12,000 to 12,000. The lowest eigenvalue is of 6000 so I get a band gap of around 12,000. $\endgroup$
    – MadLad
    Feb 17, 2021 at 20:50
  • $\begingroup$ @Ruslan the units are meV $\endgroup$
    – MadLad
    Feb 17, 2021 at 20:57
  • $\begingroup$ Do you get the same set of eigenvalues in each of the $K$ points of the BZ? $\endgroup$
    – Ruslan
    Feb 17, 2021 at 21:15

2 Answers 2


I can give you some workarounds in order to plot your electron dispersion along high symmetry paths.

The most precise thing to do would be to calculate the absolute value of your k-vector and plot the eigenenergies against it. I had some issues with that, which is why I can give you an even simpler solution to your problem.

Let us consider your hexagonal brillouin zone in 2D. Let us furthermore think of the high symmetry path $\Gamma$-$\mathrm{K}$-$\mathrm{M}$-$\Gamma$ which I drew into the hexagon below for better imagination.

enter image description here

From basic geometry we know that the distance between the center of the hexagon ($\Gamma$) and the corner ($K$) is as long as the sides of the hexagon. That basically gave us all information about the relative lengths between our high symmetry points.

The only thing that is left for you to do is to choose the number of points between two high symmetry points in such a way that they fit the corresponding relative lengths of the high symmetry lines. With that you can plot your computed eigenvalues against the k-index, so basically the natural number of the point on the high symmetry line.

Say between $\Gamma$-$\mathrm{K}$ you choose $100$ points, than $\mathrm{K}$-$\mathrm{M}$ would be half of it, so $50$ points. $\mathrm{M}$-$\Gamma$ is $\frac{\sqrt{3}}{2}$ of the length of $\Gamma$-$\mathrm{K}$, so you use $\frac{\sqrt{3}}{2} \cdot 100$ points for that symmetry line.

I would also recommend you to use matplotlib.pyplot.scatter for plotting instead of the ordenary plt.plot(). Sometimes plt.plot() does not work properly for bandstructures, as you find many energies at one k-point. Another advantage of matplotlib.pyplot.scatter is that you can plot other values, such as the spin's z-expectation value, as a color on top of the energies. I therefore attached the documentation of this matplotlib function for you.


  • $\begingroup$ Hi franz. Thank you for your answer. Just to be clear: I would need to obtain the magnitude of the k-vector between the $\Gamma-K$ point, and then scale the other magnitudes by this wave vector, right? Then if I choose 100 points for this path, then the other paths will rescale corresponding to their ratios. I will try this now. $\endgroup$
    – MadLad
    Feb 17, 2021 at 15:58
  • 1
    $\begingroup$ Correct. Keep in mind though that every point that you use in your calculations increases the total number of points by one. So that the next k-point after you arrived after $100$ steps at $\mathrm{K}$ is labeled $101$. So you basically add up the distance of the path you walked through the Brillouin zone. $\endgroup$
    – franz
    Feb 17, 2021 at 17:12
  • $\begingroup$ So I did this, and I still don't quite get what I want. Could you concretely explain how to construct the path? $\endgroup$
    – MadLad
    Feb 17, 2021 at 20:20
  • $\begingroup$ Depends. Do you use reduced momentum variables that are proportional to $k_i \cdot a_i$ or wavevectors $k \sim 10^{10} \, \mathrm{m}^{-1}$? $\endgroup$
    – franz
    Feb 17, 2021 at 23:27
  • 1
    $\begingroup$ I don't really see that in your $k_x$ & $k_y$ lists, but that's basically the same (reduced momentum variables vs wave vectors). The numbers of points seem correct to me, according to the relative lenghts of the k-paths. I don't understand why your $k_y$-value for the first route is negative but this might be due to your definition of $K^2$ or the layout of your $k_x$-$k_y$ coordinate system. Another thing you should check is whether the $k_y$-value remains constant along the $\mathrm{K}$-$\mathrm{M}$ path, as well as the $k_x$-value is constant for $\mathrm{M}$-$\Gamma$. $\endgroup$
    – franz
    Feb 19, 2021 at 13:15

This doesn't look like a problem with plotting. Look closer (ideally, zoom it in) at this part of your plot:

Notice how in the direction $K-\Gamma$ you have two periods of some oscillatory behavior of the yellow and gray curves, and then they simply diverge away from the center. This looks more like you have set up your Hamiltonian, passing $\vec k$ outside of the first Brillouin zone, and, since the Hamiltonian is a finite-dimensional matrix, rather than the full infinite-dimensional operator, you've pushed it beyond its limits of applicability.

The two oscillation periods also suggest that you may be going too fast across the Brillouin zone. In fact, this might be the sole problem with your code: if you stop at the half of the first period, it might appear to already be the $\Gamma$ point:

  • $\begingroup$ Hi Ruslan! Thanks for your answer. Yes I had absolutely noticed these oscillations and then the divergence, and thought that something was off with my Hamiltonian. However, I thought that all of the high symmetry paths would fall inside of the Brillouin zone, so I am not sure how I would fix my matrix to correctly include these k-points. $\endgroup$
    – MadLad
    Feb 17, 2021 at 16:00
  • $\begingroup$ @MadLad you haven't given any details on your Hamiltonian, on the exact values of $\vec k$ and how you incorporate them into the Hamiltonian, so I don't know where you should look to fix this. Please update your question with the relevant information. $\endgroup$
    – Ruslan
    Feb 17, 2021 at 16:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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