i'm trying to reproduce this density plot for the Berry curvature in the Brillouin zone of graphene from this website.

Berry Curvature density

In order to do this I am attempting to use this equation for the berry curvature Berry Curvature equation. In which $\vec{\mathcal{R}}$ represents a set of parameters the Hamiltonian $\mathcal{H}$ depends on, $|{n(\vec{\mathcal{R}})>}$ is an eigenstate of the Hamiltonian with energy $\epsilon_n$ and $\omega^n_{\mu \nu} (\vec{\mathcal{R}})$ is the Berry curvature.

For now, I'm using the non-Haldane Hamiltonian for graphene, with nearest-neighbour hopping only to test my calculation. This is given by

Graphene hamiltonian,


Energy Eigenstate

are $\epsilon$ and $|n>$ and the parameters $\vec{\mathcal{R}}=(k_x,k_y)$. Where i have taken $|n>$ and $|n'>$ to be the states corresponding to the conduction and valence bands of graphene. I think I must be interpreting / using this equation incorrectly as i can't seem to get any reasonable results out when i try to do this calculation in mathematica. Here is the code i have tried to use to calculate the curvature over the Brillouine zone so far (note:vx and vy are simply the derivatives of E w.r.t. to kx and ky):

kx1 = Range[0, 2*Pi, 2*Pi/99];
kx2 = Range[0, 2*Pi, 2*Pi/99];

f[kx_, ky_] := 
  Exp[I*ky / Sqrt[3]] + Exp[-I*kx/2 - I*ky/(2*Sqrt[3])] + 
   Exp[I*kx/2 - I*ky/(2*Sqrt[3])];

Energy[sign_, kx_, ky_] := 
   Sqrt[3 + 2*Cos[Sqrt[3]*kx/2 + ky/2] + 2*Cos[Sqrt[3]*kx/2 - ky/2] + 

u[n_, kx_, ky_] := {1, -n*f[kx, ky]/Abs[f[kx, ky]]};

vx[kx_, ky_] := 
  3*(-Sqrt[3]*Sin[Sqrt[3]*kx/2 + ky/2] - 
      Sqrt[3]*Cos[Sqrt[3]*kx/2 - ky/2]) / Energy[+1, kx, ky];

vy[kx_, ky_] := 
  3*(-Sin[Sqrt[3]*kx/2 + ky/2] + Sin[Sqrt[3]*kx/2 - ky/2] - 
      2*Sin[ky]) / Energy[+1. kx, ky];

vxnm[n_, m_, kx_, ky_] := 
  Conjugate[u[n, kx, ky]] . vx[kx, ky] . u[m, kx, ky];

vynm[n_, m_, kx_, ky_] := 
  Conjugate[u[n, kx, ky]] . vy[kx, ky] . u[m, kx, ky];

berry[n_, m_, kx_, ky_] := 
 vxnm [n, m, kx, ky]* 
  vynm [n, m, kx, ky]/ (Energy[n, kx, ky] - Energy[m, kx, ky])^2

del = Table[
     If[p != q, berry[p, q, kx[[i]], ky[[j]]], 0], {p, -1, 1}, {q, -1,
       1}]], {i, 100}, {j, 100}];

Could someone help me see where i've went wrong / misunderstood something? I've spent a few days on this and am genuinely stumped. Sorry if this question is very basic.

  • $\begingroup$ Did you find a way to do it? Also did you first try implementing problems whose analytical solutions are available? $\endgroup$
    – johnny123
    Jul 17, 2023 at 12:21

1 Answer 1


If I'm understanding your explanation and code correctly, you are calculating the interband velocity elements vxnm and vynm with the derivative of the energies. This is wrong. You should think of $ \partial H / \partial k_i $ as an operator, a 2x2 matrix in band space for your case. The derivative of the energies corresponds to the diagonal component $<n(k) \vert \partial H(k) / \partial k_i \vert n(k)>$.

[this is not completely trivial - to show this: $\partial \epsilon_n(k) / \partial k_i = \partial [ <n(k) \vert H(k) \vert n(k)>] / \partial k_i = <n(k) \vert \partial H(k) / \partial k_i \vert n(k)>$, where we use the fact that n(k) are eigenstates of H and integration by parts to show that the terms $<\partial_{k_i} n(k) \vert H(k) \vert n(k)> + <n(k) \vert H(k) \vert \partial_{k_i} n(k)>$ cancel]

The proper vx_nm and vy_nm elements should be calculated from the interband matrix elements $<n(k) \vert \partial_{k_i} H(k) \vert m(k)>$. From your eqs. (9), (10) you should be able to calculate the matrix $\partial_{k_i} H(k)$ .


Your Answer

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