Derivation of the Berry Curvature and Bloch Magentic Moment in Graphene

(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_xtau_zsigma_x+q_ysigma_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_zfrac{3ea^2{Delta}t^2}{4hbar({Delta^2}+3q^2a^2t^2)}$$
$$Omega(q)=tau_zfrac{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^{ikcdot{r}}$$
I have included $e^{ikcdot{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^2begin{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_qtimesleftlangle
u_n(q)
middle|
i{nabla}_q
middle|
u_n(q)
rightrangle$$
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)=ileftlangle
frac{{partial}u_n(q)}{{partial}q_x}
middle|
frac{{partial}u_n(q)}{{partial}q_y}
rightrangle-leftlangle
frac{{partial}u_n(q)}{{partial}q_y}
middle|
frac{{partial}u_n(q)}{{partial}q_x}
rightrangle, 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^2sqrt{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_zfrac{3a^2{Delta}t^2}{2({Delta^2}+3q^2a^2t^2)^{3/2}}=tau_zfrac{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{leftlangle
nmiddle|frac{{partial}H}{{partial}R^{mu}}
middle|n’rightrangle
leftlangle
n’middle|frac{{partial}H}{{partial}R^{nu}}
middle|nrightrangle-
leftlangle
nmiddle|frac{{partial}H}{{partial}R^{nu}}
middle|n’rightrangle
leftlangle
n’middle|frac{{partial}H}{{partial}R^{mu}}
middle|nrightrangle}{(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|nrightrangle = 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’rightrangle = 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 leftlangle n’ middle|n’rightrangle=1 rightarrow B=left(1+frac{(Delta + gamma)^2}{eta^* eta}right)^{-1}$$
$$A leftlangle n middle|nrightrangle=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:
$$leftlangle
nmiddle|frac{{partial}H}{{partial}R^{mu}}
middle|n’rightrangle
leftlangle
n’middle|frac{{partial}H}{{partial}R^{nu}}
middle|nrightrangle-
leftlangle
nmiddle|frac{{partial}H}{{partial}R^{nu}}
middle|n’rightrangle
leftlangle
n’middle|frac{{partial}H}{{partial}R^{mu}}
middle|nrightrangle = 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.