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

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

• Why are you computing matrix elements for each term in $H$ separately? It seems like you could speed things up by computing each matrix element of $H$ with a single call to NIntegrate. Commented May 7, 2023 at 4:49
• Did you read my answer on physics.stackexchange.com/questions/753794/… ? Commented May 7, 2023 at 6:55
• @Ghoster, thanks! I do this to make it easier to find mistakes Commented May 7, 2023 at 8:38
• @Alex Trounev, thanks for the comment! I have a question, unlike Slater Type Orbitals, Gaussian basis functions depend on coefficients $\alpha_i$ and $\beta_i$ (in the case of cylindrical system) which are generally unknown, so it is not entirely clear how to apply the minimization method you have shown to this basis. If it's not difficult for you, then please demonstrate. Commented May 7, 2023 at 8:50
• @MamMam Why not to discuss this question under my answer? :) Commented May 7, 2023 at 16:03

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$$.

• Yes, that's right! It's better to discuss here ). My question was related to excited states, i.e. and other eigenvalues of the Hamiltonian. To get them, I need to build a Hamiltonian matrix already on a known function (which was obtained in your code) and diagonalize it? Or can the remaining eigenvalues be obtained more quickly already from your code? I'm sorry for asking, but I don't see it clearly yet. Commented May 7, 2023 at 16:30
• @MamMam This question is not so clear. Could you add a00, b00 to compute excited states? Then we can test how it works. Commented May 8, 2023 at 2:35
• I try to reformulate question. Your code finds such parameters c[i] that the energy of the system with the given parameters a[i] and b[i] is minimal. After the minimization process, we have the minimum energy value and its wave function (c[i]). We can also find the minimum energy and corresponding wave function by diaganizing the matrix. In this case, we will have several eigenvalues and eigenvectors (depending on the number of basis functions), the second largest eigenvalue is the energy of the 1st excited state, etc. Commented May 8, 2023 at 6:45
• @MamMam Actually, my answer is completed. I can recommend you this report arxiv.org/pdf/2009.13556.pdf . After reading, please open new topic with question about excited states. Commented May 8, 2023 at 7:33
• Very nice and complete answer! @MamMam Please avoid the reformulation of questions. It makes it very difficult for an uninvolved but interested reader to follow. Another small comment on the scope of your questions. You try to learn physics, but this is a Mathematica forum. I encourage you to discuss this type of questions with your scientific advisor and focus here solely on the technicalities. Commented May 8, 2023 at 9:21

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.

• Thanks! But your answer seems to be highly theorized, perhaps it answers my question, but as a beginner researcher it is difficult for me to understand this answer. I would like a more direct answer to my question. In my case, I just add one additional function to the basis and the answer is wrong, while the approach remains the same. Commented May 6, 2023 at 19:53
• Background of methods added Commented May 7, 2023 at 6:33
• I think that the potential is $1/\sqrt{r^2 + z^2}$ and not $1/r$ as you assumed. Commented May 8, 2023 at 9:32