Calculating integration using for loops

My question is about 2D square Brillouin zone integrations.

Problem Summary:

I am rederiving one of the results of this publication. The mathematical expressions for and are given in two different forms, Eq 10 and Eq 12. Expected results are . However, when I calculate using Eq.10 and Eq.12, I get confilicting results. I have attached my codes both for Eq10.m and Eq12.m.

Equation 10:

The mathematical expression for are

--- Eq. 10a

--- Eq. 10b

here

• 
• 
• 
• 
• with as 2-by-2 matrices (given in code)
• is th eigenvalue of matrix H
• is eigenvalues of H
• are constants
• I use approximation

Equation 12:

The mathematical expression for are

--- Eq. 12a

--- Eq. 12b

here

• 
• 
• are constants

My Results:

So far, for Eq.10, I get (which are expected results). For Eq.12, I get opposite results . I don't know what exactly I am missing here. Please help!

My Codes:

Eq10.m

% ============= Eq10.m =============
clear; clc;
% ===== inputs ====
NBZ = 100; %number of BZ points
eta = 1e-1; %for delta function approximation
EF = 100; %E_F


% ===== constants =====
me = 9.1e-31/(1.6e-19*(1e10)^2);
aR = 2.95; 
muBx = 0.1; 
m = 0.32*me;
hbar = 6.5911e-16;
Gamma = 1.9268e+13;
e = 1.6e-19;
sx = [0 1;1 0];
sy = [0 -1i;1i 0];
sz = [1 0;0 -1];

% ===== k-points =====
kx1 = -pi; kx2 = pi;
ky1 = -pi; ky2 = pi;

kxs = linspace(kx1,kx2,NBZ+1); kxs(end) = [];
kys = linspace(ky1,ky2,NBZ+1); kys(end) = [];
dkx = kxs(2)-kxs(1);
dky = kys(2)-kys(1);

