When finding the eigenvalues and eigenfunctions of the system Hc[r, z] using NDEigensystem, the following issue arises: When changing the NDEigensystem region, the eigenvalues change periodically with a large spread. I'm only interested in the first two minimum eigenvalues (eigenvals12). Below I will try to describe the problem in more detail.

rmax={8, 9, 10, 11, 12}
z∈[-rmax,rmax]; r∈[0,rmax]

eigenvals12={-11.9385, -1.59129}

eigenvals12={-16.2499, -1.67943}

eigenvals12={-11.9385, -1.54313}

eigenvals12={-16.2443, -1.66582}

eigenvals12={-11.9385, -1.53491}

It can be seen that the eigenvalues alternate with a successive increase in the region (eigenvalues(rmax=8)≈eigenvalues(rmax=10)≈eigenvalues(rmax=12) and eigenvalues(rmax=9)≈eigenvalues(rmax=11)). \

Question: Why do the eigenvalues of the system alternate with the successive increase in the region? Especially huge changes are observed for the minimum eigenvalue. How could this be fixed?

In the code I used a cylindrical coordinate system and renamed ρ≡r:

h = 10^(-27);
eV = 1.6*10^(-12);
kB = 1.4*10^(-16);
m0 = 9.1*10^(-28);
ee = 4.8*10^(-10);
A = 10^(-8);
mv = 10^(-3)*eV;
cc = 3*10^10;
Ry0 = 2.41532928`*^-11;
e0list = 21;
ehflist = 4.7;
melist = 0.18*m0;
mhlist = 0.35*m0;
Elolist = 20*mv;
p = melist*mhlist/(melist + mhlist);
le = Sqrt[h^2/2/melist/Elolist];
lh = Sqrt[h^2/2/mhlist/Elolist];
RyX0 = Ry0*p/m0/e0list^2;
aB0 = h^2*e0list/p/ee^2;

rmax = Table[rr, {rr, 8, 12, 1}]

(*the system*)
vc[r_] = -(ee^2/(e0list*r)) - (1/ehflist - 1/e0list)*(ee^2/
      r)*((mhlist/(mhlist - melist))*
       Exp[-r/lh] - (melist/(mhlist - melist))*Exp[-r/le]);
Hc[r_, z_] := 
       Laplacian[u[r, z], {r, \[Theta], z}, "Cylindrical"]/aB0^2 + 
      vc[Sqrt[r^2 + z^2]*aB0]*u[r, z])/mv];

(*finding eigenfunctions and eigenvalues using NDEigensystem*)
energ[rr_] := 
   eigensys = 
    NDEigensystem[{Hc[r, z] + u[r, z]}, 
     u[r, z], {r, 0, rmax[[rr]]}, {z, -rmax[[rr]], rmax[[rr]]}, 100, 
     Method -> {"SpatialDiscretization" -> {"FiniteElement", \
{"MeshOptions" -> {"MaxCellMeasure" -> 0.05}}}, 
       "Eigensystem" -> {"Arnoldi", "MaxIterations" -> 10000}}];
   ord = Ordering[eigensys[[1]]];
   {eigensys[[1]], eigensys[[2]]} = {eigensys[[1]][[ord]], 
   Return[{eigensys[[1]][[1 ;; 2]] - 1, eigensys[[2]][[1 ;; 2]]}]];

valfunc = Table[energ[i], {i, 1, 5}];


(*two minimum eigenvalues, rmax={8,9,10,11,12}*)
Table[valfunc[[i, 1]], {i, 1, 5}]
  • $\begingroup$ What makes you believe the eigenvalues are not correct? After all you change the size of the region. $\endgroup$
    – user21
    Jul 18, 2023 at 4:53
  • $\begingroup$ @user21, thanks for the question! Eigenvalues should not change periodically, which behavior cannot be. They should change evenly, if the region is chosen with an inappropriate size, or not change at all in the case of any region that is suitable in size (because the eigenfunctions in this case are localized near zero). $\endgroup$
    – Mam Mam
    Jul 18, 2023 at 8:26
  • $\begingroup$ The way this is now, it's too complicated to debug. Simplify this to a single case (B=0), remove the ordering, the output you have in text form is much better than the many outputs you have in the code. Remove the ParallelTable etc. $\endgroup$
    – user21
    Jul 18, 2023 at 8:56
  • $\begingroup$ As a side note: reference.wolfram.com/language/ref/… $\endgroup$
    – user21
    Jul 18, 2023 at 8:56
  • $\begingroup$ I would suggest using the natural units of your system instead of the physical units. In this way, you might avoid numerical errors. $\endgroup$
    – Jam
    Feb 13 at 15:59

1 Answer 1


Probably numerically form of level crossing, a common phenomenon e.g. in the optical phonon spectrum in solids where the spectrum depends nearly continously linear on the wave vector with different slopes for different oscillation modes of lattice content.

The eigenvalues scale somehow as rational functions in the different cordinate directions, linear on the circle, quadratic in the box, fourth order in Coulomb systems. If one fine-tunes the crossings they evolve into neraly touching hyperbolas. Its a theorem in perturbation theory, that true linear level crossing is impossible.

Variation of the spectrum of a model Hamiltonian with respect to variation of parameters was a field of theoretical interest in the 1950-70ties when no numerical methods existed to control interaction of variation of model data and their measurable spectral consequences.

  • $\begingroup$ Thanks a lot! But the question is purely mathematical, I am just trying to solve a problem on eigenvalues with a certain potential. So my question is related to the code, why Mathematica gives such different eigenvalues that repeat when the region changes. How to understand which of them is reliable and how to make the program work correctly, so that the obtained eigenvalues would be unambiguously defined in the whole region? $\endgroup$
    – Mam Mam
    Jul 15, 2023 at 17:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.