0
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Did you check that your matrix passes HermitianMatrixQ[]? $\endgroup$ – J. M. is away Apr 8 '18 at 6:11
  • $\begingroup$ 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 =. $\endgroup$ – Patrick Stevens Apr 8 '18 at 6:41
  • $\begingroup$ 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. $\endgroup$ – Patrick Stevens Apr 8 '18 at 6:49
1
$\begingroup$

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

enter image description here

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.