% ===== for loops to calculate \alpha =====
alpha_x = 0;
alpha_y = 0;
for ikx = 1:NBZ %loop over kx
    kx = kxs(ikx);
    for iky = 1:NBZ %loop over ky
        ky = kys(iky);

        % ===== constructing H matrix ==== 
        k = sqrt(kx^2+ky^2);
        
        HR = hbar^2*k^2/(2*m)*eye(2) + aR*(kx*sy - ky*sx);
        Hzee = muBx*sx;
        H = HR + Hzee;
        
        % ===== diagonalizing H =====
        [evecs,evals] = eig(H);
        
        % ===== vx, vy matrix ====
        vx = 1/hbar*[     0, -aR*1i
            aR*1i,      0];
        
        vy = 1/hbar*[   0, -aR
            -aR,   0];

        % ===== Sx, Sy matrix =====
        Sx = hbar * (evecs' * sx * evecs); 
        Sy = hbar * (evecs' * sy * evecs); 
        
        % ===== X, Y, and J matrix =====
        Xx = 0.5*(Sx*vx + vx*Sx);  
        Xy = 0.5*(Sx*vy + vy*Sx); 
        
        Yx = 0.5*(Sy*vx + vx*Sy);
        Yy = 0.5*(Sy*vy + vy*Sy);
        
        Jx = -e*vx;
        Jy = -e*vy;
        
        % ===== loop over n =====
        for n = 1:2
            Un = evecs(:,n); %nth eigenvector
            En = evals(n,n); %nth eigenvalue
            
            delta_fun = 1/(eta*sqrt(2*pi))*exp(-0.5*((En-EF)^2/(eta^2)));
            Xx_k = Un'*Xx*Un; %value of X_x at given k 
            Xy_k = Un'*Xy*Un;
            Yx_k = Un'*Yx*Un;
            Yy_k = Un'*Yy*Un;
            Jx_k = Un'*Jx*Un;
            Jy_k = Un'*Jy*Un;
            
            crossPx = Xx_k*Jy_k - Xy_k*Jx_k; %cross product for \alpha_x
            crossPy = Yx_k*Jy_k - Yy_k*Jx_k; %cross product for \alpha_y
            
            alpha_x = alpha_x + 1/(2*Gamma)*crossPx*delta_fun*dkx*dky;
            alpha_y = alpha_y + 1/(2*Gamma)*crossPy*delta_fun*dkx*dky;
     
        end
        
    end
end


alpha_x = alpha_x/(1.6e-19); % 1.6e-19 factor comes from units (you can ignore it)
alpha_y = alpha_y/(1.6e-19); 


[alpha_x, alpha_y]
%it gives:  [0.0000 + 0.0000i  -5.0997 - 0.0000i]

Eq12.m

% ============= Eq12.m =============
clear; clc;
% ===== inputs ====
NBZ = 100; %number of BZ points
EF = 100; %E_F


% ===== constants =====
me = 9.1e-31/(1.6e-19*(1e10)^2);
aR = 2.95; 
muBx = 0.1; 
m = 0.32*me;
hbar = 6.5911e-16;
Gamma = 1.9268e+13;
e = 1.6e-19;
sx = [0 1;1 0];
sy = [0 -1i;1i 0];
sz = [1 0;0 -1];

% ===== k-points =====
kx1 = -pi; kx2 = pi;
ky1 = -pi; ky2 = pi;

kxs = linspace(kx1,kx2,NBZ+1); kxs(end) = [];
kys = linspace(ky1,ky2,NBZ+1); kys(end) = [];
dkx = kxs(2)-kxs(1);
dky = kys(2)-kys(1);


% ===== for loops to X_x, X_y, Y_x, Y_y, J_x, J_y at each k and n=====
for ikx = 1:NBZ %loop over kx
    kx = kxs(ikx);
    for iky = 1:NBZ %loop over ky
        ky = kys(iky);

        % ===== constructing H matrix ==== 
        k = sqrt(kx^2+ky^2);
        
        HR = hbar^2*k^2/(2*m)*eye(2) + aR*(kx*sy - ky*sx);
        Hzee = muBx*sx;
        H = HR + Hzee;
        
        % ===== diagonalizing H =====
        [evecs,evals] = eig(H);
        
        % ===== vx, vy matrix ====
        vx = 1/hbar*[     0, -aR*1i
            aR*1i,      0];
        
        vy = 1/hbar*[   0, -aR
            -aR,   0];

        % ===== Sx, Sy matrix =====
        Sx = hbar * (evecs' * sx * evecs); 
        Sy = hbar * (evecs' * sy * evecs); 
        
        % ===== X, Y, and J matrix =====
        Xx = 0.5*(Sx*vx + vx*Sx);  
        Xy = 0.5*(Sx*vy + vy*Sx); 
        
        Yx = 0.5*(Sy*vx + vx*Sy);
        Yy = 0.5*(Sy*vy + vy*Sy);
        
        Jx = -e*vx;
        Jy = -e*vy;
        
        % ===== loop over n =====
        for n = 1:2
            Un = evecs(:,n); %nth eigenvector
            En = evals(n,n); %nth eigenvalue
            
            Xx_k(ikx,iky,n) = Un'*Xx*Un; %value of X_x at given k 
            Xy_k(ikx,iky,n) = Un'*Xy*Un;
            Yx_k(ikx,iky,n) = Un'*Yx*Un;
            Yy_k(ikx,iky,n) = Un'*Yy*Un;
            Jx_k(ikx,iky,n) = Un'*Jx*Un;
            Jy_k(ikx,iky,n) = Un'*Jy*Un;
            E_k(ikx,iky,n) = En;
        end
        
    end
end
[X,Y] = meshgrid(kxs,kys);

% === finding curl for each n
for n=1:2
    [curlX(:,:,n),~]=curl(X,Y,Xx_k(:,:,n),Xy_k(:,:,n)); 
    [curlY(:,:,n),~]=curl(X,Y,Yx_k(:,:,n),Yy_k(:,:,n));
end
%% now evaluating integration over omega_x and omega_y

omega_x = 0;
omega_y = 0;
for n = 1:2 
    
    SumBZomega = zeros(1,2);
    for ikx = 1:NBZ %over kx
        kx = kxs(ikx);
        
        for iky = 1:NBZ %over kx
            ky = kys(iky);
           
            En = E_k(ikx,iky,n); %En at each kx,ky, and n
            
            if En<=EF %applying En<=EF condition
                SumBZomega(1) = SumBZomega(1) + curlX(ikx,iky,n)*dkx*dky;
                SumBZomega(2) = SumBZomega(2) + curlY(ikx,iky,n)*dkx*dky;
            end
            
        end
    end
    omega_x = omega_x + SumBZomega(1); 
    omega_y = omega_y + SumBZomega(2);
end

% ==== finally alpha_x, alpha_y ====
alpha_x = e/(2*Gamma*hbar*4*pi^2)*omega_x/(1.6e-19); %1.6e-19 coming from units
alpha_y = e/(2*Gamma*hbar*4*pi^2)*omega_y/(1.6e-19); 


[alpha_x, alpha_y]
%it gives [0.4840 - 0.0000i  -0.0000 + 0.0000i]