# Problems while numerically computing band structure using $\mathbf k\cdot \mathbf p$ theory

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.

• Would Computational Science be a better home for this question? Feb 19, 2015 at 10:12
• Maybe, but the background is quite specifically physics. So I thought it might have the best chance of getting an answer here, where potentially some k.p guys are around.
– DaP
Feb 19, 2015 at 10:14
• Please post your Hamiltonian (or how exactly you construct it) so that others could actually try to reproduce your results and see what's wrong. Feb 23, 2015 at 10:48
• As for distinguishing spin eigenstates, have you considered splitting their energies by some small amount? E.g. if you hard-code some initial $E(\vec k)$ into the Hamiltonian, you could use $E(\vec k)$ for spin-up states and $E(\vec k)+\delta(k)$ for spin-down ones, where $\delta/E\sim10^{-6}$ or like that. Then you'll have almost degenerate result, but will be able to distinguish the final states. Feb 23, 2015 at 11:47
• @Ruslan Thanks for the tip! I updated the question with the Hamiltonian I used and my current progress.
– DaP
Feb 25, 2015 at 8:46

The way to handle the degeneracy problem is to re-map the eight vector basis set given in Eq.(4) of your reference paper. For any direction a re-mapping will reduce the 8 x 8 Hamiltonian to two 4 x 4 blocks. For example, for the Y direction, re-map to the following eight vectors (and normalize all by dividing by the square root of 2): $$u_1 + iu_2$$, $$u_4+iu_5$$, $$-iu_3-u_6$$, $$-u_7+iu_8$$, $$u_1-iu_2$$, $$u_4-iu_5$$, $$iu_3+u_6$$, $$-u_7-iu_8$$