# Continuum model Hamiltonian for TMDS

I am trying to implement the continuum model Hamiltonian for transition metal dichalcogenide systems as described in the reference: PRL 122, 086402 (2019)

   (* Define the constants*)
(* twist angle*)
\[Theta] = 1.2 Degree;
(* For MoTe2 V, psi, w*)
V = 8 (* meV*);
psi = -89.6 ;
w = -8.5 (*meV*);
Vz = -2 (* displacement field meV*);
meff = 0.62 (* valence band effective mass*);
a0 = 3.472 (*angstrom*);
NN = 6(*truncate range*);

\[Theta] =
1.2 Degree; (*\[Theta]\[Pi]/180conversion from degress to radian*)

V = V/1000;
w = w/1000;
Vz = Vz/1000;
psi = psi Degree;
hbar = 4.135667*10^-15/2/\[Pi];
me = 0.51099895*10^6;
meff = meff*me;
kin = -hbar^2/2/meff*9*10^16;
a0 = 3.472 ;
a0 = a0*10^-10;
a1 = a0 {1, 0};
a2 = a0 {Cos[\[Pi]/3], Sin[\[Pi]/3]};
(* Reciprocal lattice vectors*)
G1 = (4 \[Pi])/((Sqrt[3]) a0) {0, 1};
G2 = (4 \[Pi])/((Sqrt[3]) a0) {-Sin[\[Pi]/3], Cos[\[Pi]/3]};
(*G2=(4\[Pi])/((Sqrt[3])a0){Cos[\[Pi]/3],Sin[\[Pi]/3]};*)

g1 = \[Theta] {G1[[2]], -G1[[1]]};
g2 = \[Theta] {G2[[2]], -G2[[1]]};
kD = g1[[1]]/Sqrt[3];

Kb = {-kD*Sqrt[3]/2, -kD/2}; (* K botom layer *)
Kt = {-kD*Sqrt[3]/2,
kD/2} ;   (* K top layer *)

(* Define the lattice*)
n = 6;
nlayer = (2 n + 1)^2;
L = Table[0, {i, 1, 2 (2 n + 1)^2}, {j, 1, 2}];
invL = Table[0, {i, 1, 2 n + 1}, {j, 1, 2 n + 1}];
k = 1;
count = 0;
Do[Do[
(*bottom layer*)
L[[k, 1]] = i; L[[k, 2]] = j;
invL[[i + n + 1, j + n + 1]] = count;
count = count + 1;

(* top layer*)
L[[k + nlayer, 1]] = i; L[[k + nlayer, 2]] = j;
k = k + 1;
, {j, -n, n}], {i, -n, n}]

(* check the form of inVL*)
MatrixForm[invL];

(* Hamiltonian implementation *)
Clear[kx, ky]

HHHH[kx_, ky_] :=
Module[{ h, NN},

NN = 6;
siteN = (2 NN + 1)^2;
h = ConstantArray[0, {2 siteN, 2 siteN}];
Do[
ix = L[[i, 1]];
iy = L[[i, 2]];

(* diagnal elements in the Hamiltonian*)

ax = kx - Kb[[1]] + ix g1[[1]] + iy g2[[1]];
ay = ky - Kb[[2]] + ix g1[[2]] + iy g2[[2]];
h[[i, i]] += kin (ax^2 + ay^2) + Vz/2;

ax = kx - Kt[[1]] + ix g1[[1]] + iy g2[[1]];
ay = ky - Kt[[2]] + ix g1[[2]] + iy g2[[2]];
h[[i + siteN, i + siteN]] += kin (ax^2 + ay^2) - Vz/2;

h[[i + siteN, i]] += w;
h[[i, i + siteN]] += w;

If[ix != NN,
j = invL[[ix + 1 + NN + 1, iy + NN + 1]];
h[[j + 1, i]] += V Exp[I*psi];
h[[j + siteN + 1, i + siteN]] += V Exp[-I*psi];
];

If [iy != -NN,
j = invL[[ix + NN + 1, iy - 1 + NN + 1]];
h[[j + 1, i]] += V Exp[I psi];
h[[j + siteN + 1, i + siteN]] += V*Exp[-I psi];
h[[j + 1, i + siteN]] += w;
];

If [ix != -NN && iy != NN,

j = invL[[ix - 1 + NN + 1, iy + 1 + NN + 1]];
h[[j + 1, i]] += V*Exp[I*psi];
h[[j + siteN + 1, i + siteN]] += V*Exp[-I*psi];
h[[j + siteN + 1, i]] += w;
];

If [ix != -NN,
j = invL[[ix - 1 + NN + 1, iy + NN + 1]];
h[[j + 1, i]] += V*Exp[-I*psi];
h[[j + siteN + 1, i + siteN]] += V*Exp[I*psi];
];
If [iy != NN,
j = invL[[ix + NN + 1, iy + 1 + NN + 1]];
h[[j + 1, i]] += V*Exp[-I*psi];
h[[j + siteN + 1, i + siteN]] += V*Exp[I*psi];
h[[j + siteN + 1, i]] += w;
];

If [ix != NN && iy != -NN,

j = invL[[ix + 1 + NN + 1, iy - 1 + NN + 1]];
h[[j + 1, i]] += V*Exp[-I*psi];
h[[j + siteN + 1, i + siteN]] += V*Exp[I*psi];
h[[j + 1, i + siteN]] += w;];

, {i, 1, siteN}];

Return[h]

]

