I have a triangular lattice model. In $k$-space, it is written as: $$ H = \sum_k\sum_\sigma \Psi_{k\sigma}^+ h_k\Psi_{k\sigma} $$ where $\Psi_{k\sigma} = [a_{k\uparrow}, b_{k\uparrow}, c_{k\uparrow},a_{k\downarrow}, b_{k\downarrow}, c_{k\downarrow}]^T$; $a_{k\sigma},b_{k\sigma},c_{k\sigma}$ are sublattice in the unitcell, and $h_k$ is $6\times6$ matrix.
We can numerically diagonalize $h_k$ and calculate band-structure, I did it in MATLAB and got (showing here only the lowest band)
Now, I want to calculate spin texture, for example, the magnitude of say z-polarized spin at each point in k-space. One example is given in Figure 8 of Ref: arXiv:2008.10815 as
All in all, I want to know how to calculate $S_z$ at each point in $k$-space when we have wavefunctions and eigenvalues?