Derivation of the Berry Curvature and Bloch Magnetic Moment in Graphene

Question

(I found a workable solution, skip to the "Solution" part to see it) I am attempting to derive equations 2 and 6 from Xiao et al. paper "Valley contrasting physics in graphene" (Link to paper). The Hamiltonian for graphene with a staggered sublattice potential (in other words, a potential energy difference between the A and B sublattices) is given by: $$H = \frac{\sqrt{3}}{2}at(q_x\tau_z\sigma_x+q_y\sigma_y)+\frac{\Delta}{2}\sigma_z$$ Where $a=\text{lattice constant}$, $t=\text{hopping energy}\approx2.82eV$, $\sigma_x,\sigma_y,\sigma_z$ are our Pauli matrices and $\Delta$ is the energy difference between the A and B sublattice with the A sublattice having an additional on site energy of $\frac{\Delta}{2}$ and the B sublattice an additional on site energy of $\frac{-\Delta}{2}$. Also, $\tau_z$ means $\pm1$ and is used to select between the two valleys/Dirac points in the Brillouin zone.

From this Hamiltonian we should be able to get everything that we want to know about the Berry curvature and Bloch magnetic moment. These equations are given in the text (but not derived there or in any of their cited texts) as follows: ($m(k) \ \text{is the magnetic moment and,} \ \Omega(q) \ \text{is the Berry curvature}$): $$m(k)=\tau_z\frac{3ea^2{\Delta}t^2}{4\hbar({\Delta^2}+3q^2a^2t^2)}$$ $$\Omega(q)=\tau_z\frac{3a^2{\Delta}t^2}{2({\Delta^2}+3q^2a^2t^2)^{3/2}}$$ It should be noted that both the magnetic moment of the Bloch electron and the Berry curvature are vector quantities and that $q^2$ is the magnitude of the crystal momentum. Also, the Berry curvature equation listed above is for the conduction band. I should also mention at this point that Xiao has a habit of switching between k and q, with q being the crystal momentum measured relative to the valley in graphene.

With this information in hand I will attempt to derive these equations with the help of Mathematica. However, I will fail and this is where I am hoping to get some help from all of you.

Step 1: First we need to determine the Eigenvalues/Dispersion relationship near the valleys (aka Dirac points of graphene). To do this all we need to do is find the eigenvalues of our Hamiltonian, which has the following form when represented as a matrix: $$H=\begin{pmatrix} \frac{\Delta}{2} & \frac{\sqrt{3}}{2}at({\tau_z}q_x-iq_y) \\ \frac{\sqrt{3}}{2}at({\tau_z}q_x+iq_y) & \frac{-\Delta}{2} \\ \end{pmatrix} $$ Using Mathematica (or I suppose you could also do this by hand relatively easily) we find that the eigenvalues of this matrix are: $${\pm}{\frac{1}{2}}\sqrt{\Delta^2+3a^2t^2q_x^2+3a^2t^2q_y^2}={\pm}{\frac{1}{2}}\sqrt{\Delta^2+3a^2t^2q^2}$$ This dispersion relationship (which is only valid near the Dirac points/Valleys) has the expected properties in that it opens up a gap of $\Delta$ between the valence and conduction bands and is also linear, as any Dirac cone should be.

Step 2: Next we need to find the eigenvectors. For this I followed the normal procedure to find the eigenvector of a matrix: (This is the eigenvector of the '+' valley and conduction band so we use $+q_x$ and the positive eigenvalue) $$\begin{bmatrix} 1 \\ \frac{-\Delta+{\sqrt{\Delta^2+3a^2t^2q_x^2+3a^2t^2q_y^2}}}{\sqrt{3}at(q_x-iq_y)} \\ \end{bmatrix}e^{ik\cdot{r}}$$ I have included $e^{ik\cdot{r}}$ since this eigenvector is supposed to be Bloch solution to the Hamiltonian near the valley. This is one of the places where I think I might have made a mistake. Now we need to normalize our eigenvectors since we are dealing with quantum mechanics: $$A^2\begin{bmatrix} 1 & \frac{-\Delta+{\sqrt{\Delta^2+3a^2t^2q_x^2+3a^2t^2q_y^2}}}{\sqrt{3}at(q_x+iq_y)} \\\end{bmatrix}\begin{bmatrix} 1 \\ \frac{-\Delta+{\sqrt{\Delta^2+3a^2t^2q_x^2+3a^2t^2q_y^2}}}{\sqrt{3}at(q_x-iq_y)} \\ \end{bmatrix}=1$$ $$A^2=\frac{1}{1+\frac{\bigr(\Delta-\sqrt{\Delta^2+3a^2t^2(q_x^2+q_y^2)}\bigr)^2}{3a^2t^2(q_x^2+q_y^2)}}$$

