Hermitian Matrix giving non-real eigenvalues

Question

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.

Answer 1

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}]]