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, 
{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?

  • $\begingroup$ Welcome to Mathematica StackExchange! How are you plotting your wavefunctions? Note that they are not radially symmetric, because your domain is a square! Look at these 3D plots. You see that I get the same function (minus the negative sign, but that doesn't matter). Compare also with circular domain. $\endgroup$
    – Domen
    Commented Nov 6, 2023 at 14:15
  • 3
    $\begingroup$ Aren't the eigenfunctions only unique up to an arbitrary factor ? In contrast to the eigenvalues which are unique! $\endgroup$ Commented Nov 6, 2023 at 14:41
  • $\begingroup$ @Domen-Thanks. This event also occurs for ` "MaxCellMeasure" -> 0.007` in circular domain see i.sstatic.net/IV7k5.png. $\endgroup$
    – AminD
    Commented Nov 6, 2023 at 16:57
  • $\begingroup$ @AminD, which version of Mathematica are you using? Are you sure you are not changing n in between the runs? Can you start with a fresh kernel? $\endgroup$
    – Domen
    Commented Nov 6, 2023 at 17:02
  • $\begingroup$ @UlrichNeumann- Yes, but these eigenfunctions up to an arbitrary factor are not correct. $\endgroup$
    – AminD
    Commented Nov 6, 2023 at 17:07

1 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,

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]

Numerical ground state

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])]


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


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

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

Boundary error

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

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

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$


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.