Why does NDEigensystem not show the minimum eigenvalue for a certain parameter range in the cylindrical coordinate system?

In my previous question Why NDEigensystem does not show the minimum eigenvalue?, I asked why the NDEigensystem does not show the minimum eigenvalue for the following system:
$$H[u(\rho,z)]=-\frac{1}{2}\Delta u(\rho,z)+Vu(\rho,z)$$, where $$V=-\frac{1}{\sqrt{\rho^2+z^2}}-\frac{3e^{-3\sqrt{\rho^2+z^2}}+e^{-7\sqrt{\rho^2+z^2}}}{\sqrt{\rho^2+z^2}}$$

In order for NDEigensystem to find the minimum eigenvalue, as it turned out in this case, it helps to reformulate the problem from a cylindrical coordinate system to a spherical one. Indeed, it is clear that the problem has spherical symmetry.

Here I would like to consider the following system:
$$H[u(\rho,z)]=-\frac{1}{2}\Delta u(\rho,z)+Vu(\rho,z)$$, where $$V=-\frac{1}{\sqrt{\rho^2+z^2}}-\frac{3e^{-3\sqrt{\rho^2+z^2}}+e^{-7\sqrt{\rho^2+z^2}}}{\sqrt{\rho^2+z^2}}+eps*\rho^2$$

This system already has cylindrical symmetry. NDEigensystem can find the minimum eigenvalue only for sufficiently large values of eps (eps > 1/10) while for eps < 1/10 NDEigensystem returns eigenvalues starting from the first excited one.

