Optimal basis set of Gaussian functions for describing a quantum system (part 1)

Question

This question arose during the discussion of the previous questions: Why does the minimum eigenvalue change dramatically when one basis function is added to the basis set? and Finding excited states using the condition of wave functions orthogonality

Brief description of the problem:

Consider the following system with the Hamiltonian: $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. Physically, this Hamiltonian describes a hydrogen atom in a magnetic field equal to 5 (in dimensionless units).
In order to solve this problem by the matrix method, it is necessary to choose a set of basis functions. To describe ground and excited states with the z component of the angular momentum m = 0 (1s, 2s, 3d (m = 0), 3s ...), a Gaussian basis set is well suited. Such basis functions have the following form: $\psi_j=e^{-b_{j} z^2}e^{-a_{j} \rho^2}$, where $a_{j}$ and $b_{j}$ are parameters. The task is reduced to the correct finding of parameters.

The parameters describing the ground state (1s) with high accuracy were found in this work (see here https://arxiv.org/abs/1709.05553 or here https://pubs.aip.org/aip/jcp/article/147/24/244108/195534/Accurate-and-balanced-anisotropic-Gaussian-type). But the basis presented in this article does not describe excited states (2s, 3d (m = 0), 3s ...) well, since there are not enough basis functions to describe them. The basis set parameters presented in this article has a complex structure, so increasing the basis set functions in this case is a difficult task.

It seems that a more general way to find the basis set of parameters is presented in this article (see here https://drive.google.com/file/d/1q74hAn0UAdNd8DtPkoCsdYPkr61-xhr2/view?usp=sharing or here https://onlinelibrary.wiley.com/doi/abs/10.1002/pssb.2220930140). The sets of parameters that the authors use are geometric progressions:
$\alpha_j=\alpha_{j-1}(\frac{\alpha_N}{\alpha_1})^{1/(N-1)}$

$\beta_j=\alpha_{j-1}(\frac{\beta_M}{\beta_1})^{1/(M-1)}$,
where $\alpha_1$, $\alpha_N$, $\beta_1$ and $\beta_M$ are the first and the last parameters. About finding these parameters, the authors write the following: "the parameters $\alpha_1$, $\alpha_N$, $\beta_1$ and $\beta_M$ were varied in a four dimensional optimization search which minimized the sum of the eigen-values of interest" and "We have performed the calculations presented in this paper with a basis with ten Gaussians in the $\rho$ direction and twelve Gaiissians in the z direction, for a total of 120 different basis functions.". I really don't understand what that means?

Energy values from the articles:
ground state (1s) = - 1.3793 (exact value = -1.380398866427)
the first excited state (2s) = -0.19335 (exact value = -0.193746709717)
the second excited state (3d (m=0)) = -0.07365
the third excited state (3s) = -0.03835

All these states should have the same structure of the basis set function ($\psi_j=e^{-b_{j} z^2}e^{-a_{j} \rho^2}$) due to the m=0 (the z component of the angular momentum).

I have a few questions:
I really can't understand what the authors mean in varied optimization search of the parameters $\alpha_1$, $\alpha_N$, $\beta_1$ and $\beta_M$, could you please explain to me what that means?

Why do the parameters ($\alpha_1$, $\alpha_N$, $\beta_1$ and $\beta_M$) have different indexes, shouldn't their number be the same? Since they are paired in the basis function.

I am well aware that this question may not be very typical, but it should be considered in the context of the previous two questions, since this question is a continuation of the reasoning that is demonstrated there. On the other hand, a possible answer will be demonstrated in Mathematica, which has optimization and minimization tools.

I would be glad to any comments on solving the problem.

I would like to draw attention to the fact that to get the energy from the Aldrich, C., & Greene, R. L article it's needed the value from the table 1 to subtract the field value and divide by 2 (this is due to the fact that the authors of this article do not take into account the spin). For example for 2s state: value from the table 1 for magnetic field = 5 is equal 4.6133 so the energy of the 2s state = (4.6133 - 5)/2 = -0,19335

Answer 1

To compute exited states $\Psi_i$ with some basis functions $\Phi_j$ we need to solve equation $H\Psi=\lambda U\Psi$, where $H_{ij}=\langle\Psi_i|\hat{H}|\Psi_j\rangle$, and $U_{ij}=\langle\Phi_i|\Phi_j\rangle$. Last matrix is not inversible in general case, so we have some numerical problem here. For example, exited states of the ground state computed here can be calculated with Eigensystem[] as follows

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];
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];
mat = (Kx + Px) /. 
    Join[Table[a[n] -> a0j[[n]], {n, nmax}], 
     Table[b[n] -> b0j[[n]], {n, nmax}]] // Chop;

