Discontinuities in eigenvalues plotting with Plot3D

I have a 6x6 matrix, depending on 2 variables kx and ky and defined as the sum of the two following matrices (t1, t2 and b are some numerical parameters)

H1[kx_,ky_]:=t1*{{0, Exp[I*kx*b], 0, Exp[I*ky*b], 0, 0},
{Exp[-I*kx*b], 0, Exp[I*ky*b], 0, 0, 0},
{0, Exp[-I*ky*b], 0, Exp[-I*kx*b], 0, 0},
{Exp[-I*ky*b], 0, Exp[I*kx*b], 0, 0, 0},
{0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0}}

H2[kx_, ky_] := t2*{{0, 0, 0, 0, Exp[I*(-kx*Sqrt[3]/2 + ky/2)*b], Exp[I*(kx/2 - ky*Sqrt[3]/2)*b]},
{0, 0, 0, 0, Exp[I*(kx*Sqrt[3]/2 + ky/2)*b], Exp[I*(-kx/2 - ky*Sqrt[3]/2)*b]},
{0, 0, 0, 0, Exp[I*(kx*Sqrt[3]/2 - ky/2)*b], Exp[I*(-kx/2 + ky*Sqrt[3]/2)*b]},
{0, 0, 0, 0, Exp[I*(-kx*Sqrt[3]/2 - ky/2)*b], Exp[I*(kx/2 + ky*Sqrt[3]/2)*b]}, {Exp[-I*(-kx*Sqrt[3]/2 + ky/2)*b], Exp[-I*(kx*Sqrt[3]/2 + ky/2)*b], Exp[-I*(kx*Sqrt[3]/2 - ky/2)*b], Exp[-I*(-kx*Sqrt[3]/2 - ky/2)*b], 0, 0},
{Exp[-I*(kx/2 - ky*Sqrt[3]/2)*b], Exp[-I*(-kx/2 - ky*Sqrt[3]/2)*b], Exp[-I*(-kx/2 + ky*Sqrt[3]/2)*b], Exp[-I*(kx/2 + ky*Sqrt[3]/2)*b], 0, 0}}

H[kx_, ky_] := H1[kx,ky] + H2[kx,ky]


The matrix is Hermitian, and I need to obtain the 6 real eigenvalues as functions of kx and ky. There's some physical meaning behind those eigenvalues, and they should be continuous functions of kx and ky. So I do

ev[kx_, ky_] := Eigenvalues[H[kx, ky]]


However, when I plot the eigenvalues with Plot3D, I clearly get some discontinuities and very messy plots

Plot3D[{ev[kx, ky][[4]]}, {kx, -2*Pi, 2*Pi}, {ky, -2*Pi, 2*Pi}]


I do not think there are typos in the matrix expressions as I have other ways to check the numerical values of some eigenvalues, and the matrix is hermitian as expected.

Following a similar answered question here on SE, I've tried to use

ev[kx, ky][[4]] // Re


but it does not fix the problem. Am I doing something wrong, or how can I fix that?

• You need to fix H2. I don't get a 6 by 6 matrix. Also, give us some good values for b, t1, t2 and y. Finally are the surfaces formed by the eigenvalues intersecting? If so how do we know how to define a surface? – Hugh Sep 23 '16 at 18:45
• Woops, actually I messed up when formatting the H2, in my notebook H2 is ok. It should be a 6x6 matrix now. For t1 and t2, typical values would be -1 and -1. For b, it needs to be 1/(1+sqrt(3)). kx and ky do not matter, eigenvalues will be periodic functions. Concerning your last question, actually I dont now. When I try to plot several of the eigenvalues the plot becomes very messy with all the discontinuities. The point is that the eigenvalues represent the energy levels of some physical systems, so they need to be continuous surfaces. – Castou Sep 26 '16 at 1:17
• Related: How to plot several functions without jumping? (multiple eigenvalues of a system as functions of 2 parameters). I believe this is the key problem: "...during the plotting of the eigenvalues, for certain regions of the parameter space, Mathematica organizes the 8 values differently, so the different sheets swap places, so we get tons of random garbage vertical bits between sheets." Not sure there is a foolproof solution, but you can take a look at the answers there. – Rahul Sep 26 '16 at 2:59

1 Answer

Using the parameters you provided I played around with the calculations. Here are your functions again.

ClearAll[vals, H1, H2];
vals = {t1 -> -1, t2 -> -1, b -> 1/(1 + Sqrt[3])};
H1[kx_, ky_] :=
t1*{{0, Exp[I*kx*b], 0, Exp[I*ky*b], 0, 0}, {Exp[-I*kx*b], 0,
Exp[I*ky*b], 0, 0, 0}, {0, Exp[-I*ky*b], 0, Exp[-I*kx*b], 0,
0}, {Exp[-I*ky*b], 0, Exp[I*kx*b], 0, 0, 0}, {0, 0, 0, 0, 0,
0}, {0, 0, 0, 0, 0, 0}} /. vals;

H2[kx_, ky_] :=
t2*{{0, 0, 0, 0, Exp[I*(-kx*Sqrt[3]/2 + ky/2)*b],
Exp[I*(kx/2 - ky*Sqrt[3]/2)*b]}, {0, 0, 0, 0,
Exp[I*(kx*Sqrt[3]/2 + ky/2)*b],
Exp[I*(-kx/2 - ky*Sqrt[3]/2)*b]}, {0, 0, 0, 0,
Exp[I*(kx*Sqrt[3]/2 - ky/2)*b],
Exp[I*(-kx/2 + ky*Sqrt[3]/2)*b]}, {0, 0, 0, 0,
Exp[I*(-kx*Sqrt[3]/2 - ky/2)*b],
Exp[I*(kx/2 + ky*Sqrt[3]/2)*b]}, {Exp[-I*(-kx*Sqrt[3]/2 + ky/2)*
b], Exp[-I*(kx*Sqrt[3]/2 + ky/2)*b],
Exp[-I*(kx*Sqrt[3]/2 - ky/2)*b],
Exp[-I*(-kx*Sqrt[3]/2 - ky/2)*b], 0,
0}, {Exp[-I*(kx/2 - ky*Sqrt[3]/2)*b],
Exp[-I*(-kx/2 - ky*Sqrt[3]/2)*b],
Exp[-I*(-kx/2 + ky*Sqrt[3]/2)*b],
Exp[-I*(kx/2 + ky*Sqrt[3]/2)*b], 0, 0}} /. vals;

H[kx_, ky_] := H1[kx, ky] + H2[kx, ky];
ev[kx_, ky_] := Eigenvalues[H[kx, ky]];


If we plot the eigenvalues then things are messed up since they are not ordered satisfactory.

Plot3D[ev[kx, ky], {kx, -2*Pi, 2*Pi}, {ky, -2*Pi, 2*Pi},
BoxRatios -> {1, 1, 2}]


However there are clearly constant eigenvalues at 0 and 2. The rest seem to appear in bands.

I then decided to sort them and did this

a = Table[
Sort@ev[kx, ky], {kx, -2*Pi, 2*Pi, Pi/10.}, {ky, -2*Pi, 2*Pi,
Pi/10.}];

Table[ListPlot3D[a[[All, All, n]]], {n, 6}]


Giving

This looks much better to my eyes. Does this help?