I want to use k.p theory to numerically compute the band structure of a bulk semiconductor. The band I like to include are the lowest conduction band (cb), the heavy-hole (hh), the light-hole (lh) and the split-off (so) band, each twice spin-degenerate. To do this, I have to find the eigenvalues and eigenvectors of an 8x8 Hamiltonian. I know there are plenty of different types for the very same problem depending on the definitions of the states at $k=0$, so I do not show this Hamiltonian here.
The lattice periodic wave function of band $\lambda$ can be written as an expansion into the states at $k=0$:
$ u_{\lambda,{\bf k}} = \sum_\nu c_{\lambda,\nu,{\bf k}} u_{\lambda,{\bf 0}} $.
To numerically get the band structure and the eigenvectors which contain the expansion coefficients $c_{\lambda,\nu,{\bf k}}$, I use an eigenproblem solver like zgeev (Lapack) or ComplexEigenSolver (c++ eigen library) for each k-point (in a certain direction). These routines give you the eigenvalues in no particular order so I sort them (and the corresponding eigenvectors) in descending order (for each k also). I can then say, that the first two values are the two spin degenerate cb, the next two the hh, then lh, then so - I just do not know which one is spin up and which is spin down (since they are degenerate). The band structure I get looks good (the spin degenerate bands have the same energies):
My problem now are the expansions coefficients, i.e. the eigenvectors, from which I get 8 and each have 8 entries. These entries have a particular order, depending on how my Hamiltonian is build. Lets assume the eigenvector ${\bf a}_{\lambda,{\bf k}}$ for band $\lambda$ and momentum $\bf k$ has the following order:
${\bf a}_{\lambda,{\bf k}} = (c_{\lambda,cb\uparrow,{\bf k}},c_{\lambda,cb\downarrow,{\bf k}},c_{\lambda,hh\uparrow,{\bf k}},c_{\lambda,hh\downarrow,{\bf k}},c_{\lambda,lh\uparrow,{\bf k}},c_{\lambda,lh\downarrow,{\bf k}},c_{\lambda,so\uparrow,{\bf k}},c_{\lambda,so\downarrow,{\bf k}})$
Since I only know the two vectors that are either $\lambda=cb\uparrow$ or $\lambda=cb\downarrow$, I can (and will) confuse the two spin coefficients. How can I sort them correctly? I will need some additional info for that and I do not know one.
Another problem is how my (unsorted) coefficients look like (shown are the 4 possible cb to cb coefficients):
Since they are not sorted by spin, I would not assume a smooth function here but when I plot only points, I would assume to see something like a smooth function if I interchange the colors (black and blue as well as red and green are interchangable). But, this is far from being smooth, and whats about this single k point with an imaginary part?
Long story short: How can I sort my spin degenerate eigenvectors and eigenvalues correctly? And what could be a reason for my eigenvectors to have such a weird behavior? I checked my implementation really hard and double checked against another from a colleague. Any tips/hints appreciated.
EDIT
The construction of the Hamiltonian and the appropriate lattice periodic basis functions can be found e.g. on this page. While this describes what I want to do in the end (Quantum Well Structure with strain effects etc.) I need to understand the easier bulk problem first completely. So the basis set is given by Eq. (4), the Hamiltonian by Eq. (5), and its subparts by Eqs. (6), (8), (17), (18), (20), (22). Using their notation for the basis states, that would give the eigenvector the following order:
${\bf a}_{\lambda,{\bf k}} = (c_{\lambda,u_1,{\bf k}},c_{\lambda,u_2,{\bf k}},c_{\lambda,u_3,{\bf k}},c_{\lambda,u_4,{\bf k}},c_{\lambda,u_5,{\bf k}},c_{\lambda,u_6,{\bf k}},c_{\lambda,u_7,{\bf k}},c_{\lambda,u_8,{\bf k}})$
where the pairs $u_1,u_2$, $u_3,u_6$, $u_4,u_5$ and $u_6,u_7$ are doubly degenerate. I tried to add a slight energy $\epsilon=10^{-6}eV$ to the first one of each pair to split the degeneracy (as suggested by Ruslan), but it did not work properly. BUT: If I add $\epsilon$ or multiple of $\epsilon$ not necessarily to the first one of each pair but to both, it works! This is what I get:
Here, I plot the real part of the coefficients (imaginary part is zero) so that I also see the sign. Since the eigenvector can be flipped due to numerics, I always multiplied the eigenvector ${\bf a}_{\lambda,\nu,k}$ whenever its component ${\bf a}_{\lambda,\lambda,k}<0$. This is actually the assumption that its "primary" component is always positive, is this a valid assumption? The description of the states in my figure follows $|l,s,j,m_j>$. I added the following $\epsilon$ to the diagonal parts (following the notation from the page given and in superscript the indices of the matrix):
$-\epsilon$ to $H_{cc}^{11}$
$-\epsilon$ to $H_{ss}^{22}$
$-5\epsilon$ to $H_{vv}^{22}$
$-10\epsilon$ to $H_{vv}^{33}$
$-\epsilon$ to $H_{vv}^{44}$
This seems quite arbitrary to me, no real system behind it. My only remaining question then is, can I be sure that this combination of adding $\epsilon$ will work for other material systems (I used GaAs parameters) and directions (I plotted $\Gamma \to X$)? As far as I have tested it for other parameters it seems to work.
