Why does using different basis functions in the matrix method give different ground state energies?

Question

I would like to find the ground state energy of the following Hamiltonian:

$$H=-\frac{1}{2}\Delta -\frac{1}{r}-\frac{1}{2r}e^{-1.5r}$$where $-\frac{1}{r}-\frac{1}{2r}e^{-1.5r}$ is a Potential energy.

In order to do it I have chosen matrix method and have done it with two basis functions.

1. The first basis functions which I used were Coulomb eigenfunctions (l=0 for ground state):

$$\text{Psi}(\text{r},\text{n})=\frac {2} {n^2} \sqrt{\frac{n!}{(n-1)!}} e^{-\frac{r}{n}} \, _1F_1\left(1-n;2;\frac{2 r}{n}\right)$$ where $_1F_1\left(1-n;2;\frac{2 r}{n}\right)$ is a Confluent hypergeometric function, $n = 1, 2, 3, 4 ,5...$

After constructed the matrix and finding eigenvalues I got the following energy minimum:

Emin = -0.670837 $\quad$ (using 6 basis functions)
Emin = -0.671109 $\quad$ (using 11 basis functions)
Emin = -0.671174 $\quad$ (using 16 basis functions) \

We can see that an increase in the number of the basis functions doesn't change the energy much, therefore the ground state energy is $\approx- 0.67$

2. Now move to the second basis functions which I used, these were 3D oscillator eigenfunctions (again l=0 for ground state): $$\text{Psi}(\text{r},\text{n})\text{=}(-1)^n e^ {-\frac{r^2}{2} }\sqrt{\frac{2 n!}{\Gamma \left(n+\frac{3}{2}\right)}} L_n^{\frac{1}{2}}\left(r^2\right)$$ where $L_n^{\frac{1}{2}}\left(r^2\right)$ is a Laguerre polynomial, $n = 0, 1, 2, 3, 4 ...$

Again after constructed the matrix and finding eigenvalues I got the following energy minimum:

Emin = -0.675029 $\quad$ (using 6 basis functions)
Emin = -0.697315 $\quad$ (using 11 basis functions)
Emin = -0.705009 $\quad$ (using 16 basis functions)
Emin = -0.709063 $\quad$ (using 21 basis functions) \

From this we can understand that the ground state energy $\approx -0.71$

Why the ground state energies received from different basis functions in 1. and 2. are different? While they must be the same.

Below you can see codes in Wolfram Mathematica, with the help of which I received the energies values.

In points (I) and (III) I have put the codes that are clear and simple, but these codes are calculated slowly. That's why in points (II) and (IV) I have put upgraded code that are calculated fast.
(I) and (II) are the codes based on Coulomb basis functions, (III) and (IV) are the codes based on 3D oscillator basis functions.

(I)

In[317]:= ClearAll["Global`*"]

nmax = 16;

(*Coulomb basis l=0*)
Psi[r_, n_] = (2 E^(-(r/n)) Sqrt[n!/(-1 + n)!] Hypergeometric1F1[
      1 - n, 2, (2 r)/n])/n^2;
(*kinetic energy*)
Kk[r_, n1_, n2_] := 
  FullSimplify[
   Psi[r, n2]*
    Laplacian[Psi[r, n1], {r, \[Theta], \[Phi]}, "Spherical"]];
Kx[n1_, n2_] := -1/2*Integrate[Kk[r, n1, n2]*r^2, {r, 0, \[Infinity]}];
KK = Table[Kx[n1, n2], {n1, 1, nmax}, {n2, 1, nmax}];

(*potential energy*)
VH1[r_] := -1/r;
VH2[r_] := -(1/(2*r))*Exp[-r*1.5];
Px1[n1_, n2_] := 
  Integrate[Psi[r, n2]*VH1[r]*Psi[r, n1]*r^2, {r, 0, \[Infinity]}];
Px2[n1_, n2_] := 
  Integrate[
   Psi[r, n2]*FullSimplify[VH2[r]]*Psi[r, n1]*r^2, {r, 
    0, \[Infinity]}];
Px[n1_, n2_] := Px1[n1, n2] + Px2[n1, n2];

PP = Table[Px[n1, n2], {n1, 1, nmax}, {n2, 1, nmax}];
EE = Min[Eigenvalues[KK + PP]]

Out[329]= -0.671174

(II)

In[250]:= ClearAll["Global`*"]
nmax = 16;

(*Coulomb basis l=0*)
(*Psi[r_,n_]=(2 E^(-(r/n)) Sqrt[n!/(-1+n)!] \
Hypergeometric1F1[1-n,2,(2 r)/n])/n^2;*)

