# Hermitian Matrix giving non-real eigenvalues

I am using the following code to find the eigenvalues of a tightbinding matrix

ClearAll["Global*"] (*Clear all variables*)
h := 6.63*10^-34; (* J.s, Planck's constant *)
hbar := (663*10^-34)/(100*2 \[Pi]); (* J.s, Reduced Planck's constant *)
q := 16/10*10^-19 ;  (* Coulomb, electron charge *)
m0 := 91/10*10^-31; (* kg free electron rest mass *)
a := 5.43*10^-10;
vss\[Sigma] := 2.03;
vsp\[Sigma] := 2.55;
vpp\[Sigma] := 4.55;
vpp\[Pi] := 1.09;
esa := -13.55;
esb := -13.55;
epa := -6.52;
epb := -6.52;
n1 := a/4*{1, 1, 1};
n2 := a/4*{-1, -1, 1};
n3 := a/4*{-1, 1, -1};
n4 := a/4*{1, -1, -1};
k = {k1, k1, k1};
g0[k_] = Exp[I*(k.n1)] + Exp[I*(k.n2)] + Exp[I*(k.n3)] + Exp[I*(k.n4)];
g1[k_] = Exp[I*(k.n1)] - Exp[I*(k.n2)] - Exp[I*(k.n3)] + Exp[I*(k.n4)];
g2[k_] = Exp[I*(k.n1)] - Exp[I*(k.n2)] + Exp[I*(k.n3)] - Exp[I*(k.n4)];
g3[k_] = Exp[I*(k.n1)] + Exp[I*(k.n2)] - Exp[I*(k.n3)] - Exp[I*(k.n4)];
v0 = vss\[Sigma];
v1 = 1/Sqrt[3]*vsp\[Sigma];
v2 = 1/3*vpp\[Sigma] - 2/3*vpp\[Pi];
v3 = 1/3*vpp\[Sigma] + 1/3*vpp\[Pi];

H [k_] = {{esa, 0, 0, 0, -v0*g0[k], v1*g1[k], v1*g2[k], v1*g3[k]},
{0, epa, 0, 0, -v1*g1[k], v2*g0[k], v3*g3[k], v3*g2[k]},
{0, 0, epa, 0, -v1*g2[k], v3*g3[k], v2*g0[k], v3*g1[k]},
{0, 0, 0, epa, -v1*g3[k], v3*g2[k], v3*g1[k], v2*g0[k]},
{-v0*g0[k]\[Conjugate], -v1*g1[k]\[Conjugate], -v1*
g2[k]\[Conjugate], -v1*g3[k]\[Conjugate], esb, 0, 0, 0},
{v1*g1[k]\[Conjugate], v2*g0[k]\[Conjugate], v3*g3[k]\[Conjugate],
v3*g2[k]\[Conjugate], 0, epb, 0, 0},
{v1*g2[k]\[Conjugate], v3*g3[k]\[Conjugate], v2*g0[k]\[Conjugate],
v3*g1[k]\[Conjugate], 0, 0, epb, 0},
{v1*g3[k]\[Conjugate], v3*g2[k]\[Conjugate], v3*g1[k]\[Conjugate],
v2*g0[k]\[Conjugate], 0, 0, 0, epb}};


Now, I want that k1 should be real, so that the eigenvalues of H are real. How do I put the condition that k1 is real, and therefore plot the 8 eigenvalues as a function of k1.

• Did you check that your matrix passes HermitianMatrixQ[]? Commented Apr 8, 2018 at 6:11
• You have what is probably a bug: you should use H[k_]:=…, not H[k_]=…. This is because you've defined k above to be {k1,k1,k1}. Additionally you've got some weird use of := in the lines above. Rule of thumb: if it's a function definition, use :=. Otherwise, use =. Commented Apr 8, 2018 at 6:41
• To be clear, you don't use the physical constants defined at the top, do you? The first variable you define which you actually use is a. Commented Apr 8, 2018 at 6:49

Get rid of the exponential terms with some rules.

rule1 = k1 -> j 1000000000000/543;

rule2 = First@Solve[k1 == (k1 /. rule1), j];

rule3 = {r1 -> E^(-((3 I j)/4)), r2 -> E^((3 I j)/4),
r3 -> E^(-(( I j)/4)), r4 -> E^(( I j)/4), r5 -> E^( -I j),
r6 -> E^( I j)};

rule4 = Reverse[#] & /@ rule3;

H[j_] = Rationalize[{{esa, 0, 0, 0, -v0*g0[k], v1*g1[k], v1*g2[k],
v1*g3[k]}, {0, epa, 0, 0, -v1*g1[k], v2*g0[k], v3*g3[k],
v3*g2[k]}, {0, 0, epa, 0, -v1*g2[k], v3*g3[k], v2*g0[k],
v3*g1[k]}, {0, 0, 0, epa, -v1*g3[k], v3*g2[k], v3*g1[k],
v2*g0[k]}, {-v0*g0[k]\[Conjugate], -v1*g1[k]\[Conjugate], -v1*
g2[k]\[Conjugate], -v1*g3[k]\[Conjugate], esb, 0, 0,
0}, {v1*g1[k]\[Conjugate], v2*g0[k]\[Conjugate],
v3*g3[k]\[Conjugate], v3*g2[k]\[Conjugate], 0, epb, 0,
0}, {v1*g2[k]\[Conjugate], v3*g3[k]\[Conjugate],
v2*g0[k]\[Conjugate], v3*g1[k]\[Conjugate], 0, 0, epb,
0}, {v1*g3[k]\[Conjugate], v3*g2[k]\[Conjugate],
v3*g1[k]\[Conjugate], v2*g0[k]\[Conjugate], 0, 0, 0, epb}},
0] /. rule1 // FullSimplify[#, j \[Element] Reals] &;


Now eigenvalues can be found easily.

ev[k1_] = Eigenvalues[H[j] /. rule4] /. rule3 /. rule2 // FullSimplify

(*   {1/100 (-652 - Sqrt[192506 - 92650 Cos[(543 k1)/1000000000000]]),
1/100 (-652 - Sqrt[192506 - 92650 Cos[(543 k1)/1000000000000]]),
1/100 (-652 + Sqrt[192506 - 92650 Cos[(543 k1)/1000000000000]]),
1/100 (-652 + Sqrt[192506 - 92650 Cos[(543 k1)/1000000000000]]),
Root[55675094863665015431373457525894578076350 +
116105129973774631996274824320583157882289 Cos[(543 k1)/
1000000000000] -
3596590815713159387599933503345760461615 Cos[(543 k1)/
500000000000] + (82327817623765464574496363654696641471000 +
23835350206662433354238825284654354098200 Cos[(543 k1)/
1000000000000]) #1 + \
(21517418806718149491518110004046239465000 +
612313303756608508515855923570789780000 Cos[(543 k1)/
1000000000000]) #1^2 +
1763955906068083714962627142616247000000 #1^3 +
43945089837271642126622499816050000000 #1^4 &, 1],
Root[55675094863665015431373457525894578076350 +
116105129973774631996274824320583157882289 Cos[(543 k1)/
1000000000000] -
3596590815713159387599933503345760461615 Cos[(543 k1)/
500000000000] + (82327817623765464574496363654696641471000 +
23835350206662433354238825284654354098200 Cos[(543 k1)/
1000000000000]) #1 + \
(21517418806718149491518110004046239465000 +
612313303756608508515855923570789780000 Cos[(543 k1)/
1000000000000]) #1^2 +
1763955906068083714962627142616247000000 #1^3 +
43945089837271642126622499816050000000 #1^4 &, 2],
Root[55675094863665015431373457525894578076350 +
116105129973774631996274824320583157882289 Cos[(543 k1)/
1000000000000] -
3596590815713159387599933503345760461615 Cos[(543 k1)/
500000000000] + (82327817623765464574496363654696641471000 +
23835350206662433354238825284654354098200 Cos[(543 k1)/
1000000000000]) #1 + \
(21517418806718149491518110004046239465000 +
612313303756608508515855923570789780000 Cos[(543 k1)/
1000000000000]) #1^2 +
1763955906068083714962627142616247000000 #1^3 +
43945089837271642126622499816050000000 #1^4 &, 3],
Root[55675094863665015431373457525894578076350 +
116105129973774631996274824320583157882289 Cos[(543 k1)/
1000000000000] -
3596590815713159387599933503345760461615 Cos[(543 k1)/
500000000000] + (82327817623765464574496363654696641471000 +
23835350206662433354238825284654354098200 Cos[(543 k1)/
1000000000000]) #1 + \
(21517418806718149491518110004046239465000 +
612313303756608508515855923570789780000 Cos[(543 k1)/
1000000000000]) #1^2 +
1763955906068083714962627142616247000000 #1^3 +
43945089837271642126622499816050000000 #1^4 &, 4]}   *)


The plot shows the 6 independent solutions

Plot[Evaluate[ev[k1]], {k1, -2 10^10, 2 10^10},
PlotStyle -> Table[Hue[.8 j/7], {j, 0, 7}]]


• I didn't get what you did there. What do the different rules do? Commented Apr 8, 2018 at 14:44
• rule1 = k1 -> j 1000000000000/543; What is this for? Commented Apr 8, 2018 at 14:47
• To simplify expressions. to get a better overview. Commented Apr 8, 2018 at 19:20
• Substituting exponential terms helps the Eigenvalues ` command a lot. Commented Apr 8, 2018 at 19:24
• The program is taking indefinitely long to run. Commented Apr 9, 2018 at 12:24