Twisted Bilayer Graphene Mathematica Plots
Mathematica Asked on August 27, 2021
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]}]
One 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.
Correct answer by b3m2a1 on August 27, 2021
Add your own answers!
Ask a Question
Get help from others!
Recent Answers
- Jon Church on Why fry rice before boiling?
- Joshua Engel on Why fry rice before boiling?
- Lex on Does Google Analytics track 404 page responses as valid page views?
- haakon.io on Why fry rice before boiling?
- Peter Machado on Why fry rice before boiling?
Recent Questions
- How can I transform graph image into a tikzpicture LaTeX code?
- How Do I Get The Ifruit App Off Of Gta 5 / Grand Theft Auto 5
- Iv’e designed a space elevator using a series of lasers. do you know anybody i could submit the designs too that could manufacture the concept and put it to use
- Need help finding a book. Female OP protagonist, magic
- Why is the WWF pending games (“Your turn”) area replaced w/ a column of “Bonus & Reward”gift boxes?