3
$\begingroup$

I am trying to plot band structure of twisted bilayer graphene. This is a ten bands model studied in this work reference. Everything is working fine but the bands are taking jump among each other for no reason. I have tried some solutions but none of them have worked. If anyone have a nicer approach that also takes optimum time to plot these bands. Because with these bad visualization it takes almost one minute to plot it. This is the code. The other approaches could be calling, Method -> {"Arnoldi", "Criteria" -> "RealPart"}. But this only works for 8 bands not for 10 bands. enter image description hereenter image description here

taz = 0;
tp = 0.003;
tpp = 0.004;
tmp = 0.004;
tk = 0;
tkb = 0;
tpppz = 0.016;
tmppz = 0.016;
tpkp = 0.016;
tmkp = -0.016;
\[Mu]az = -0.1 - 6*taz;
\[Mu]aa = 0.11 - 4 (tk + tkb);
\[Mu]dw = 0.11 - 4 (tk + tkb);
A = 0.110;
b = 0.033;
c = 0.033;
d = 0.573;
g = \[Pi]/8;
te = I;
AA[px_, py_, a_] := 
te*Exp[-I*g] (1 + \[Phi]1[px, py, a, 1] + \[Phi]2[px, py, a, -1]);
AA1[px_, py_, a_] := Conjugate[AA[px, py, a]]
\[Phi]1[px_, py_, a_, l_] := Exp[-I*l*Sqrt[3]/2 a (-px + Sqrt[3] py)]
\[Phi]2[px_, py_, a_, l_] := Exp[-I*l*Sqrt[3]/2 a (px + Sqrt[3] py)]
\[Phi]3[px_, py_, a_, l_] := Exp[-I*l*3 a*py]
W[l_] := Exp[I*l*(2*\[Pi])/3]
\[Xi][l_] := Exp[I*l*(2*\[Pi])/6]
H1[px_, py_, a_] := (\[Phi]2[px, py, a, 1] + \[Phi]3[px, py, a, 1] + \[Phi]1[px, 
py, a, 1] + \[Phi]2[px, py, a, -1] + \[Phi]3[px, py, 
a, -1] + \[Phi]1[px, py, a, -1])
R1[px_, py_, a_] := {0, 0, 0, 0, 0, 
0, -(W[1] + \[Phi]3[px, py, a, 1]*W[-1] + \[Phi]1[px, py, a, 
   1]) \[Xi][-1]*
A, (1 + \[Phi]3[px, py, a, 1]*W[-1] + \[Phi]1[px, py, a, 1]*
  W[1]) \[Xi][-1]*
A, (W[1] + \[Phi]2[px, py, a, 1]*W[-1] + \[Phi]3[px, py, a, 
  1]) \[Xi][1]*
A, -(W[1] + \[Phi]2[px, py, a, 1] + \[Phi]3[px, py, a, 1]*
   W[-1]) \[Xi][1]*A}
R2[px_, py_, a_] := {0, 0, 0, 0, 0, 
0, (W[-1] + \[Phi]3[px, py, a, 1]*W[1] + \[Phi]1[px, py, a, 
  1]) \[Xi][1]*
b, (1 + \[Phi]3[px, py, a, 1] + \[Phi]1[px, py, a, 
  1]) c, (W[-1] + \[Phi]2[px, py, a, 1]*W[1] + \[Phi]3[px, py, a, 
  1]) \[Xi][-1]*
b, (1 + \[Phi]2[px, py, a, 1] + \[Phi]3[px, py, a, 1]) c}
R3[px_, py_, a_] := {0, 0, 0, 0, 0, 
0, (1 + \[Phi]3[px, py, a, 1] + \[Phi]1[px, py, a, 
  1]) c, (1 + \[Phi]3[px, py, a, 1]*W[1] + \[Phi]1[px, py, a, 1]*
  W[-1]) \[Xi][1]*
b, (1 + \[Phi]2[px, py, a, 1] + \[Phi]3[px, py, a, 
  1]) c, (W[-1] + \[Phi]2[px, py, a, 1] + \[Phi]3[px, py, a, 1]*
  W[1]) \[Xi][-1]*b}
