# Twisted Bilayer Graphene Mathematica Plots

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.

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


## 1 Answer

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


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


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.

• Thanks a lot@b3m2a1. Jul 23, 2020 at 6:04
• @b3m2a1 how can we correct the avoided crossings issue in this plot? I have not idea about. Your help will be appreciated. Thanks Jun 17, 2023 at 20:45
• @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 Jun 24, 2023 at 0:20