# Numerically solving radial Schrödinger equation with Yukawa potential

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

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

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.

• I think the shift is everything, because but it seems like NDEigensystem looks for eigenvalues closest to zero first, then the next closest to zero, then the next, etc. Since, for the hydrogen atom, the eigenvalues bunch up near zero, it's never going to look for the ground state first! This is the purpose of the shift. It's there to shift the eigenvalues up enough so that the smallest eigenvalue (the ground state energy) is closest to zero, and so it will look for those first. Dec 5, 2023 at 21:28
• With your code, if you look for 32 eigenvalues instead of 3, you'll find that the eigenvalues go from -0.960029 all the way up to 0.25 (which last corresponds to the first excited state). However, if you look for 33 eigenvalues, it will then find the value 1, which corresponds to the ground state. Dec 5, 2023 at 21:28
• With that in mind, I suggest that this is a duplicate of 105298, since the fix of shifting the eigenvalues is shown there. Dec 5, 2023 at 21:29
• Further to @march's point: Note that in the top answer to the second question you link to, the code actually finds the 20 eigenvalues closest to zero and then selects the seven lowest-energy ones. The full result of vals in that code includes several positive eigenvalues, which are discarded by the Take command. Also note that the top answer doesn't actually find the ground state; it only gets the $n = 2$ through $n = 8$ states. Dec 5, 2023 at 22:03

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.} *)


• This answers my main question, but I should note that there was a sign error -- the eigenvalues you should get should be {1, 0.25, 0.1111}, etc. However, I suppose it's better to include the minus signs so -1 is the lowest eigenvalue. Thanks! Dec 6, 2023 at 8:27
• Ah, I see. If you want to use the original sign, then you'll need to switch the sign of lambda0 as well; it should be greater than the upper bound on the spectrum of $-H$. Dec 6, 2023 at 12:29
• This is really a great explanation of the 'Shift' functionality. If you do not mind, I'd like to adopt this for the documentation. Is that OK? Dec 6, 2023 at 12:30
• @user21: Be my guest! Dec 6, 2023 at 12:32

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.