eps = 2/10, eigenvalues:
{-4.79274, 0.252694, 0.331186, 0.516154, 0.534339, 0.578607, 0.584517, 0.599544, 0.603317, 0.626888,...

the minimum eigenvalue returned (eps = 2/10) = -4.79274

eps = 8/100, eigenvalues:
{0.0176122, 0.133175, 0.28278, 0.307919, 0.345772, 0.353845, 0.366874, 0.372319, 0.393895, 0.404982,...

the minimum eigenvalue returned (eps = 8/100) = 0.0176122

How to get the minimum eigenvalues for eps < 1/10 using NDEigensystem?

In the code I renamed $$\rho≡r$$

ClearAll["Global*"]
rmax = 30;
zmax = 30;
eps = 8/100;
V[r_, z_] := -(1/Sqrt[r^2 + z^2]) - (
3 E^(-3 Sqrt[r^2 + z^2]) + E^(-7 Sqrt[r^2 + z^2]))/Sqrt[r^2 + z^2] +
eps*r^2
H = Simplify[-1/2*
Laplacian[u[r, z], {r, \[Theta], z}, "Cylindrical"] +
V[r, z]*u[r, z]];
{vals, funs} =
NDEigensystem[{H + u[r, z]},
u[r, z], {r, 0, rmax}, {z, -zmax, zmax}, 100,
Method -> {"SpatialDiscretization" -> {"FiniteElement",{"MeshOptions" -> {"MaxCellMeasure" -> 0.05}}},
"Eigensystem" -> {"Arnoldi", "MaxIterations" -> 10000}}];
Sort[vals - 1]

(*{0.0176122, 0.133175, 0.28278, 0.307919, 0.345772, 0.353845, \
0.366874, 0.372319, 0.393895, 0.404982, 0.437912, 0.453121, 0.494824, \
0.513963, 0.563591, 0.586744, 0.643667, 0.671056, 0.734496, 0.766569, \
0.834287, 0.87273, 0.914979, 0.964405, 0.98457, 1.02342, 1.0778, \
1.10505, 1.1153, 1.13243, 1.1516, 1.15974, 1.17009, 1.17874, 1.19899, \
1.21414, 1.21703, 1.24727, 1.25582, 1.27378, 1.30444, 1.32846, \
1.36032, 1.37677, 1.39599, 1.42214, 1.45572, 1.48835, 1.51769, \
1.54832, 1.5687, 1.59619, 1.64993, 1.67898, 1.69014, 1.75062, 1.7534, \
1.78336, 1.81212, 1.84621, 1.87005, 1.89194, 1.91728, 1.92068, \
1.93658, 1.95466, 1.95792, 1.96202, 1.97229, 1.98234, 2.00268, \
2.01963, 2.03, 2.04612, 2.06667, 2.06758, 2.08391, 2.11124, 2.13287, \
2.14269, 2.17302, 2.18454, 2.20754, 2.23827, 2.25266, 2.27148, \
2.29537, 2.33371, 2.34156, 2.35559, 2.38996, 2.40919, 2.45431, \
2.47053, 2.49862, 2.4993, 2.55017, 2.56405, 2.5894, 2.59755}*)

• Oh! It is very interesting, that for eps>1/10 we can compute ground state using the Arnoldi algorithm. For eps=0 we return to the previous question discussed on mathematica.stackexchange.com/questions/286303/… Commented Jun 14, 2023 at 3:20
• @Alex Trounev, thanks! Yes, it's strange also because 1/10 is not such a small value to neglect (or not?) Commented Jun 14, 2023 at 6:48
• Box {r, 0, rmax}, {z, -zmax, zmax} with rmax = 30; zmax = 30; is too large for the ground state. Try rmax = 12; zmax = 10; eps = 1/100; and see how ground state looks like. Commented Jun 14, 2023 at 15:03
• @Alex Trounev, thanks a lot! Yes, for these sizes rmax = 12; zmax = 10, NDEigensystem can get the ground state. I have plotted a few first eigenfunctions, but the features are not very clear. Could you explain please Commented Jun 14, 2023 at 22:54

Using small box we can compute ground state even for eps=10^-3 as follows

ClearAll["Global*"]
rmax = 12;
zmax = 10;
eps = 1/1000;
V[r_, z_] := -(1/Sqrt[r^2 + z^2]) - (3 E^(-3 Sqrt[r^2 + z^2]) +
E^(-7 Sqrt[r^2 + z^2]))/Sqrt[r^2 + z^2] + eps*r^2
H = Simplify[-1/2*
Laplacian[u[r, z], {r, \[Theta], z}, "Cylindrical"] +
V[r, z]*u[r, z]];
{vals, funs} =
NDEigensystem[{H + u[r, z]},
u[r, z], {r, 0, rmax}, {z, -zmax, zmax}, 100,
Method -> {"SpatialDiscretization" -> {"FiniteElement", \
{"MeshOptions" -> {"MaxCellMeasure" -> 0.05}}},
"Eigensystem" -> {"Arnoldi", "MaxIterations" -> 10000}}];

Sort[vals - 1] // First

(*Out[]= -3.58121*)


The problem is that ground state eigenfunction localized around origin r=0, z=0 so that in a large box we can't compute it with a given mesh. To see this feature let us visualize this function as

vf = SortBy[Table[{vals[[i]], funs[[i]]}, {i, Length[vals]}], First];

Table[
Plot3D[vf[[i]][[2]], {r, 0, rmax}, {z, -zmax, zmax},
PlotTheme -> "Minimal", PlotRange -> All, BoxRatios -> Automatic,
Boxed -> False, Axes -> False, PlotLabel -> vf[[i]][[1]] - 1,
ColorFunction -> Hue], {i, 6}]


• Thank you very much! Do I understand correctly that the mesh is adjusted here {"MaxCellMeasure" -> 0.05}? Commented Jun 15, 2023 at 12:55
• @MamMam Yes, we can play with mesh, but there is also restriction on this way. Commented Jun 15, 2023 at 13:23
• Thanks! Could you please give an example of restrictions? Commented Jun 15, 2023 at 13:45
• @MamMam I can recommend to use exact solution to compare with. For example the ground state for Coulomb is -1/2. Try to reproduce this value with different mesh. Also very supportive the parametric research as in my answer on mathematica.stackexchange.com/questions/286303/… Commented Jun 16, 2023 at 1:53
• @MamMam You are welcome! Commented Jun 16, 2023 at 7:30