Why does the minimum eigenvalue change dramatically when one basis function is added to the basis set?

Question

I have a basis set that describes with high accuracy the ground state of this system:$$H=-\frac{1}{2}\Delta-\frac{1}{r}+\frac{25}{8}\rho^2-5/2$$where $r=(\rho,z,\phi)$ is a coordinate in the cylindrical system.

Ground state from the literature $E_{\min}=-1.3803975$

The anisotropic Gaussian basis set which I use to solve this task:$$\psi_j=e^{-b_{0j}\;z^2}e^{-a_{0j}\;\rho^2}$$where $a_{0j}$ and $b_{0j}$ are parameters

$a_{0j}$={1.25, 1.25003, 1.25006, 1.251393, 1.252726, 1.2587955, 1.264865, 1.260499, 1.303186, 1.546232, 2.337185, 4.507307, 9.426298, 20.270036, 43.977048, 95.675215, 208.32642, 453.736619, 988.32236, 2152.803005, 4689.357171, 10214.647096}

$b_{0j}$={0.026284, 0.0417685, 0.057253, 0.09098300000000001, 0.124713,0.19818550000000001, 0.271658,0.417292, 0.61588,1.008925, 1.949878, 4.24735, 9.251849, 20.15297, 43.898488,95.622497, 208.291042, 453.712878`, 988.306428, 2152.792314,4689.349997, 10214.642282}

This basis gives the energy of the ground state $E_{\min}=-1.3803971$, which agrees perfectly with the value from the literature.

Now if we add only one element to the list $a_{0j}$ and list $b_{0j}$ (1.2626819999999999 and 0.344475 respectively), then the absolute value of minimum energy becomes huge ($E_{\min}=-7.1824042$).
Why is this happening?

It is known from theory that an increase in the number of basis functions should either increase the accuracy or not change the values of the eigenvalues.

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

ClearAll["Global`*"]
b00 = {0.026284`, 0.0417685`, 0.057253`, 0.09098300000000001`, 
   0.124713`, 0.19818550000000001`, 0.271658`,(*0.344475`,*)0.417292`,
    0.61588`, 1.008925`, 1.949878`, 4.24735`, 9.251849`, 20.15297`, 
   43.898488`, 95.622497`, 208.291042`, 453.712878`, 988.306428`, 
   2152.792314`, 4689.349997`, 10214.642282`};
a00 = {1.25`, 1.25003`, 1.25006`, 1.251393`, 1.252726`, 1.2587955`, 
   1.264865`,(*1.2626819999999999`,*)1.260499`, 1.303186`, 1.546232`, 
   2.337185`, 4.507307`, 9.426298`, 20.270036`, 43.977048`, 
   95.675215`, 208.32642`, 453.736619`, 988.32236`, 2152.803005`, 
   4689.357171`, 10214.647096`};

nmax = Dimensions[b00][[1]];

a0 [j_] := a00[[j]];
b0 [j_] := b00[[j]];

Psi[r_, z_, j_] := Exp[-b0[j]*z^2]*Exp[-a0[j]*r^2];


Kk[r_, z_, j1_, j2_] := 
  FullSimplify[
   Psi[r, z, j2]*
    Laplacian[Psi[r, z, j1], {r, \[Theta], z}, "Cylindrical"]*r*2*Pi];
Kx[j1_, j2_] := -(1/2)*
   NIntegrate[
    Kk[r, z, j1, j2], {r, 0, Infinity}, {z, -Infinity, Infinity}];
Kxx = Table[Kx[j1, j2], {j1, 1, nmax}, {j2, 1, nmax}];


Ka[r_, z_, j1_, j2_] := 
  FullSimplify[Psi[r, z, j2]*25/8*r^2*Psi[r, z, j1]*r*2*Pi];
KA[j1_, j2_] := 
  KA[j1, j2] = 
   KA[j2, j1] = 
    NIntegrate[
     Ka[r, z, j1, j2], {r, 0, Infinity}, {z, -Infinity, Infinity}];
KAx = Table[KA[j1, j2], {j1, 1, nmax}, {j2, 1, nmax}];


Ks[r_, z_, j1_, j2_] := 
  FullSimplify[Psi[r, z, j2]*(-5/2)*Psi[r, z, j1]*r*2*Pi];
KS[j1_, j2_] := 
  KS[j1, j2] = 
   KS[j2, j1] = 
    NIntegrate[
     Ks[r, z, j1, j2], {r, 0, Infinity}, {z, -Infinity, Infinity}];
KSx = Table[KS[j1, j2], {j1, 1, nmax}, {j2, 1, nmax}];

VP1[r_, z_] := -(1/Sqrt[r^2 + z^2]);
Px1[j1_, j2_] := 
  Px1[j1, j2] = 
   Px1[j2, j1] = 
    NIntegrate[
     Psi[r, z, j2]*VP1[r, z]*Psi[r, z, j1]*r*2*Pi, {r, 0, 
      Infinity}, {z, -Infinity, Infinity}];
Pxx1 = Table[Px1[j1, j2], {j1, 1, nmax}, {j2, 1, nmax}];

Bb[j1_, j2_] := 
  Bb[j1, j2] = 
   Bb[j2, j1] = 
    NIntegrate[
     Psi[r, z, j2]*Psi[r, z, j1]*r*2*Pi, {r, 0, 
      Infinity}, {z, -Infinity, Infinity}];
Bxx = Table[Bb[j1, j2], {j1, 1, nmax}, {j2, 1, nmax}];

EE1 = Eigenvalues[{Kxx + Pxx1 + KAx + KSx, Bxx}];
Min[EE1]

Out[48]= -1.3804

Answer 1

We use code from our answer here. With a small modifications we have

ClearAll["Global`*"]

a0j = {1.25, 1.25003, 1.25006, 1.251393, 1.252726, 1.2587955, 
   1.264865, 1.260499, 1.303186, 1.546232, 2.337185, 4.507307, 
   9.426298, 20.270036, 43.977048, 95.675215, 208.32642, 453.736619, 
   988.32236, 2152.803005, 4689.357171, 10214.647096};

b0j = {0.026284, 0.0417685, 0.057253, 0.09098300000000001, 0.124713, 
   0.19818550000000001, 0.271658, 0.417292, 0.61588, 1.008925, 
   1.949878, 4.24735, 9.251849, 20.15297, 43.898488, 95.622497, 
   208.291042, 453.712878`, 988.306428, 2152.792314, 4689.349997, 
   10214.642282};
nmax = Length[a0j]; var = 
 Join[Table[a[n], {n, nmax}], Table[b[n], {n, nmax}], 
  Table[c[n], {n, nmax}]];
Psi[r_, z_, j_] := Exp[-b[j]*z^2]*Exp[-a[j]*r^2];



(*kinetic energy*)
Kk[r_, z_, n1_, n2_] := 
  FullSimplify[
   Psi[r, z, n2]*
    Laplacian[Psi[r, z, n1], {r, \[Theta], z}, "Cylindrical"]];

ss = Integrate[
  Kk[r, z, l1, l2] r, {r, 0, \[Infinity]}, {z, -Infinity, Infinity}, 
  Assumptions -> {a[l1] > 0, b[l1] > 0, a[l2] > 0, b[l2] > 0}];

Kx = -1/2 2 Pi Sum[c[l1] c[l2] ss, {l1, nmax}, {l2, nmax}];

(*potential energy*)
VP1[r_, z_] := -1/Sqrt[r^2 + z^2] + 25/8 r^2 - 5/2;
Px = Integrate[
  Psi[r, z, l2]*VP1[r, z]*Psi[r, z, l1]*r , {r, 
   0, \[Infinity]}, {z, -Infinity, Infinity}, 
  Assumptions -> {a[l1] > 0, b[l1] > 0, a[l2] > 0, b[l2] > 0}]

Px = 2 Pi Sum[c[l1] c[l2] Px, {l1, nmax}, {l2, nmax}];



int = Integrate[
  Psi[r, z, l1] Psi[r, z, l2] r , {r, 0, \[Infinity]}, {z, -Infinity, 
   Infinity}, 
  Assumptions -> {a[l1] > 0, b[l1] > 0, a[l2] > 0, b[l2] > 0}];

norm = {2 Pi Sum[c[l1] c[l2] int, {l1, nmax}, {l2, nmax}] == 
   1}; cons = 
 Join[Table[a[n] == a0j[[n]], {n, nmax}], 
  Table[b[n] == b0j[[n]], {n, nmax}], norm];

With using Kx, Px we can minimize energy as follows

sol = 
 NMinimize[{Kx + Px, cons}, var, Method -> "DifferentialEvolution"]

Out[]= {-1.3804 + 0. I, {a[1] -> 1.25, a[2] -> 1.25003, 
  a[3] -> 1.25006, a[4] -> 1.25139, a[5] -> 1.25273, a[6] -> 1.2588, 
  a[7] -> 1.26487, a[8] -> 1.2605, a[9] -> 1.30319, a[10] -> 1.54623, 
  a[11] -> 2.33719, a[12] -> 4.50731, a[13] -> 9.4263, a[14] -> 20.27,
   a[15] -> 43.977, a[16] -> 95.6752, a[17] -> 208.326, 
  a[18] -> 453.737, a[19] -> 988.322, a[20] -> 2152.8, 
  a[21] -> 4689.36, a[22] -> 10214.6, b[1] -> 0.026284, 
  b[2] -> 0.0417685, b[3] -> 0.057253, b[4] -> 0.090983, 
  b[5] -> 0.124713, b[6] -> 0.198186, b[7] -> 0.271658, 
  b[8] -> 0.417292, b[9] -> 0.61588, b[10] -> 1.00893, 
  b[11] -> 1.94988, b[12] -> 4.24735, b[13] -> 9.25185, 
  b[14] -> 20.153, b[15] -> 43.8985, b[16] -> 95.6225, 
  b[17] -> 208.291, b[18] -> 453.713, b[19] -> 988.306, 
  b[20] -> 2152.79, b[21] -> 4689.35, b[22] -> 10214.6, 
  c[1] -> 0.0065233 + 0. I, c[2] -> -0.0605427 + 0. I, 
  c[3] -> 0.146483 + 0. I, c[4] -> -0.315022 + 0. I, 
  c[5] -> 0.416621 + 0. I, c[6] -> -0.360113 + 0. I, 
  c[7] -> 0.391223 + 0. I, c[8] -> 0.0120385 + 0. I, 
  c[9] -> 0.168188 + 0. I, c[10] -> 0.177786 + 0. I, 
  c[11] -> 0.173833 + 0. I, c[12] -> 0.124923 + 0. I, 
  c[13] -> 0.0822642 + 0. I, c[14] -> 0.0567521 + 0. I, 
  c[15] -> 0.0373556 + 0. I, c[16] -> 0.0261548 + 0. I, 
  c[17] -> 0.0166696 + 0. I, c[18] -> 0.0127193 + 0. I, 
  c[19] -> 0.00636482 + 0. I, c[20] -> 0.00797396 + 0. I, 
  c[21] -> -0.000160254 + 0. I, c[22] -> 0.00659502 + 0. I}}

Therefore, our answer for the ground state energy is sol[[1]]=-1.3803982710899914 while from the literature (?) we have $E_{\min}=-1.3803975$. If we increase bases by adding 1.2626819999999999 and 0.344475 to a0j,b0j respectively, then we have $E_{min}=-1.380398305740276$.

Answer 2

In the standard approach to find the spectrum of an operator

$\left\{ H=-\Delta_{r,r} + V\left(r \right), d\mu= r d r, r>0\right\}$

on the unit sphere in the Hilbert space

$\mathit H = \left\{ r \to f\left(r\right) \in \mathit C | \int_0^\infty r dr \left(f f^*\right)\left(r\right) =1 \right\}$

one takes a finite subspace on the countable infinite dimensional unit sphere, a n dimensional unit subsphere, like any circle on the 3-sphere. Since we don't know the operator very well at r=0 its a safe way to take the series of spaces of eigenvectors of the harmonic oscillator as a better orthogonal choice compared to the pure set of powers $e^{-(r^2/2 )}\ r^n$ .

In order to get rid of the first order term $\partial_r$ and the measure r, one may put a factor $e^{-r^/4} \ r^{-1/2} $ in front of the functions

$ \mathit H_n\left(c\right) = \sum_{k=0}^n c_k \mathit H_k[r]e^{-r^2/2}$

This set does not form an orthnormal system in the Hilbert space $(R_+,r dr)$ , but they have the advantage to be algebraic eigenfunctions of $-\partial_r^2 +r^2$ with eigenvalues 2n+1

((-D[#, {r, 2}] + r^2 #)/# &)[(HermiteH[5,r] Exp[-r^2/2])] // Simplify
 11

So the energy of the ocillator part is simply $ \sum_k c_k^2 (2k+1)$.

The adapted minimum problem is now

$\min_{c,|c|=1} (\sum_k (2k+1) c_k^2 - \int_0^\infty dr \left(\sum_k c_k f_k\right)^2 )$

The ground state is the limit of dimension $n\to \infty$. The minimum becomes smaller by each step, because the minimum cannot grow by enlarging the search space. No basis function must be left out, because the ground state is orthogonal to the complete rest of the unit sphere. Therefore the series of n-spheres used must be forming dense subsets in the unit sphere. This is the crux with infinite spaces, when the first n appoximations, by chance, are orthogonal to the ground state.

Spectral theory in practise one determines the ground state in an n-dimensional subsphere and takes the projection to the orthogonal complement to find the ground state in that subspace, the first exited state.

Now: The Hamilton operator given is not bounded from below, because the potential $-1/r$ in cylinder coordinates

$H = -\Delta -1/r$

This potential has a charge density

$ - \Delta 1/r =- 1/r \partial_r r \partial 1/r = - 1/r^3 $

$\int_0 r dr 1/r^3 = \infty $

has an infinite charge per z-length.

The potential energy of a wave at r=0 may be of form $r^{-1+eps}$, such that $\int\ dr\ r\ r^{-2+\varepsilon}$ is integrable, but the energy with one extra power $r^{-1}$ is infinite.

For spherical problems the measure is $dr\ r^2$, so $V =-1/r$ yields a finite negative potential energy for all square integrable functions.

If the Hamiltonian is not bounded from below at a singularity, physicist begin to speak about absorption at the center.