Step 3: Now comes the calculation. I did this part in Mathematica because I didn't think it was worth my time to work through all of the algebra by hand. Others may disagree. In any case, the general formula for the Berry curvature is given by: (again, in this case $q$ is a vector quantity in 3 space) $$\Omega_n(q)=\nabla_q\times\left\langle u_n(q) \middle| i{\nabla}_q \middle| u_n(q) \right\rangle$$ This is equation 1.27 from "Berry phase effects on electronic properties" by Xiao et al. (Link to paper) Since our function $u_n(q)$ is just our eigenvector without the added exponential (This is another part where I could be wrong) and it only has components in $q_x$ and $q_y$ the equation for the Berry curvature can be re-written as: $$\Omega_n(q)=i\left\langle \frac{{\partial}u_n(q)}{{\partial}q_x} \middle| \frac{{\partial}u_n(q)}{{\partial}q_y} \right\rangle-\left\langle \frac{{\partial}u_n(q)}{{\partial}q_y} \middle| \frac{{\partial}u_n(q)}{{\partial}q_x} \right\rangle, \ \text{where} \ u_n(q)=\begin{bmatrix} 1 \\ \frac{-\Delta+{\sqrt{\Delta^2+3a^2t^2q_x^2+3a^2t^2q_y^2}}}{\sqrt{3}at(q_x-iq_y)} \\ \end{bmatrix}A^{1/2}$$ Following this operation ought to give me the Berry curvature for +K valley in the conduction band. At this point I handed the calculation over to Mathematica and told it to simplify its result. The final formula that I obtained was: $$\Omega(q)=\frac{-3a^2t^2(q_x^2+q_y^2)\Delta-\Delta^3+\Delta^2\sqrt{\Delta^2+3a^2t^2(q_x^2+q_y^2)}}{(q_x^2+q_y^2)(\Delta^2+3a^2t^2(q_x^2+q_y^2))^{3/2}}$$ This is very close to equation 6 in Xiao's paper but it's off by a a few terms in the numerator. What I am trying to figure out now is if I made an error in my derivation or if Xiao made a simplification step. I suspect that since we are only dealing with small $q$ that the last term in the numerator can be viwed as $\Delta^3$ and since only positive $\Delta$ makes physical sense the second and third terms cancel in the numerator. I will discuss the equation for the magnetic moment after I have heard back from all of you on whether or not my derivation so far is correct. Just so you don't have to scroll all the way back up here is the equation I was hoping for: (In my equation $\tau_z=-$) $$\Omega(q)=\tau_z\frac{3a^2{\Delta}t^2}{2({\Delta^2}+3q^2a^2t^2)^{3/2}}=\tau_z\frac{3a^2{\Delta}t^2}{2({\Delta^2}+3a^2t^2(q_x^2+q_y^2))^{3/2}}$$

