# Issue with numerical evaluation of eigensystem

I am calculating eigenvalues of a Hamiltonian numerically but I am getting avoided crossings and gap in the curves (see the output of the code) which are not correct. Please help me out to resolve this issue in the attached notebook. I tried a lot to resolve it but unsuccessful. Thanks.

c = 3*10^8;
m = (0.16*5.1*10^5)/(c^2);
\[HBar] = 6.6*10^-16;
v = 5*10^5;
\[Lambda] = 0.5*1.65*10^-28;
g = 8.4;
mu = 5.78*10^-5;
\[CapitalDelta] = -0.5*g*mu*B;
\[Omega]0 = (\[HBar]*B)/(m);
\[Omega]D = Sqrt[2*\[HBar]*B*v^2];
ll = Sqrt[(2*B)/(\[HBar])];
Block[{L = 100, \[CapitalGamma] = 1.3*10^-3, i, j, one, a, adag, wa,
qd, H, Es, list, nn, x},
one = IdentityMatrix[L, SparseArray,
WorkingPrecision -> MachinePrecision];
a = SparseArray[{j_, i_} /; i - j == 1 -> Sqrt[j], {L, L}, 0.];
adag = SparseArray[{j_, i_} /; j - i == 1 -> Sqrt[i], {L, L}, 0.];
H[B_] =
ArrayFlatten[
Normal[{{(\[CapitalDelta] + 0.5*\[Omega]0)*one + qd + wa,
I*\[Omega]D*a}, {-I*\[Omega]D*
adag, (-\[CapitalDelta] + 0.5*\[Omega]0)*one + qd - wa}}]];
list[B_?NumericQ] := Sort@Eigenvalues[H[N[B]]];
list[B_?NumericQ, nn_] := list[B][[nn]];
Plot[Evaluate@
Transpose@Delete[Table[list[B, nn], {nn, 2*L}], L + 1], {B, 0, 40},
PlotRange -> {{0, 40}, {-0.1, 0.1}}, ImageSize -> 500]
]


• Does this answer your question? Tracking Eigenvalues Through a Crossing Jun 21, 2023 at 20:07
• Thanks. But this is not answer to my question. I want to solve it with Eigensystem. Is there any option to correct this code If you know please let me know? Jun 21, 2023 at 20:30
• But in the code you provided, you use Eigenvalues not Eigensystem ... Jun 21, 2023 at 20:48
• Yes, you are right. I want the correct solution with Eigenvalues. If you can do it with Eigensystem, it will also appreciated. I tried it with both Eigenvalues and Eigensystem but all my efforts were in vain. Therefore, I posted it here. Jun 21, 2023 at 20:55
• Crossposted here. Jun 24, 2023 at 23:14

Your code has fatal error in simultaneously usage Plot and Transpose. Also code is time consuming since list[B] computes 2 L times for every B while we need one time only. To simplifier problem we first compute one eigenvalue plot as follows

pl1=Block[{L = 100, \[CapitalGamma] = 1.3*10^-3, i, j, one, a, adag, wa,
qd, H, Es, list, nn, x},
one = IdentityMatrix[L, SparseArray,
WorkingPrecision -> MachinePrecision];
a = SparseArray[{j_, i_} /; i - j == 1 -> Sqrt[j], {L, L}, 0.];
adag = SparseArray[{j_, i_} /; j - i == 1 -> Sqrt[i], {L, L}, 0.];
H[B_] =
ArrayFlatten[
Normal[{{(\[CapitalDelta] + 0.5*\[Omega]0)*one + qd + wa,
I*\[Omega]D*a}, {-I*\[Omega]D*
adag, (-\[CapitalDelta] + 0.5*\[Omega]0)*one + qd - wa}}]];
list[B_?NumericQ] := Sort@Eigenvalues[H[N[B]]];
list[B_?NumericQ, nn_] := list[B][[nn]];
Plot[Evaluate@list[B, 2 L], {B, 0, 40}, PlotRange -> All,
ImageSize -> 500]]


Then we compute data for the last 2L-1 eigenvalues as

eb = Block[{L = 100, \[CapitalGamma] = 1.3*10^-3, i, j, one, a, adag,
wa, qd, H, Es, list, nn, x},
one = IdentityMatrix[L, SparseArray,
WorkingPrecision -> MachinePrecision];
a = SparseArray[{j_, i_} /; i - j == 1 -> Sqrt[j], {L, L}, 0.];
SparseArray[{j_, i_} /; j - i == 1 -> Sqrt[i], {L, L}, 0.];
H[B_] =
ArrayFlatten[
Normal[{{(\[CapitalDelta] + 0.5*\[Omega]0)*one + qd + wa,
I*\[Omega]D*a}, {-I*\[Omega]D*
adag, (-\[CapitalDelta] + 0.5*\[Omega]0)*one + qd - wa}}]];
list[B_?NumericQ] := Sort@Eigenvalues[H[N[B]]];
Table[Drop[list[B], L + 1], {B, 0, 40}]]; // AbsoluteTiming


Note, that it takes 0.24 s only. Finally we show pl1 and eb in one plot

Show[pl1,
ListPlot[eb // Transpose, PlotRange -> All, DataRange -> {0, 40},
PlotStyle -> {PointSize[Medium]}]]


As we can see there is no any crossing or gap. Also we can check region $$0\le B\le 1$$

