How to expand the existing basis set so that it becomes more complete?

Question

I have a set of anisotropic gaussian basis set which describes the ground state of the system with great accuracy.

The Hamiltonian of the system has the following form: $$H=-\frac{1}{2}\Delta-\frac{1}{r}+\frac{1}{2}\rho^2-1$$ where $r=(\rho,z,\phi)$ is a coordinate in the cylindrical system

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}$=0.500074, 0.502676, 0.510451, 0.507502, 0.52768, 0.625164, 0.935248, 1.826161, 3.845598, 8.293859, 18.014295, 39.20736, 85.3836, 185.975252, 405.095196, 882.399269, 1922.095421, 4186.829292, 9120.018288

$b_{0j}$=0.026284, 0.057253, 0.124713, 0.184064, 0.271658, 0.429143, 0.799226, 1.740925, 3.792194, 8.2604, 17.993331, 39.194226, 85.375371, 185.970096, 405.091965, 882.397245, 1922.094153, 4186.828497, 9120.01779

The exact value of the ground state energy from the literature is -1.0222124 and this basis set with great accuracy give this value (-1.02221). It's know the excited levels energies from the literature too: -0.174, -0.069, -0.036865, -0.022815, -0.015495. The basis set that I used give the follow values for the excited levels energies: -0.173374, -0.0199577, 0.312572, 1.18162... If the first excited level differs in the third decimal place, then the energies of the subsequent excited levels are completely different. This is because the basis set I am using is not complete and it is necessary to increase the number of basis functions in order to get more accurate values for the excited levels.

Is it possible, having a certain set of basic functions (in my case, it is 19 functions), to supplement it with additional functions so that it becomes more complete? How to do it? Because I would like to get the energies of the excited levels too.

In the code I rename $\rho\equiv r$

ClearAll["Global`*"]

(*the basis set*)
b00 = {0.026284, 0.057253, 0.124713, 0.184064, 0.271658, 0.429143, 
   0.799226, 1.740925, 3.792194, 8.2604, 17.993331, 39.194226, 
   85.375371, 185.970096, 405.091965, 882.397245, 1922.094153, 
   4186.828497, 9120.01779};
a00 = {0.500074, 0.502676, 0.510451, 0.507502, 0.52768, 0.625164, 
   0.935248, 1.826161, 3.845598, 8.293859, 18.014295, 39.20736, 
   85.3836, 185.975252, 405.095196, 882.399269, 1922.095421, 
   4186.829292, 9120.018288};
b0 [j_] := b00[[j]];
a0 [j_] := a00[[j]];

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

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]*1/2*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]*(-1)*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}];

EEE = Kxx + Pxx1 + KAx + KSx;
Sort[Eigenvalues[{EEE, Bxx}]](*Energies*)


Out[88]= {-1.02221, -0.173374, -0.0199577, 0.312572, 1.18162,
2.72596, 5.63718, 12.0333, 17.9735, 28.2709, 61.4758, 135.447,
298.338, 656.001, 1440.89, 3168.65, 7007.42, 15772.6, 37523.6}

Answer 1

Its a pure mathematical question.

In cylindrical coordinates with an atractive plane 1/r Coulomb potential and and a harmonic oscillator in the x,y plane, in states with plane integer angular momentum M yielding a repulsive potential

 M^2/r^2, M=0,  

the natural Hilbert space basis is

{ HermiteH[n, r] Exp[-r^2] e^i k z} 

because the system is free in z.

Of course any base in z will do, eg

 HermiteH[m] Exp[-z^2 ]

in order to have a base of normalizable orthogonal functions centered at z=0.

By my knowlege, the plane - r-> 1/r is not an electric potential. The plane central potential of a wire is

q Log (r/r0) 

with zero divergence

D[ 1/r D[r Log'[r]] ==0

while

2pi r dr * 1/r D[r (1/r)'] =2 pi 1/r dr = 2 pi r dr (1/r^2)

has a cylindrical charge density falling off with 1/r^2 on circles.