I kinda figured you'd need more info. But I didn't want to spam everyone right away!
Here is everything I have with annotations. In desperation I have also attached the file. Any tips at all would be appreciated.
P.S. I should say that the end goal of this will be to plot that value of T[n,a,d,kx,l,g] {k,-10,10} without it taking infinite time.
Remove["Global`*"];
$Assumptions = a > 0 && d > 0 && l > 0 && g > 0 && n \[Element] Integers;
Hamiltonian for region with SOI and Zeeman field:
Hsoz[qx_, l_, g_] = {{-g, 0, px + I py, -2 I l}, {0, g, 0, px + I py}, {px - I py, 0, -g,
0}, {2 I l, px - I py, 0, g}}
Evals[qx_, l_, g_] = Simplify[Eigenvalues[Hsoz[qx, l, g]]];
Ez1[qx_, l_, g_] = Evals[qx, l, g][[1]];
Ez2[qx_, l_, g_] = Evals[qx, l, g][[2]];
Ez3[qx_, l_, g_] = Evals[qx, l, g][[3]];
Ez4[qx_, l_, g_] = Evals[qx, l, g][[4]];
Evecs[qx_, l_, g_] = Simplify[Eigenvectors[Hsoz[qx, l, g]]];
Normalization below (I suspect there is a better way to do this)
n1[qx_, l_, g_] = g^2 + l^2 + Sqrt[l^4 + g^2 qx^2 + l^2 qx^2];
n2[qx_, l_, g_] = g^2 + l^2 - Sqrt[l^4 + g^2 qx^2 + l^2 qx^2];
Nz1[qx_, l_, g_] =
Simplify[1. + (Abs[l]^2*Abs[qx]^2 + Abs[l*(Ez1[qx, l, g] + g)]^2)/
Abs[Ez1[qx, l, g]*g + n2[qx, l, g]]^2 +
Abs[Ez1[qx, l, g]*(n1[qx, l, g] - g^2) +
g*(n1[qx, l, g] - g^2 - qx^2)]^2/(
Abs[Ez1[qx, l, g]*g + n2[qx, l, g]]^2*Abs[qx]^2)];
Nz2[qx_, l_, g_] =
Simplify[1. + (Abs[l]^2*Abs[qx]^2 + Abs[l*(Ez2[qx, l, g] + g)]^2)/
Abs[Ez2[qx, l, g]*g + n2[qx, l, g]]^2 +
Abs[Ez2[qx, l, g]*(n1[qx, l, g] - g^2) +
g*(n1[qx, l, g] - g^2 - qx^2)]^2/(
Abs[Ez2[qx, l, g]*g + n2[qx, l, g]]^2*Abs[qx]^2)];
Nz3[qx_, l_, g_] =
Simplify[1. + (Abs[l]^2*Abs[qx]^2 + Abs[l*(Ez3[qx, l, g] + g)]^2)/
Abs[Ez3[qx, l, g]*g + n1[qx, l, g]]^2 +
Abs[Ez3[qx, l, g]*(n2[qx, l, g] - g^2) +
g*(n2[qx, l, g] - g^2 - qx^2)]^2/(
Abs[Ez3[qx, l, g]*g + n1[qx, l, g]]^2*Abs[qx]^2)];
Nz4[qx_, l_, g_] =
Simplify[1. + (Abs[l]^2*Abs[qx]^2 + Abs[l*(Ez4[qx, l, g] + g)]^2)/
Abs[Ez4[qx, l, g]*g + n1[qx, l, g]]^2 +
Abs[Ez4[qx, l, g]*(n2[qx, l, g] - g^2) +
g*(n2[qx, l, g] - g^2 - qx^2)]^2/(
Abs[Ez4[qx, l, g]*g + n1[qx, l, g]]^2*Abs[qx]^2)];
wz1[qx_, l_, g_] =
Simplify[(1./Sqrt[Nz1[qx, l, g]]) Evecs[qx, l, g][[1]]];
wz2[qx_, l_, g_] =
Simplify[(1./Sqrt[Nz2[qx, l, g]]) Evecs[qx, l, g][[2]]];
wz3[qx_, l_, g_] =
Simplify[(1./Sqrt[Nz3[qx, l, g]]) Evecs[qx, l, g][[3]]];
wz4[qx_, l_, g_] =
Simplify[(1./Sqrt[Nz4[qx, l, g]]) Evecs[qx, l, g][[4]]];
Conservation of energy:
qz3x[kx_, l_, g_] = Sqrt[
g^2 + kx^2 - 2 Sqrt[g^2 kx^2 - g^2 l^2 + kx^2 l^2]];
qz4x[kx_, l_, g_] = Sqrt[
g^2 + kx^2 + 2 Sqrt[g^2 kx^2 - g^2 l^2 + kx^2 l^2]];
Generic wave equation constructed:
M1[x_, kx_, l_, g_] =
Simplify[wz3[qz3x[kx, l, g], l, g]*E^(I qz3x[kx, l, g] x)];
M2[x_, kx_, l_, g_] =
Simplify[wz4[qz4x[kx, l, g], l, g]*E^(I qz4x[kx, l, g] x)];
M3[x_, kx_, l_, g_] =
Simplify[wz3[-qz3x[kx, l, g], l, g]*E^(-I qz3x[kx, l, g] x)];
M4[x_, kx_, l_, g_] =
Simplify[wz4[-qz4x[kx, l, g], l, g]*E^(-I qz4x[kx, l, g] x)];
And in matrix form:
Mz[x_, kx_, l_, g_] = ({
{M1[x, kx, l, g][[1]], M2[x, kx, l, g][[1]], M3[x, kx, l, g][[1]],
M4[x, kx, l, g][[1]]},
{M1[x, kx, l, g][[2]], M2[x, kx, l, g][[2]], M3[x, kx, l, g][[2]],
M4[x, kx, l, g][[2]]},
{M1[x, kx, l, g][[3]], M2[x, kx, l, g][[3]], M3[x, kx, l, g][[3]],
M4[x, kx, l, g][[3]]},
{M1[x, kx, l, g][[4]], M2[x, kx, l, g][[4]], M3[x, kx, l, g][[4]],
M4[x, kx, l, g][[4]]}
});
Hamiltonian for region without SOI and Zeeman field:
L[x_, kx_] = 1./Sqrt[2.] ({{E^(I kx x), 0, -E^(-I kx x), 0},
{0, E^(I kx x), 0, -E^(-I kx x)},
{E^(I kx x), 0, E^(-I kx x), 0},
{0, E^(I kx x), 0, E^(-I kx x)}});
Generalising for potential barriers of width "d" every "a" units of distance
Sl[a_, kx_] = Simplify[L[2 a + d, kx].Inverse[L[a + d, kx]]];
Where it starts to really slow down
Smz[a_, d_, kx_, l_, g_] = Mz[a + d, kx, l, g].Inverse[Mz[a, kx, l, g]];
Sz[a_, d_, kx_, l_, g_] = Sl[a, kx].Smz[a, d, kx, l, g];
Impossible to compute beyond here. Even takes forever if I define every variable:
Anz[n_, a_, d_, kx_, l_, g_] := Inverse[L[((n - 1)/2) (a + d), kx]].Smz[a, d, kx, l,
g].(MatrixPower[Sz[a, d, kx, l, g], (n - 3)/2]).L[a, kx]
Tferz[n_, a_, d_, kx_, l_, g_] := ({
{1, 0, -Anz[n, a, d, kx, l, g][[1, 3]], -Anz[n, a, d, kx, l, g][[1,4]]},
{0, 1, -Anz[n, a, d, kx, l, g][[2, 3]], -Anz[n, a, d, kx, l, g][[2,4]]},
{0, 0, -Anz[n, a, d, kx, l, g][[3, 3]], -Anz[n, a, d, kx, l, g][[3,4]]},
{0, 0, -Anz[n, a, d, kx, l, g][[4, 3]], -Anz[n, a, d, kx, l, g][[4,4]]}})
Auz[n_, a_, d_, kx_, l_, g_] := Inverse[Tferz[n, a, d, kx, l, g]].{{Anz[n, a, d, kx, l, g][[1,1]]}, {Anz[n, a, d, kx, l, g][[2,1]]}, {Anz[n, a, d, kx, l, g][[3,1]]}, {Anz[n, a, d, kx, l, g][[4,1]]}}
Tuz[n_, a_, d_, kx_, l_, g_] := Abs[Auz[n, a, d, kx, l, g][[1]]]^2 + Abs[Auz[n, a, d, kx, l, g][[2]]]^2
Ruz[n_, a_, d_, kx_, l_, g_] := 1 - Tuz[n, a, d, kx, l, g]
Adz[n_, a_, d_, kx_, l_, g_] := Inverse[Tferz[n, a, d, kx, l, g]].{{Anz[n, a, d, kx, l, g][[1, 2]]}, {Anz[n, a, d, kx, l, g][[2, 2]]}, {Anz[n, a, d, kx, l, g][[3, 2]]}, {Anz[n, a, d, kx, l, g][[4, 2]]}}
Tdz[n_, a_, d_, kx_, l_, g_] := Abs[Adz[n, a, d, kx, l, g][[1]]]^2 + Abs[Adz[n, a, d, kx, l, g][[2]]]^2
Rdz[n_, a_, d_, kx_, l_, g_] := 1 - Tdz[n, a, d, kx, l, g]
Tz[n_, a_, d_, kx_, l_, g_] := (Tuz[n, a, d, kx, l, g] + Tdz[n, a, d, kx, l, g])/2.
Rz[n_, a_, d_, kx_, l_, g_] := 1 - Tz[n, a, d, kx, l, g]
Attachments: