Numerically solving radial Schrödinger equation with Yukawa potential

Question

I am trying to solve the radial Schrödinger equation using NDEigensystem but I am running into some issues. There are posts about doing this (see here for example), but I don't particularly understand the specifics such as why a shift is necessary for unbounded potentials like the Coulomb potential (here the Coulomb solutions were found without a shift).

Regardless, I've non-dimensionalized the Schrödinger equation, with a Yukawa potential, to be of the form

\begin{equation} \frac{d^2f}{dx^2} + \left[ \frac{2}{x}e^{-x/\xi} -\gamma^2 - \frac{\ell(\ell + 1)}{x^2} \right]f = 0 . \end{equation}

Here, $f$ is the radial part of the wave function multiplied by the radius, and $x$ is a nondimensionalized radius.

The energy eigenvalue is given by $\gamma^2$ and in the limit that $\xi\to \infty$, $\gamma \to 1/n$ which are the non-dimensional energy eigenvalues for a plain Coulomb potential.

As a warm up, I've simply tried to solve the Coulomb case (ignoring the possible variable substitutions suggested here). The code for this is as follows:

H[l_] = f''[x] + (2/x - (l (l +  1))/x^2)*f[x];

{vals, funs} = 
NDEigensystem[{H[0], DirichletCondition[f[x] == 0, True]}, 
f[x], {x, 0, 100}, 3, 
Method -> {"SpatialDiscretization" -> {"FiniteElement", \
{"MeshOptions" -> {"MaxCellMeasure" -> 0.001}}}}];
Take[Sort[vals]]

which outputs

{-0.0185335, -0.00331214, 0.00932567} 

This is obviously significantly different from the {1, 0.25, 0.111111} that I expect. I've assumed $\ell = 0$ for simplicity, but naturally, I would also like to recover the $\ell\geq 0$ results as well. Thanks in advance.

Answer 1

First, I think you have a sign error. The eigenvalues $\gamma^2$ being returned by the code will be proportional to $-E$ for the system, not $E$. This turns out to make a difference in what follows.

The "shift" trick is necessary because by default, NDEigensystem returns

the $n$ smallest magnitude eigenvalues and eigenfunctions for the linear differential operator $\mathcal{L}$ over the region $\Omega$.

(bolding mine) In other words, NDEigensystem is giving you the eigenvalues closest to zero that it can find.

So if you want to find a bunch of negative eigenvalues with large magnitudes, the simplest trick is to shift the operator. Specifically, if we know that the eigenvalues are bounded below by $\lambda_0$, then we can solve the problem for the operator $\mathcal{L} - \lambda_0$, and the smallest-magnitude eigenstates for $\mathcal{L} - \lambda_0$ will be the lowest-energy states for $\mathcal{L}$. Note that this doesn't have anything in particular to do with the potential being unbounded below; you'd have to do this for any system for which the ground state energy is negative (e.g. a particle in a finite square well with $V = 0$ outside the square well.)

So here's the modified code that gives the expected results. Note the use of the "Shift" suboption within the "Eigensystem" option, which performs this shift internally and then undoes it when it returns the actual eigenvalues.

H[l_] = -f''[x] - (2/x - (l (l + 1))/x^2)*f[x] ;

lambda0 = -2;
{vals, funs} = 
  NDEigensystem[{H[0], DirichletCondition[f[x] == 0, True]}, 
   f[x], {x, 0, 100}, 3, 
   Method -> {"Eigensystem" -> {"Arnoldi", "Shift" -> lambda0}, 
     "SpatialDiscretization" -> {"FiniteElement", {"MeshOptions" -> \
{"MaxCellMeasure" -> 0.001}}}}];
vals
Plot[funs, {x, 0, 40}, PlotRange -> All]

(* {-0.111111, -0.25, -1.} *)

Answer 2

The problem seems to be, that the exponential decay cannotb be fixed at $x=100$. The minimal eigenvalue decreases with reduction of the interval.

Using the physical scaling of the Schrödinge operator

   H[l_] = -(1/2) f''[x] + (l (l + 1))/(2 x^2) f[x] - 1/x f[x];

  {vals, funs} = 
     NDEigensystem[{H[0], DirichletCondition[f[x] == 0, True]}, 
                   f[x], {x, 0, 12}, 4, 
                   Method -> {"SpatialDiscretization" -> {"FiniteElement", 
                            {"MeshOptions" -> {"MaxCellMeasure" -> 0.0001}}}}];
     Sort[{vals, funs}\[Transpose], (#1[[1]] < #2[[1]] &)] 

    Plot[Evaluate[Rest /@ sorted], {x, 0, 12}, 
           PlotLegends -> Placed[First /@ sorted, Scaled[{0.6, 0.2}]]]

As you see, you will get reasonable results only, if the condition at $x=\infty$ is adapted to a small multiple of the last zero. Here the second eigenfunction does not show the exponential decay, that shifts the wight of the function to the left, lowering the eigenvalue. If one enlarges the interval, the lowest eigenvalue moves up to zero. Conclusion: Method is not adapted to infinite intervals.