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

Question

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}*)

Answer 1

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}]