R4[px_, py_, a_] := {0, 0, 0, 0, 0, 
0, -I*d*\[Phi]2[px, py, a, -1], -I*d*\[Phi]2[px, py, a, -1], I*d, 
I*d}
R5[px_, py_, a_] := {0, 0, 0, 0, 0, 0, -I*d*W[1], -I*d*W[-1], 
I*d*W[1], I*d*W[-1]}
R6[px_, py_, a_] := {0, 0, 0, 0, 0, 
0, -I*W[-1]*d*\[Phi]1[px, py, a, 1], -I*W[1]*
d*\[Phi]1[px, py, a, 1], I*W[-1]*d, I*W[1]*d}
R7[px_, py_, 
a_] := {-(W[-1] + \[Phi]3[px, py, a, -1]*W[1] + \[Phi]1[px, py, 
   a, -1]) \[Xi][
1] A, (W[1] + \[Phi]3[px, py, a, -1]*W[-1] + \[Phi]1[px, py, 
  a, -1]) b*\[Xi][-1], (1 + \[Phi]3[px, py, a, -1] + \[Phi]1[px, 
  py, a, -1]) c, I*d*\[Phi]2[px, py, a, 1], I*W[-1]*d, 
I*W[1]*d*\[Phi]1[px, py, a, -1], 0, 0, 0, 0}
R8[px_, py_, 
a_] := {(1 + \[Phi]3[px, py, a, -1]*W[1] + \[Phi]1[px, py, a, -1]*
  W[-1]) \[Xi][
1] A, (1 + \[Phi]3[px, py, a, -1] + \[Phi]1[px, py, 
  a, -1]) c, (1 + W[-1]*\[Phi]3[px, py, a, -1] + 
 W[1]*\[Phi]1[px, py, a, -1]) \[Xi][-1]*b, 
I*d*\[Phi]2[px, py, a, 1], I*W[1]*d, 
I*W[-1]*d*\[Phi]1[px, py, a, -1], 0, 0, 0, 0}
R9[px_, py_, 
a_] := {(W[-1] + \[Phi]2[px, py, a, -1]*W[1] + \[Phi]3[px, py, 
  a, -1]) \[Xi][-1] A, (W[1] + \[Phi]2[px, py, a, -1]*
  W[-1] + \[Phi]3[px, py, a, -1]) b*\[Xi][
1], (1 + \[Phi]2[px, py, a, -1] + \[Phi]3[px, py, a, -1]) c, -I*
d, -I*W[-1]*d, -I*W[1]*d, 0, 0, 0, 0}
R10[px_, py_, 
a_] := {-(W[-1] + \[Phi]2[px, py, a, -1] + \[Phi]3[px, py, a, -1]*
   W[1]) \[Xi][-1] A, (1 + \[Phi]2[px, py, a, -1] + \[Phi]3[px, 
  py, a, -1]) c, (W[1] + \[Phi]2[px, py, 
  a, -1] + \[Phi]3[px, py, a, -1]*W[-1]) b*\[Xi][1], -I*d, -I*
W[1]*d, -I*W[-1]*d, 0, 0, 0, 0}
r1[px_, py_, 
a_] := {\[Mu]az, (-I*
 tpppz (\[Phi]1[px, py, a, -1] + \[Phi]3[px, py, a, 1]*
    W[-1] + \[Phi]2[px, py, a, -1] W[1]) - (-1) I*
 tmppz (\[Phi]1[px, py, a, 1] + \[Phi]3[px, py, a, -1]*
    W[-1] + \[Phi]2[px, py, a, 1] W[1])), (-(-1) I*
 tpppz (\[Phi]1[px, py, a, 1] + \[Phi]3[px, py, a, -1]*
    W[1] + \[Phi]2[px, py, a, 1] W[-1]) - (-1) (-1) I*
 tmppz (\[Phi]1[px, py, a, -1] + \[Phi]3[px, py, a, 1]*
    W[1] + \[Phi]2[px, py, a, -1] W[-1])), 0, 0, 0, 0, 0, 0, 0}