pl2=Block[{L = 100, \[CapitalGamma] = 1.3*10^-3, i, j, one, a, adag, wa,
qd, H, Es, list, nn, x},
one = IdentityMatrix[L, SparseArray,
WorkingPrecision -> MachinePrecision];
a = SparseArray[{j_, i_} /; i - j == 1 -> Sqrt[j], {L, L}, 0.];
adag = SparseArray[{j_, i_} /; j - i == 1 -> Sqrt[i], {L, L}, 0.];
H[B_] =
ArrayFlatten[
Normal[{{(\[CapitalDelta] + 0.5*\[Omega]0)*one + qd + wa,
I*\[Omega]D*a}, {-I*\[Omega]D*
adag, (-\[CapitalDelta] + 0.5*\[Omega]0)*one + qd - wa}}]];
list[B_?NumericQ] := Sort@Eigenvalues[H[N[B]]];
list[B_?NumericQ, nn_] := list[B][[nn]];
Plot[Evaluate@list[B, 2 L], {B, 0, 1}, ImageSize -> 500,
PlotRange -> {-.1, .3}]]
vals = Block[{L = 100, \[CapitalGamma] = 1.3*10^-3, i, j, one, a,
adag, wa, qd, H, Es, list, nn, x},
one = IdentityMatrix[L, SparseArray,
WorkingPrecision -> MachinePrecision];
a = SparseArray[{j_, i_} /; i - j == 1 -> Sqrt[j], {L, L}, 0.];
SparseArray[{j_, i_} /; j - i == 1 -> Sqrt[i], {L, L}, 0.];
H[B_] =
ArrayFlatten[
Normal[{{(\[CapitalDelta] + 0.5*\[Omega]0)*one + qd + wa,
I*\[Omega]D*a}, {-I*\[Omega]D*
adag, (-\[CapitalDelta] + 0.5*\[Omega]0)*one + qd - wa}}]];
list[B_?NumericQ] := Sort@Eigenvalues[H[N[B]]];
Table[Drop[list[B], L + 1], {B, 0, 1, .05}]]; // AbsoluteTiming


Visualization

Show[pl2,
ListPlot[vals // Transpose, PlotRange -> All, DataRange -> {0, 1},
PlotStyle -> {PointSize[Medium]}]]


In this region also there is no crossing or gap.

Update 1. Nevertheless these data consist of zigzag due to Arnoldi algorithm instability. To exclude zigzag for first 100 eigenvalues, we can increase dimension of matrix up to $$400\times 400$$ and then take 100 eigenvalues only as follows

c = 3*10^8;
m = (0.16*5.1*10^5)/(c^2);
\[HBar] = 6.6*10^-16;
v = 5*10^5;
\[Lambda] = 0.5*1.65*10^-28;
g = 8.4;
mu = 5.78*10^-5;
\[CapitalDelta] = -0.5*g*mu*B;
\[Omega]0 = (\[HBar]*B)/(m);
\[Omega]D = Sqrt[2*\[HBar]*B*v^2];
ll = Sqrt[(2*B)/(\[HBar])];
eb = Block[{L = 200, \[CapitalGamma] = 1.3*10^-3, i, j, one, a, adag,
wa, qd, H, Es, list, nn, x},
one = IdentityMatrix[L, SparseArray,
WorkingPrecision -> MachinePrecision];
a = SparseArray[{j_, i_} /; i - j == 1 -> Sqrt[j], {L, L}, 0.];
SparseArray[{j_, i_} /; j - i == 1 -> Sqrt[i], {L, L}, 0.];
H[B_] =
ArrayFlatten[
Normal[{{(\[CapitalDelta] + 0.5*\[Omega]0)*one + qd + wa,
I*\[Omega]D*a}, {-I*\[Omega]D*
adag, (-\[CapitalDelta] + 0.5*\[Omega]0)*one + qd - wa}}]];
list[B_?NumericQ] := Take[Sort@Eigenvalues[H[N[B]]], 100];
Table[list[B], {B, 0, 40, 1/4}]];


Visualization

eb1 = eb // Transpose;
nn = Length[eb1]
tt = Table[
ListLinePlot[eb1[[i]], PlotRange -> All, DataRange -> {0, 40},
ColorFunction -> Hue], {i, nn}];
{Show[Take[tt, 10]], Show[Take[tt, -10]]}


• @ Alex Trounev thanks. But there is still zigzag like pattern between curves if you plot by using ListLinePlot command. How can it be removed this trend? Jun 22, 2023 at 18:41
• @user199 The reason for zigzag is that you analyze very exited states computed with Arnoldi algorithm. Actually only ground state + maybe first 20 eigenvalues can be computed with high accuracy using Arnoldi, While all other are just a garbage. Also please pay attention that Arnoldi algorithm is unstable, and it is why there is shift effect discussed on mathematica.stackexchange.com/questions/286303/… Jun 23, 2023 at 5:22
• I have made a lot of tries to remove this zigzag trend but I could not. If you have a more specific solution for it, please share. I have no idea about it. Also, I tried to understand it by using the above link but I could not remove them again. Jun 23, 2023 at 6:29
• @user199 Tell me please, how many states beginning from the ground state do you want to analyze? Jun 23, 2023 at 8:28
• I need 100 states for analyze. Jun 23, 2023 at 17:46