Numerical ground state wavefunction of Schrödinger equation with a Coulomb potential in 2D from NDEigensystem

Question

I want to numerically solve the ground state wave function of the hydrogen atom with the Coulomb potential using the NDEigensystem. Here is the code to get the ground state wave function from the NDEigensystem with "MaxCellMeasure" -> 0.01

 V[r_] := - 2/Sqrt[ r.r]
[ScriptCapitalL] = - !(
*SubsuperscriptBox[([Del]), ({x, y}), (2)](f[x, y])) + 
   V[{x, y}] f[x, y];
 d = 10; n = 1;
A = ImplicitRegion[-d <= x <= d && -d <= y <= d , {x, 
        y}];
{vals, funs} =  NDEigensystem[{[ScriptCapitalL], 
        DirichletCondition[f[x, y] == 0, True]}, 
      f, {x, y} [Element] A, n, 
   Method -> {"SpatialDiscretization" -> {"FiniteElement", 

{"MeshOptions" -> {"MaxCellMeasure" -> 0.01}}}}] ;

If I consider values very close to each other for "MaxCellMeasure", the wave function undergoes significant changes, and none of the resulting wave functions are valid. A diagram illustrating different "MaxCellMeasure" values is presented below.

![ (https://stackoverflow.com/image.jpg)

How can I find the correct wave function?

Answer 1

The NDEigensystem[...,66] finds about 8 negative eigenvalues all greater than -0.5. The true lower bound is -4. The plots of all solutions found are not ground states without nodes, they have higher angular momentum components.

   sol = NDEigensystem[{-Laplacian[f[x, y], {x, y}] - 2/Sqrt[x^2 + y^2] f[x, y], 
     DirichletCondition[  f[x, y] == 0, (  x == 10 || x == -10 || y == 10 || y == -10)]}, 
      f, {x, y} \[Element] Rectangle[{-10, -10}, {10, 10}], 66];
negs=Select[Sort[Transpose[sol],(#1[[1]]<#2[[1]]&)],(Negative[#[[1]]] &)];
First /@ negs
{-0.442905, -0.442905, -0.396269, -0.135906, -0.120594, -0.0935093, -0.0935093,
    -0.0495637}

The state of lowest energy has angular maomentum 1, that has a node and a zero at $(x,y)=0$ , supresssing the potentially negative infinite potential.

      Plot3D[negs[[1, 2]][x, y], {x, -10, 10}, {y, -10, 10}, 
                     PlotRange -> All, Axes -> None, Boxed -> False]

Define the energy functional

  Q[f_] := NIntegrate[-f Laplacian[f, {x, y}]-2/Sqrt[x^2 + y^2] f^2, 
          {x, -10, 10}, {y, -10, 10}] /
           NIntegrate[f^2, {x, -10, 10}, {y, -10, 10}]
Q[E^(-2 Sqrt[x^2 + y^2])]
  -4


Q[negs[[1, 2]][x, y]] // Quiet
   -0.454647

The boundary condition is very good approximated by the free ground state.

   Plot[E^(-2 Sqrt[x^2 + 10^2]), {x, -10, 10}]

The algebraic ground state $e^{-2 r }$ without a boundary

  ( -Laplacian[#, {r, \[Phi]}, "Polar"] - 2/r # == -4 # &)[  E^(-2 r)] 
    True 

So for this problem the numeric eigensystem algorithm is useless. If one uses the mesh conditions above that miss the boundary as an integer multiple, the results are complete chaos if plotted 3d.

The coulomb problem in 2D polar coordinates with Hamilton operator

$$-\partial_{r,r} \ \psi -\partial_{r} \ \psi - r^{-2}\partial_{\phi,\phi} \ \psi -2/r \ \psi $$

is solved in the same way as the 3D problem, with a different dimension term $r^{1-d} \partial_r$ only. Solution of the eigenvalue problem by ansatz $$\psi=e^{i m \phi} \ r^n \ \chi(r) $$ yields

$$-\partial_{r,r} \ \chi + m^2 r^{-2} \ \chi -2/r \ \chi +1/p^2 \chi=0 $$

The algebraic eigenvectors are confluent hypergeometric functions, the regular solutions at $r=0$ are Laguerre functions that are square integrable, if they have the form $$\chi = e^{-r/p} \ L^\alpha_n$$

The ground state has $m=0, n=0,\ L_0(r)=\text{const}$ and for all numerical approximations this solution is zero on the square $x,y=\pm10$

Nevertheless, there is a problem, the Hamitonian on the half line for $m=0$ is not a selfadjoint operator. Problem: The image vector $$H \quad\left( r^m \ e^{i m \phi} \ e^{-r/p}\right) $$ has to be square integrable on $\mathbb R$ with measure of integration $d\mu=r \ dr \ d\phi$. This excludes m=0. The particle may hit the center and disapear. In 3D with measure $d\mu=r ^2 dr$. the case $m=0$ is included.

The question is, if, for any $m$ in the sub-Hilbert space, the norm integral at $r=0$ converges and allows for partial integration there to obtain a selfadjoint operator, bounded from below, yielding a minimal energy eigenvector by the Riesz variation principle.

$$-\infty < \int_0^\infty \left(\ ( \ \partial_r\ (r^m e^{-r}))^2 + r^{2m} e^{-2r} \ m^2 \ r^{-2} - r^{2m} e^{-2r} r^{-2} \right) \ r^{\text{dim}-1} \ dr < \infty$$ such that the potential term is balanced by the kinetic term at least in mean.

Bottom line: The Coulomb problem in 2D is not the quantum image of the Kepler problem reduced by fixing the L-axis. This naive appoache simply degenerates to nonsense for L=0.

The set of Kepler ellipses of fixed energy and fixed greater axis degenerates into a spherical line bundle of line segments emerging from the center in in R^3, meaning the particle is reflected at the origin, but he direction of the reflected ray is not fixed by any physical law

(Pauli amd Dirac mended this classical hole by the spin. There is no j=0 massive stable particle. A spinning classical particle is reflected back to it's incoming ray by conservation of total angular momentum).

The spherical harmonics times radial waves degenerate to a spherical system of radial polynomial waves with a polynomial factor generating the nodes.

This concept of building tensor product spaces of the Hilbert space functions in 1d by separation of variables for $n>=3$. For $n=2$ a hard core repulsing force ar $r=0$ is necessary to prevent falling into the pit for $m=0$