Numerically calculating the berry curvature for graphene

Question

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

In order to do this I am attempting to use this equation for the berry curvature . 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

,

and

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_] := 
  3*sign*
   Sqrt[3 + 2*Cos[Sqrt[3]*kx/2 + ky/2] + 2*Cos[Sqrt[3]*kx/2 - ky/2] + 
     2*Cos[ky]];

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[
   N[Sum[
     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.

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)$ .