Num = 50;
KXX1 = N@Subdivide[g1[[1]]/2, 0, Num - 1];
KXX2 = N@Subdivide[0, -g1[[1]]/2, Num - 1];
KXX3 = N@Subdivide[-g1[[1]]/2, -g1[[1]]/2, Num - 1];
KXX4 = N@Subdivide[-g1[[1]]/2, g1[[1]]/2, (2 Num - 1)];
KX = Join[KXX1, KXX2, KXX3, KXX4];

KYY1 = N@Subdivide[-g1[[1]]/(Sqrt [3]/2), 0, Num - 1];
KYY2 = N@Subdivide[0, g1[[1]]/(Sqrt[3]/2), Num - 1];
KYY3 = N@
Subdivide[g1[[1]]/(Sqrt[3]/2), -g1[[1]]/(Sqrt [3]/2), Num - 1];
KYY4 = N@
Subdivide[-g1[[1]]/(Sqrt [3]/2), -g1[[1]]/(Sqrt [3]/2), (2 Num -
1)];
KY = Join[KYY1, KYY2, KYY3, KYY4];

Eigen = Table[0, {i1, Length[KX]}, {j1, 2 siteN}];

Do[Eigen[[k]] = Eigenvalues[HHHH[KX[[k]], KY[[k]]]] // Sort, {k, 1,
Length[KX]}];

out = ParallelTable[
Eigenvalues[HHHH[KX[[k]], KY[[k]]]] // Sort, {k, 1, Length[KX]}];
ListLinePlot[Transpose[out], PlotRange -> {{0, 250}, {0, 0.025}}]


The resulting band structure is look like this

I recently came across the github code which reproduce the band structure from the reference. (https://github.com/zihaophys/twisted_bilayer_graphene/blob/master/TMD/homobilayer/moireTMD.py)

I have studied and trying to reproduce in exact same way. Can someone please point out how to fix the band structure.

Thanks

• Welcome to the Mathematica Stack Exchange. The code you have uploaded does not evaluate or plot anything (yet). This stack site is about the technical software called Mathematica and the associated Wolfram Language. Stack sites are not project assistance sites. Please present a minimal example and ask an on-topic question that can be answered. Where exactly are you having difficulty?
– Syed
Mar 17, 2023 at 4:49
• There several posts about twisted bilayer graphene - see mathematica.stackexchange.com/search?q=graphene Mar 18, 2023 at 6:50
• Dear All, I was trying to get the continuum band structure for TMDs. So I did implement the Hamiltonian. My main problem is getting the band structure along the given high symmetry line. The band structure getting from this code is different than the one in the reference. I'm thinking may be the way I define the high symmetry line is not correct. Will look into you. Yes I did refer the examples about twisted bilayer system. Meantime I would really appreciate if someone can give me a comment who is familiar with continuum models of TMDs. Thanks Mar 20, 2023 at 0:32