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

Question

It's a part 2 of this question Optimal basis set of Gaussian functions for describing a quantum system (part 1)

There, the answer was given to a question related to finding geometric progressions of the coefficients ($\alpha_i$ and $\beta_i$) of a Gaussian basis set in the general case using minimization. But there was still the second question regarding why the authors of the article use different indexes for the parameters $N$ and $M$ for $\alpha$ and $\beta$ respectively.

Brief description of the problem: The system has the following 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).

To find the ground and excited states of the system by the matrix method, one can use the Gaussian set of basis functions $\psi_j=e^{-b_{j} z^2}e^{-a_{j} \rho^2}$, where $a_{j}$ and $b_{j}$ are parameters. This approach is described in the article C. Aldrich & R.L. Greene (https://drive.google.com/file/d/1q74hAn0UAdNd8DtPkoCsdYPkr61-xhr2/view).

The authors use as sets of parameters geometric progressions in the following form:

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

Why do authors use different indices for parameters (N and M for $\alpha$ and $\beta$ respectively)? Aren't they paired in the basis function?

Answer 1

In the paper "Hydrogen-Like Systems in Arbitrary Magnetic Fields. A Variational Approach" the method to optimize basis functions has been proposed. They described basis functions dependent on 4 quantum numbers in a form

$\psi_{ijq}=z^q e^{-a_i r^2} e^{-b_j z^2}, \psi_m=r^{|m|}e^{i m\theta}$

where $a_i,b_j$ are parameters to be optimized. The variational wave functions in terms of a set of basis functions can be written as follows $\Psi_{mq}=\sum_{ij}\psi_{ijq}\psi_m$

The code to optimize energy is not so differ from what we discussed here:

nmax = 3; jmax = 4;
B = 5; Lz = m = 0; q = 0; ms = -1/2;
VP1[r_, z_] := -1/Sqrt[r^2 + z^2] + 1/8 B^2 r^2 + B/2 (m + 2 ms);

Psi[r_, z_, i_, j_, q_] := z^q Exp[-b[j]*z^2]*Exp[-a[i]*r^2]; 
Psim[r_, \[Theta]_, m_] := r^Abs[m] Exp[I m \[Theta]];
(*kinetic energy*)
Kk = 1/2 FullSimplify[
    Psi[r, z, i2, j2, q] Psim[r, \[Theta], -m]*
      Laplacian[
       Psi[r, z, i1, j1, q] Psim[r, \[Theta], m], {r, \[Theta], z}, 
       "Cylindrical"] + 
     Psi[r, z, i1, j1, q] Psim[r, \[Theta], m]*
      Laplacian[
       Psi[r, z, i2, j2, q] Psim[r, \[Theta], -m], {r, \[Theta], z}, 
       "Cylindrical"]];

ss = Integrate[Kk r, {r, 0, \[Infinity]}, {z, -Infinity, Infinity}, 
   Assumptions -> {a[i1] > 0, b[i1] > 0, a[i2] > 0, b[i2] > 0, 
     a[j1] > 0, b[j1] > 0, a[j2] > 0, b[j2] > 0}];
Kx = - 1/2 2 Pi Sum[
    c[i1, j1] c[i2, j2] ss, {i1, nmax}, {i2, nmax}, {j1, jmax}, {j2, 
     jmax}];


(*potential energy*)
Px = Integrate[
   Psi[r, z, i2, j2, q] Psim[r, \[Theta], -m]*VP1[r, z]*
    Psi[r, z, i1, j1, q] Psim[r, \[Theta], m]*r, {r, 
    0, \[Infinity]}, {z, -Infinity, Infinity}, 
   Assumptions -> {a[i1] > 0, b[i1] > 0, a[i2] > 0, b[i2] > 0, 
     a[j1] > 0, b[j1] > 0, a[j2] > 0, b[j2] > 0}];

Px = 2 Pi Sum[
    c[i1, j1] c[i2, j2] Px, {i1, nmax}, {i2, nmax}, {j1, jmax}, {j2, 
     jmax}];



int = Integrate[ 
   Psi[r, z, i2, j2, q] Psim[r, \[Theta], -m]*
    Psi[r, z, i1, j1, q] Psim[r, \[Theta], m] r, {r, 
    0, \[Infinity]}, {z, -Infinity, Infinity}, 
   Assumptions -> {a[i1] > 0, b[i1] > 0, a[i2] > 0, b[i2] > 0, 
     a[j1] > 0, b[j1] > 0, a[j2] > 0, b[j2] > 0}];
norm = {2 Pi Sum[
      c[i1, j1] c[i2, j2] int, {i1, nmax}, {i2, nmax}, {j1, 
       jmax}, {j2, jmax}] == 1};

 
U = 2 Pi ArrayFlatten[
    Table[ int, {i1, nmax}, {i2, nmax}, {j1, jmax}, {j2, jmax}], 2];

 
var = Join[{a1, b1, qa, qb}, 
  Flatten[Table[c[i, j], {i, nmax}, {j, jmax}]]]; a[1] = a1; 
b[1] = b1; Do[a[i] = a[i - 1] qa, {i, 2, nmax}]; Do[
 b[j] = b[j - 1] qb, {j, 2, jmax}]; con = 
 Join[norm, {a1 > 0, b1 > 0, qa > 0, qb > 0}];

 sol = NMinimize[{Kx + Px, con}, var]



(*Out[]= {-1.37834 + 0. I, {a1 -> 1.39282, b1 -> 19.9993, 
  qa -> 3.4335, qb -> 0.252195, c[1, 1] -> -0.00902664, 
  c[1, 2] -> 0.0164657, c[1, 3] -> -0.317789, c[1, 4] -> -0.375104, 
  c[2, 1] -> 0.0377631, c[2, 2] -> -0.209612, c[2, 3] -> -0.11075, 
  c[2, 4] -> 0.0505215, c[3, 1] -> -0.158495, c[3, 2] -> 0.0175225, 
  c[3, 3] -> 0.0210703, c[3, 4] -> -0.0143571}}*)

Please, pay attention that we use 3 a and 4 b to prepare system of $3\times 4=12$ basis functions. Eigenvalues can be computed as follows

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

 Hij = 
 Last@CoefficientArrays[(Kx + Px) /. abq, 
   Flatten[Table[c[i, j], {i, nmax}, {j, jmax}]], "Symmetric" -> True]

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

 val // Sort // Chop

(*Out[]= {-1.37834, 0.542836, 5.79876, 6.72021, 7.97674, 13.9054, \
32.4115, 34.3706, 35.4432, 40.6693, 41.9055, 67.6762}*)

The ground state is -1.37834, and it is not bad compared to exact value -1.380398866427.