r2[px_, py_, 
a_] := {(I*
 tpppz (\[Phi]1[px, py, a, 1] + \[Phi]3[px, py, a, -1]*
    W[1] + \[Phi]2[px, py, a, 1] W[-1]) - 
I*tmppz (\[Phi]1[px, py, a, -1] + \[Phi]3[px, py, a, 1]*
    W[1] + \[Phi]2[px, py, a, -1] W[-1])), 
tp*H1[px, py, a] + \[Mu]aa, 
tpp (\[Phi]1[px, py, 
   a, -1] + \[Phi]3[px, py, a, 1] W[-1] + \[Phi]2[px, py, 
    a, -1] W[1]) + 
tmp (\[Phi]1[px, py, a, 
   1] + \[Phi]3[px, py, a, -1] W[-1] + \[Phi]2[px, py, a, 1] W[
    1]), tpkp*\[Phi]2[px, py, a, 1] - tmkp*\[Phi]3[px, py, a, 1], 
tpkp*\[Phi]3[px, py, a, 1]*W[1] - tmkp*W[1], 
tpkp*W[-1] - tmkp*\[Phi]2[px, py, a, 1]*W[-1], 0, 0, 0, 0}
r3[px_, py_, 
a_] := {(-I*
 tpppz (\[Phi]1[px, py, a, -1] + \[Phi]3[px, py, a, 1]*
    W[-1] + \[Phi]2[px, py, a, -1] W[1]) - (-1) I*
 tmppz (\[Phi]1[px, py, a, 1] + \[Phi]3[px, py, a, -1]*
    W[-1] + \[Phi]2[px, py, a, 1] W[1])), 
tpp (\[Phi]1[px, py, a, 
   1] + \[Phi]3[px, py, a, -1] W[1] + \[Phi]2[px, py, a, 
    1] W[-1]) + 
tmp (\[Phi]1[px, py, 
   a, -1] + \[Phi]3[px, py, a, 1] W[1] + \[Phi]2[px, py, 
    a, -1] W[-1]), tp*H1[px, py, a] + \[Mu]aa, 
tpkp*\[Phi]3[px, py, a, 1] - tmkp*\[Phi]2[px, py, a, 1], 
tpkp*W[-1] - tmkp*\[Phi]3[px, py, a, 1]*W[-1], 
tpkp*\[Phi]2[px, py, a, 1]*W[1] - tmkp*W[1], 0, 0, 0, 0}
r4[px_, py_, a_] := {0, 
tpkp*\[Phi]2[px, py, a, -1] - tmkp*\[Phi]3[px, py, a, -1], 
tpkp*\[Phi]3[px, py, a, -1] - tmkp*\[Phi]2[px, py, a, -1], \[Mu]dw, 
0, 0, 0, 0, 0, 0}
r5[px_, py_, a_] := {0, 
tpkp*\[Phi]3[px, py, a, -1]*W[-1] - tmkp*W[-1], 
tpkp*W[1] - tmkp*\[Phi]3[px, py, a, -1]*W[1], 0, \[Mu]dw, 0, 0, 0, 
0, 0}
r6[px_, py_, a_] := {0, tpkp*W[1] - tmkp*\[Phi]2[px, py, a, -1]*W[1], 
tpkp*\[Phi]2[px, py, a, -1]*W[-1] - tmkp*W[-1], 0, 0, \[Mu]dw, 0, 0,
0, 0}
r7[px_, py_, a_] := {0, 0, 0, 0, 0, 0, 0, 0, AA[px, py, a], 0}
r8[px_, py_, a_] := {0, 0, 0, 0, 0, 0, 0, 0, 0, AA[px, py, a]}
r9[px_, py_, a_] := {0, 0, 0, 0, 0, 0, AA1[px, py, a], 0, 0, 0}
r10[px_, py_, a_] := {0, 0, 0, 0, 0, 0, 0, AA1[px, py, a], 0, 0}
hh[px_, py_, a_] := {R1[px, py, a] + r1[px, py, a], 
R2[px, py, a] + r2[px, py, a], R3[px, py, a] + r3[px, py, a], 
R4[px, py, a] + r4[px, py, a], R5[px, py, a] + r5[px, py, a], 
R6[px, py, a] + r6[px, py, a], R7[px, py, a] + r7[px, py, a], 
R8[px, py, a] + r8[px, py, a], R9[px, py, a] + r9[px, py, a], 
R10[px, py, a] + r10[px, py, a]}

ev4 = Plot[
Eigenvalues[
hh[Cos[\[Theta]], Sin[\[Theta]], 1.4]], {\[Theta], -0.5 \[Pi], 
1.5 \[Pi]}]
$\endgroup$

1 Answer 1

7
$\begingroup$

I can provide an answer for the fixing the visual. We'll use some dynamically scoped caching + add a function that only evaluates once Plot is ready to go:

Block[{getEv},
 getEv[θ_?NumericQ] :=
  getEv[θ] = Sort@Eigenvalues[hh[Cos[θ], Sin[θ], 1.4]];
 getEv[θ_?NumericQ, n_] :=
  getEv[θ][[n]];
 Plot[Evaluate@Table[getEv[θ, n], {n, 10}], {θ, -0.5 π, 1.5 π}]
 ]

enter image description here

The jumps happen because Mathematica returns the eigenvalues ordered by magnitude, i.e. it ignores the sign. You can plot the sorted eigenvalues to smooth that out

Block[{getEv},
 getEv[θ_?NumericQ] :=
  getEv[θ] = Sort@Eigenvalues[hh[Cos[θ], Sin[θ], 1.4]];
 getEv[θ_?NumericQ, n_] :=
  getEv[θ][[n]];
 Plot[Evaluate@Table[getEv[θ, n], {n, 10}], {θ, -0.5 π, 1.5 π}]
 ]

enter image description here

Note that this leads to avoided crossings, where the bands could actually switch their energy ordering. I'll leave it up to you how you want to resolve this.

$\endgroup$
3
  • $\begingroup$ Thanks a lot@b3m2a1. $\endgroup$ Jul 23, 2020 at 6:04
  • $\begingroup$ @b3m2a1 how can we correct the avoided crossings issue in this plot? I have not idea about. Your help will be appreciated. Thanks $\endgroup$
    – user199
    Jun 17, 2023 at 20:45
  • $\begingroup$ @user199 you'll need to either detect them after the fact or note that they're coming from working in an adiabatic picture and you need a transformation to a diabatic picture, which will require updating the underlying model to include non-adiabatic couplings $\endgroup$
    – b3m2a1
    Jun 24, 2023 at 0:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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