(*kinetic energy*)

ffk[a_, n_] = 
  Integrate[r^(n + 1)*Exp[-a*r], {r, 0, Infinity}, 
   Assumptions -> {a > 0, n > 0}];
coeffk22[n_] := 2*Sqrt[n!/(-1 + n)!]/n^2;
coeffk11[n_] := 2/n^(7/2);
Kx1[n1_, n2_] := 
  Kx1[n1, n2] = 
   If[n1 > n2, Kx1[n2, n1], 
    coeffk22[n2]*coeffk11[n1]*(1)*
     Total[Map[ffk[1/n1 + 1/n2, #[[1, 1]] + 1]*#[[2]] &, 
       CoefficientRules[
        Hypergeometric1F1[1 - n1, 2, 2*r/n1]*
         Hypergeometric1F1[1 - n2, 2, 2*r/n2], r]]]];

fk[a_, n_] = 
  Integrate[r^n*Exp[-a*r], {r, 0, Infinity}, 
   Assumptions -> {a > 0, n > 0}];
coeffk2[n_] := 2*Sqrt[n!/(-1 + n)!]/n^2;
coeffk1[n_] := (-1)*4/n^(3/2);
Kx2[n1_, n2_] := 
  Kx2[n1, n2] = 
   If[n1 > n2, Kx2[n2, n1], 
    coeffk2[n2]*coeffk1[n1]*(1)*
     Total[Map[fk[1/n1 + 1/n2, #[[1, 1]] + 1]*#[[2]] &, 
       CoefficientRules[
        Hypergeometric1F1[1 - n1, 2, 2*r/n1]*
         Hypergeometric1F1[1 - n2, 2, 2*r/n2], r]]]];
Kx[n1_, n2_] := -1/2*(N[Kx1[n1, n2]] + N[Kx2[n1, n2]])
KK = Table[Kx[n1, n2], {n1, 1, nmax}, {n2, 1, nmax}];

(*potential energy*)

f[a_, n_] = 
  Integrate[r^n*Exp[-a*r], {r, 0, Infinity}, 
   Assumptions -> {a > 0, n > 0}];
coeff[n_] := 2*Sqrt[n!/(-1 + n)!]/n^2;
(*-1/r*)
Px1[n1_, n2_] := 
  Px1[n1, n2] = 
   If[n1 > n2, Px1[n2, n1], 
    coeff[n1]*coeff[n2]*(-1)*
     Total[Map[f[1/n1 + 1/n2, #[[1, 1]] + 1]*#[[2]] &, 
       CoefficientRules[
        Hypergeometric1F1[1 - n1, 2, 2*r/n1]*
         Hypergeometric1F1[1 - n2, 2, 2*r/n2], r]]]];
(*-(1/(2*r))*Exp[-1.5r]*)
Px2[n1_, n2_] := 
  Px2[n1, n2] = 
   If[n1 > n2, Px2[n2, n1], 
    coeff[n1]*coeff[n2]*(-(1/2))*
     Total[Map[f[1.5 + 1/n1 + 1/n2, #[[1, 1]] + 1]*#[[2]] &, 
       CoefficientRules[
        Hypergeometric1F1[1 - n1, 2, 2*r/n1]*
         Hypergeometric1F1[1 - n2, 2, 2*r/n2], r]]]];

Px[n1_, n2_] := N[Px1[n1, n2]] + N[Px2[n1, n2]];

PP = Table[Px[n1, n2], {n1, 1, nmax}, {n2, 1, nmax}];

EE = Min[Eigenvalues[KK + PP]]

Out[268]= -0.671174 

(III)

In[1]:= ClearAll["Global`*"]

nmax = 10;

(*3d oscillator basis (l=0)*)
Psi[r_, n_] := (-1)^n *Exp[-1/2 *r^2] *Sqrt[2 *n!/Gamma[n + 3/2]]* 
   LaguerreL[n, 1/2, r^2];
(*kinetic energy*)
Kk[r_, n1_, n2_] := 
  FullSimplify[
   Psi[r, n2]*
    Laplacian[Psi[r, n1], {r, \[Theta], \[Phi]}, "Spherical"]];
Kx[n1_, n2_] := -1/2*Integrate[Kk[r, n1, n2]*r^2, {r, 0, \[Infinity]}];
KK = Table[Kx[n1, n2], {n1, 0, nmax}, {n2, 0, nmax}];

(*potential energy*)
VH1[r_] := -1/r;
VH2[r_] := -(1/(2*r))*Exp[-r*1.5];
Px1[n1_, n2_] := 
  Integrate[Psi[r, n2]*VH1[r]*Psi[r, n1]*r^2, {r, 0, \[Infinity]}];
Px2[n1_, n2_] := 
  Integrate[
   Psi[r, n2]*FullSimplify[VH2[r]]*Psi[r, n1]*r^2, {r, 
    0, \[Infinity]}];
Px[n1_, n2_] := Px1[n1, n2] + Px2[n1, n2];

PP = Table[Px[n1, n2], {n1, 0, nmax}, {n2, 0, nmax}];
EE = Min[Eigenvalues[KK + PP]]

Out[13]= -0.697315

(IV)

In[358]:= ClearAll["Global`*"]

nmax = 10;

(*3d oscillator basis*)
(*Psi[r_,n_]:=(-1)^n *Exp[-1/2 *r^2] *Sqrt[2 *n!/Gamma[n+3/2]]* \
LaguerreL[n,1/2,r^2];*)

(*kinetic energy*)
nmx = nmax;
hm = ConstantArray[0, {nmx + 1, nmx + 1}];
Do[hm[[1 + n1, 1 + n2]] = 
  KroneckerDelta[n1 - n2]*(n1 + 3/4) - 
   Sqrt[(2 (n1)^2 + (n1))/8] KroneckerDelta[n1 - n2 - 1];
 hm[[1 + n2, 1 + n1]] = hm[[1 + n1, 1 + n2]];, {n1, 0, nmx}, {n2, 0, 
  n1}]
KK = N[hm];

(*potential energy*)

f0[a_, n_] = 
  Integrate[r^n*Exp[-r^2]*Exp[-a*r], {r, 0, Infinity}, 
   Assumptions -> {a > 0, n > 0}];
f[a_, n_] := f[a, n] = N[f0[a, n], 100];
coeff[n_] := (-1)^n*Sqrt[2*n!/Gamma[n + 3/2]];
(*Coulomb part VH1[r]*)
Px1[n1_, n2_] := 
  Px1[n1, n2] = 
   If[n1 > n2, Px1[n2, n1], 
    coeff[n1]*coeff[n2]*(-1)*
      Total[Map[Function[{x}, f[0, x[[1, 1]] + 1]*x[[2]]], 
        CoefficientRules[
         LaguerreL[n1, 1/2, r^2]*LaguerreL[n2, 1/2, r^2], r]]] // N];

(*VH2[r]*)
Px2[n1_, n2_] := 
  Px2[n1, n2] = 
   If[n1 > n2, Px2[n2, n1], 
    coeff[n1]*coeff[n2]*(-(1/2))*
      Total[Map[Function[{x}, f[-1.5, x[[1, 1]] + 1]*x[[2]]], 
        CoefficientRules[
         LaguerreL[n1, 1/2, r^2]*LaguerreL[n2, 1/2, r^2], r]]] // N];

Px[n1_, n2_] := Px1[n1, n2] + Px2[n1, n2];
PP = Table[Re[Px[n1, n2]], {n1, 0, nmax}, {n2, 0, nmax}];
EE = Min[Eigenvalues[KK + PP]]

Out[371]= -0.697315

Answer 1

The matrix method is an approximation. The Hamiltonian you want to diagonalize acts on a hilbert space of infinite dimensionality. By choosing a finite basis, you make an approximation. Thus it is only natural that the eigenvalues that you recover depend on your choice of approximation.

Your question indicates that you might be aware of that, but still think that you took a sufficient amount of basis functions. You write:

We can see that an increase in the number of the basis functions doesn't change the energy much, therefore the ground state energy is ≈−0.67

The logic that you would like to apply here is that, having taken enough basis functions, you should eventually converge to the true ground state energy. That is right: The limit of your matrix-method result (in using the complete Basis) should yield the true ground state energy. However:

• The values that you postet don't show a sign of convergence at all. Look at your numbers: You can't conclude from any of them that you have reached a value that won't change anymore.
• You can't (similar argument) argue that the changes for bigger n are small - They are small compared to what? Remember that your state space is infinite-dimensional.
• Even if your calculation would have reached a convergence point for a finite n, let's say at n = 100, you can't conclude from this alone that you have reached the groundstate. It might still be that at n = 1000, you will get different values again. This will heavily depend on the set of groundstates that you use.

Long story short: You know that your calculation converges to the true ground-state energy, but only in the limit of fully exhausting your hilbert-space. However you can't conclude from your calculation that you have already reached this convergence.