Derivation of the Berry Curvature and Bloch Magentic Moment in Graphene

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}}$$