Solution: (at least the one I used) So after posting this question I eventually found a way to solve the problem but I never listed the solution here online, which in retrospect was a tad rude of me. In any case, here's how I went about solving the problem. First, recall that the general equation for calculating the Berry curvature is as follows: $$\Omega_{ \mu \nu}^n(\overrightarrow{R})=i \sum_{n' \neq n} \frac{\left\langle n\middle|\frac{{\partial}H}{{\partial}R^{\mu}} \middle|n'\right\rangle \left\langle n'\middle|\frac{{\partial}H}{{\partial}R^{\nu}} \middle|n\right\rangle- \left\langle n\middle|\frac{{\partial}H}{{\partial}R^{\nu}} \middle|n'\right\rangle \left\langle n'\middle|\frac{{\partial}H}{{\partial}R^{\mu}} \middle|n\right\rangle}{(\epsilon_n - \epsilon_n')^2}$$ Remember that the n and n' kets are our normalized eigenvectors of the Hamiltonian with $\mu$ and $\nu$ representing our vector components.These vector components are of course $q_x$ and $q_y$ and so when we take the derivative of the Hamiltonian, we are taking it with respect to those quasi-momentum vectors which are measured relative to the valley/Dirac point position. The denominator is the difference between the conduction and valance band squared. Next, let's define the above variables in a bit more detail: $$H=\frac{\sqrt{3}}{2}at(q_x \tau_x \sigma_x + q_y \sigma_y)+\frac{\delta}{2}\sigma_z$$ $$(\epsilon_n - \epsilon_n')^2 = \delta^2 + 3a^2 t^2 (q_x^2 + q_y^2)$$ $$\left|n\right\rangle = A^{1/2} \left[\frac{-i \left(\Delta+\sqrt{\Delta^2 + 3 a^2 t^2 q_y^2 + 3a^2 t^2 q_x^2 \tau_z^2}\right)}{\sqrt{3} at(q_y-i q_x \tau_z)}\right]$$ $$\left|n'\right\rangle = B^{1/2} \left[\frac{-i \left(\Delta+\sqrt{\Delta^2 + 3 a^2 t^2 q_y^2 + 3a^2 t^2 q_x^2 \tau_z^2}\right)}{\sqrt{3} at(q_y-i q_x \tau_z)}\right]$$ In this case, $A^{1/2}$ and $B^{1/2}$ are the normalization factors for our eigenvectors. Several symbols will not be introduced to simplify this derivation to make it less cumbersome in additional to a a definition of the normalization factors: $$\gamma = \sqrt{\Delta^2 + 3 a^2 t^2 q_y^2 + 3a^2 t^2 q_x^2 \tau_z^2}$$ $$\eta = \sqrt{3}at(q_y + iq_x \tau_z)$$ $$B \left\langle n' \middle|n'\right\rangle=1 \rightarrow B=\left(1+\frac{(\Delta + \gamma)^2}{\eta^* \eta}\right)^{-1}$$ $$A \left\langle n \middle|n\right\rangle=1 \rightarrow B=\left(1+\frac{(\Delta - \gamma)^2}{\eta^* \eta}\right)^{-1}$$ The numerator of the Berry curvature equation can then be shown to have the following form: $$\left\langle n\middle|\frac{{\partial}H}{{\partial}R^{\mu}} \middle|n'\right\rangle \left\langle n'\middle|\frac{{\partial}H}{{\partial}R^{\nu}} \middle|n\right\rangle- \left\langle n\middle|\frac{{\partial}H}{{\partial}R^{\nu}} \middle|n'\right\rangle \left\langle n'\middle|\frac{{\partial}H}{{\partial}R^{\mu}} \middle|n\right\rangle = \frac{2i \Delta \tau_z \sqrt{\Delta^2 + 3a^2 t^2 (q_y^2 + q_x^2 \tau_z^2)}}{q_y^2 + q_x^2 \tau_z^2}(A^{1/2}B^{1/2})^2$$ $$\frac{2i \Delta \tau_z \sqrt{\Delta^2 + 3a^2 t^2 (q_y^2 + q_x^2 \tau_z^2)}}{q_y^2 + q_x^2 \tau_z^2}(A^{1/2}B^{1/2})^2 = \frac{-6i \Delta \tau_z \gamma a^2 t^2}{\eta^* \eta}(A^{1/2}B^{1/2})^2$$ $$(A^{1/2}B^{1/2})^2 = \frac{\gamma^2 - \Delta^2}{4 \gamma^2}$$ Altogether then, our equation now reads as follows: $$\Omega_{q_x q_y}^{n}\left(\overrightarrow{R} \right) = i \frac{-6i \Delta \gamma \tau_z a^2 t^2}{\eta^* \eta}\left(\frac{\gamma^2 - \Delta^2}{4 \gamma^2}\right)\left(\frac{1}{\gamma^2}\right)$$ This simplifies to: (after substituting in our original variables) $$\Omega_{q_x q_y}^{n} = \frac{3 a^2 t^2 \Delta \tau_z}{2 \left(\Delta^2 + 3 a^2 t^2 (q_y^2 + q_x^2)\right)^{3/2}}$$ Which is the same as equation 6 in Xiao's paper on valley contrasting physics in graphene (minus a few changes related to switching to the momentum magnitude). Knowing the Berry curvature near the valleys it is now possible to also calculate the Bloch magnetic moment using an un-numbered equation equation from Xiao's paper: $$m(k)=\left(\frac{e}{\hbar}\right)\epsilon(k)\Omega(k)$$ $$m(k)=\frac{3 e a^2 t^2 \Delta \tau_z}{4 \hbar \left(\Delta^2 + 3 a^2 t^2 (q_y^2 + q_x^2)\right)}$$ This is the same as equation 2 in Xiao's paper. This is as far as I've gotten with this derivation and it's all I needed for the work that I was doing.

Answer 1

In addition to your solution which explicitly derives the eigenstates in details, I would recommend a different method which uses the polar parametrisation of the $2\times 2$ Hamiltonian. This is actually shown in the second paper you cited (Xiao et. al., 2009), in section I.A.3 on two level systems. In particular, a $2\times 2$ Hamiltonian is written as $H=\mathbf{h}(\mathbf{R}).\boldsymbol{\sigma}$, where $\mathbf{h}=h(\sin \theta \cos \phi, \sin \theta \sin \phi, \cos \theta)$ is the spherical parametrisation. After some derivations, the Berry curvature for the general case $\mathbf{h}(\mathbf{R})$ is given by equation (1.19) in the paper \begin{equation} \Omega_{R_1R_2} = \frac{1}{2}\frac{\partial (\phi, \cos \theta)}{\partial (R_1, R_2)} \end{equation} For the case of Cartesian parametrisation $\mathbf{h}=(h_x,h_y,h_z)$, you can convince yourself that \begin{equation} \Omega_{h_x,h_y} = \frac{1}{2}\frac{h_z}{h^3} \end{equation}

In particular, for a hexagonal lattice with staggered potential, $\mathbf{h}=\left(v_F\tau_zq_x,v_Fq_y,\frac{\Delta}{2}\right)$ where $v_F=\frac{\sqrt{3}}{2}ta$. We have $h=\sqrt{h_x^2+h_y^2+h_z^2} = \sqrt{v_F^2 + \frac{\Delta^2}{4}}$ \begin{equation} \Omega(\mathbf{q}) := \Omega_{q_x,q_y} = \tau_z v_F^2 \Omega_{h_x,h_y}=\tau_z v_F^2 \frac{\frac{\Delta}{2}}{\left(v_F^2 + \frac{\Delta^2}{4}\right)^{3/2}} \end{equation} which is equal to the expression for the Berry phase in the first cited paper (Xiao et. al., 2007). If you want to obtain the eigenstates of H, the change of variables from $(h, \theta, \phi)$ to $(q_x,q_y,\Delta)$ can also be made to equation (1.16) of Xiao et. al., 2009.

Answer 2

as a matter of fact I was going through the paper that you have cited above a few days back and I was faced with the same problem. Our task here is to calculate $\Omega(q)=\partial_{\theta}A_{q}-\partial_{q}A_{\theta}$.The reason you are getting the expression for $\Omega(q)$. It can be seen from the expressions for $\Omega(q)$ that in the end, only $A_{\theta}$ will contribute to the Berry curvature. To see this , observe that the expressions for $A_{\theta}$ and $A_{q}$ are independent of theta. Therefore, in the end, you are left with the simple expression for $\Omega(q)$, which is $-\partial_{q}A_{\theta}$.

As to the expression that you are getting, I believe the error stems form the fact that at first sight, one may be misled(as I was) into assuming that $\partial_{\theta}$ should in fact be $\frac{\partial_{\theta}}{r}$ on dimensional grounds. However, it should be noted that $A_{q}$ and $A_{\theta}$ have different dimensions and the expression for the Berry curvature is dimensionally consistent. This is the origin of the additional terms in your expression.. Hope this helps