Help with translating Hamiltonian into matrix

Question

Eq. 19 in this paper gives the following Hamiltonian:

$\sigma_a, \tau_a, \eta_a$ are respectively the spin, sublattice pseudospin and valley pseudospin respectively.

Normally, I would have chosen a basis convention such as $\sigma_a, \tau_a, \eta_a$ and plugged in the $2\times 2$ identity matrix for terms that don't appear above. Then I would have taken the Kronecker Delta product. However, in this case, the authors explicitly define $(s_z,\tau_z\,\eta_z)=(\pm 1, \pm 1, \pm 1)$ for $($up/down, sublattice A/B, valley K/K'$)$ respectively. So, as a relatively inexperienced person in the field, I feel like there is some ambiguity. I am having trouble translating the above into a matrix that I can play with numerically. Can someone help me understand the convention here?

I tried two interpretations in Python. np.kron takes the Kronecker product. First attempt:

# Choose basis: (tau_xy, sigma_xyz), replace sz, tauz, etaz by s, tau, eta = +/- 1
# Gives 4-band model, with eye = identity and s_xyz as Pauli matrices
eta = 1; # instead of using eta_z as the Pauli z matrix
h1=h*vf*(eta*kx*np.kron(sx,eye)+ky*np.kron(sy,eye)) # e-e hoppings
h2=lSO*eta*np.kron(eye,sz)*tau # SO coupling
h3=m*np.kron(eye,sz)*tau # break TRS (AFM exchange mag.)
H = h1+h2+h3

Second:

# Choose basis: (eta, sigma, tau) using Pauli matrices for all, but no +/- 1 values.
# Gives 8-band model.
h1=h*vf*(kx*np.kron(sz,np.kron(eye,sx))+ky*np.kron(eye,np.kron(eye,sy))) # e-e hoppings
h2=lSO*np.kron(sz,np.kron(sz,sz)) # SO coupling
h3=m*np.kron(eye,np.kron(sz,sz))# break TRS (AFM exchange mag.)
H = h1+h2+h3

Both of these seem to give me very weird Berry curvature fields when plotted. So, I figured I am doing something wrong. Note: I set $v_F = 1, \hbar=1$.

Answer 1

For this type of calculation, I find MATLAB/Octave to be easier, at least for demonstrating what to do. You can add the np. everywhere afterwards if you want to use Python.

Assuming you have the following matrices already defined in your workspace:

$\eta_z$ = eta
$\tau_x$ = taux
$\tau_y$ = tauy
$\tau_z$ = tauz
$m_z$ = tauz
$\sigma_z$ = sigmaz,

and the following scalars:

$\hbar$ = hbar
$v_F$ = v
$\lambda_{\textrm{SO}}$ = lambda,
$k_x$ = kx
$k_y$ = ky,

the matrix would typically be calculated in the following manner:

eta=kron(kron(kron(eta,eye(length(mz)),eye(length(sigmaz)),eye(length(tauz)); 
mz=kron(kron(eye(length(eta)),mz),eye(length(mz)+length(sigmaz)+length(tauz)))));    
sigmaz=kron(kron(eye(length(eta)+length(mz)),sigmaz),eye(length(tauz)));    
tauz=(kron(eye(length(eta)+length(mz)+length(sigmaz)),tauz)));
taux=(kron(eye(length(eta)+length(mz)+length(sigmaz)),taux)));
tauy=(kron(eye(length(eta)+length(mz)+length(sigmaz)),tauy)));

h1 = hbar*v*(eta*kx*taux + ky*tauy);
h2 = lambda*eta*sigmaz*tauz;
h3 = mz*sigmaz*tauz;
H = h1+h2+h3;

The is because typically:

$\tau_x,\tau_z,$ and $\tau_y$ would all lie in the same Hilbert space, and
$\eta_z,m_z$ and $\sigma_z$ would each lie in a different Hilbert spaces to everything else.

There's a lot of left and right parentheses which you'd have to be careful about matching, and since there's 6 operators involved across 4 different Hilbert spaces, this is not such a simple example. If you wanted to give a simpler example with fewer operators and types of operators, or to explicitly give each matrix, I could try to simplify my above code. The length(operator) functions can probably all be replaced by 2, but my code works more generally for arbitrary spin (not just spin-1/2 which has eigenvalues $\pm$1, but also spin-1 which has eigenvalues -1,0,1, etc.).