matrix = 
 Last@CoefficientArrays[mat, Table[c[n], {n, nmax}], 
   "Symmetric" -> True]

U = Table[
   2 Pi int /. 
    Join[Table[a[n] -> a0j[[n]], {n, nmax}], 
     Table[b[n] -> b0j[[n]], {n, nmax}]], {l1, nmax}, {l2, nmax}];

mat = Inverse[U] . matrix; {val, vec} = Eigensystem[mat]; 

Therefore, the ground state and exited states are given by

val // Sort

(*Out[]= {-1.3804, -0.193692, -0.0640854, 0.0576317, 0.278105, \
0.641224, 1.27678, 2.35197, 4.03324, 5.46049, 9.06797, 16.7551, \
33.9461, 71.7555, 154.926, 337.651, 738.55, 1617.97, 3553.12, \
7850.89, 17666., 42031.}*) 

First and second state are in a good agreement with ground state (1s) = -1.380398866427, and the first excited state (2s) = -0.193746709717. But all other eigenvalues are differ from exact values. Apparently we need to extend our basis with constrain on U so that Inverse[U] exists.

Update 1. To implement main idea from the paper "Hydrogen-Like Systems in Arbitrary Magnetic Fields. A Variational Approach" we use code

nmax = 5;

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}]; U = 
 Table[2 Pi int, {l1, nmax}, {l2, nmax}]; 


var = Join[{a1, b1, qa, qb}, Table[c[i], {i, nmax}]]; a[1] = a1; 
b[1] = b1; Do[{a[i] = a[i - 1] qa, b[i] = b[i - 1] qb}, {i, 2, nmax}];

 sol = 
 NMinimize[{Kx + Px, {norm == 1., a1 > 0, b1 > 0, qa > 0, qb > 0}}, 
  var]


(*Out[]= {-1.37514 + 0. I, {a1 -> 1.24159, b1 -> 0.455892, 
  qa -> 2.12986, qb -> 3.77828, c[1] -> 0.49116, c[2] -> 0.318804, 
  c[3] -> 0.127836, c[4] -> 0.0669821, c[5] -> 0.0506232}}*)

Note, that we compute the ground state as sol[[1]]=-1.37514, and it is not so bad compared to exact value -1.380398866427. Other eigenvalues computed as follows

abq = Take[sol[[2]], 4];

 matrix = 
 Last@CoefficientArrays[(Kx + Px) /. abq, Table[c[n], {n, nmax}], 
   "Symmetric" -> True]

 mat = Inverse[U /. abq] . matrix; {val, vec} = 
 Eigensystem[mat];

val // Sort

(*Out[]= {-1.37514, 2.20185, 11.2994, 43.117, 159.438}*)  

The problem of this method is that even for nmax=8 we have badly conditioning matrix U.

Answer 2

Maxwell: There are three forms of fields: |E|=|B| light, |E|>|B| electric fields with a bit of movement of charges. |E|<|B| magetic fields of strong neutral current or worse spin fields on surfaces, with a bit of extra charges. They are easly discriminated by their 4-vector potential and these sectors are Lorentz invariant.

The Hamiltonian of the H-atom in a magnetic field has the following forms with B M = B_z L_z the constant magnetic energy of the angular momentum in state M

Cylindrical coordinates

$$ H = -r^{-1} \ \partial_r \left(\ r \ \partial_ r \psi \right) - \partial_{zz}\psi + B M \ \psi + C^2 \psi \ r^2 - \frac{1}{\sqrt{r^2+z^2}} \psi $$

Spherical coordinates

$$ H = -\frac{1}{r^2} \partial_r \left( r^2 \partial_r \psi\right) - \frac{1}{r^2 \sin \theta} \partial_{ \theta} (\sin \theta \partial_\theta \psi) - \frac{1}{r} \psi \ + \ B M \psi + C^2 r^2 \psi$$

This is a problematic Hamiltonian for a simple reason: Its good for small magnetic fields everywhere smaller than the electric field. It's ok if the central charge is not a point but is closed at $r=0$ by the oscillator potential of a constant charge density in the nucleus.

In a strong magnetic field, if there is a region where B > E eg on the surface of a neutron star, the electric field can locally transformed away yielding a magnetic field of helix form. The effect is known from northern aura, the wave functions concentrates along any magnetic line between energetic boundaries eg by a surface.

In this case, with $L_z = \pm M \hbar$ and no limit by $L^2$ , there is no lower bound. By emission of a photon and aquiring one more wavelength on the circle the energy can go down to -oo.

Its essential therefore to exactly fixing the numerical constants first. The rest can be seen in

Liberman Hydrogen Atom