# Computing solutions to the radial Schroedinger equation?

• A

## Main Question or Discussion Point

Hi all, I'm trying to compute the solutions to a general case for a Schroedinger equation with a radial potential but I'm stuck on a rather small detail that I'm not sure about. It's well known that I can perform the change of variables to spherical coordinates and express the radial part of the wavefunction as:

$\left(-\frac{\hbar^2}{2m}\frac{\partial^2}{\partial r^2}+V(r)+\frac{\hbar^2l(l+1)}{2mr^2}\right)\chi(r) = E\chi(r)$

with

$\chi(r) = rR(r)$

and $R(r)$ radial part of the solution. Now this works wonderfully for me because the only thing I'm interested in is the radial density, so basically $\chi^2$, which I think means I don't even have to worry too much with issues of low precision for small $r$. I already have a Numerov solver for the 1D equation so I thought I'd apply that here. I compute the solutions for a potential $V_{rot}(r) = V(r) + \frac{\hbar^2l(l+1)}{2mr^2}$ imposing as conditions that $\chi(r) = 0$ both at zero and infinity (actually some high but finite r value). What I'm a bit perplexed about though is the degeneracy and the role of $l$. My understanding of the problem is that I have to compute these solutions for various values of $l$, and each of them is going to provide a number of states, all with degeneracy $2l+1$, which then I can put together and sort by energy. I also assume I can apply a cutoff on $l$ (for example only compute solutions for $l\leq5$) since high $l$ also brings high energy and I only care for the low energy solutions. Is this the right way to proceed? It sounds a bit iffy because it uses different quantum numbers from the typical hydrogen atom solution, so I can't tell how that would work out exactly. Thanks!

Related Atomic and Condensed Matter News on Phys.org
BvU
Homework Helper
It's well known that I can perform the change of variables to spherical coordinates
Sure, but shouldn't you come up with a different equation ? Looks to me the $\partial^2\over dr^2$ isn't all there should be ...
 oversight on my part ? $\chi = rR$

BvU
Homework Helper
Sorry I got you off the 'unanswered' list -- but aside from that: for the hydrogen atom I remember the normalizability requires a series solution to break off and that imposes a limit on $|l| \le n$. How does that work in your general case ?

Sorry I got you off the 'unanswered' list -- but aside from that: for the hydrogen atom I remember the normalizability requires a series solution to break off and that imposes a limit on $|l| \le n$. How does that work in your general case ?
Good question, I don't know. I was wondering exactly about whether something like this was happening and if so how can it be dealt with numerically without any assumptions on $V(r)$.

BvU
Homework Helper
Seems like a tall order to me. Maybe @blue_leaf77 has an idea ?
@A. Neumaier mentions a famous package, but I suppose you need to give it $V(r)$
Ah, no, he suggests Gaussian but for what I know that's simply a molecular DFT package. I know many of those but they aren't generic solvers, they solve chemical systems, using specific potentials representing atoms. Plus they work with a 3D grid, which... I could do I guess, it's just that solving radially in 1D is that much more efficient if I can get the method right!

blue_leaf77
Homework Helper
Good question, I don't know. I was wondering exactly about whether something like this was happening and if so how can it be dealt with numerically without any assumptions on $V(r)$.
There is no guarantee that the radial Schroedinger equation is solvable analytically for any form of $V(r)$. As you have mentioned, you can try numerical solution but that requires you to input the explicit form of $V(r)$. I would suggest that you use finite difference method which transform your problem into a matrix eigenvalue problem and hence can be solved using the available eigenvalue subroutines. Before you do this of course you need to set $l$ which I think should only be integer numbers (i.e. no half-integer ones) since this will ensure that $2\pi$ rotation will return your original wavefunction (if half-integer $l$ is used, $2\pi$ rotation will give a negative sign to the rotated wavefucntion). You can start from $l=0$, and calculate the energies and eigenfunctions. For each $l$ you input, you should get the number of $n$ that can accomodate this value of $l$. As I recall, the cutoff for $l$ for each $n$ is the same as that for hydrogen atom but they are in general no longer degenerate.

BvU
You can start from $l=0$, and calculate the energies and eigenfunctions. For each $l$ you input, you should get the number of $n$ that can accomodate this value of $l$. As I recall, the cutoff for $l$ for each $n$ is the same as that for hydrogen atom but they are in general no longer degenerate.
Thanks! This is the key point for me. You mean you think I should only take the first eigenstate for $l=0$, the first two for $l=1$ and so on? Do you know how can I test for the validity of this or if there is somewhere I can consult as a reference? About the solution existing, I believe that should be the case as the potential that I'm working with is confining (either a quartic or higher order polynomial well).

blue_leaf77
Homework Helper
I never tried it myself but quoting the book "Physics of Atoms and Molecules" by Bransden and Joachain
For example, we shall see in chapter 7 and 8 that many properties of alkali atoms can be understood in terms of the motion of a 'single' valence electron in a potential which is central, but which deviate from $1/r$ Coulomb behaviour because of the presence of the 'inner' electrons. As a result, the energy of this valence electron does depend on $l$ and the degeneracy with respect to $l$ is removed, leading to $n$ distinct levels $E_{nl}$ for a given principal quantum number $n$.
The red part indicates that for a given $n$ there indeed should be $n$ numbers of $l$ which most likely also start from $l=0$.

I don't think the boundary condition χ(0) = 0 is valid. Consider hydrogenic s-orbitals.

I also assume I can apply a cutoff on l (for example only compute solutions for l≤5) since high l also brings high energy
For the hydrogen atom, at least, the energy levels are independent of l and are determined solely by n (in the basic, non-relativistic analysis).

blue_leaf77
Homework Helper
I don't think the boundary condition χ(0) = 0 is valid. Consider hydrogenic s-orbitals.
$\chi(r)$ is the actual radial wavefunction times $r$.

I don't think the boundary condition χ(0) = 0 is valid. Consider hydrogenic s-orbitals.
$\chi(r)$ is the actual radial wavefunction times $r$.
Yeah, this. If you had to recover the full wavefunction I'd have to compute $R(r) = \chi(r)/r$, but the problem with that is that the precision goes to hell for small values of $r$. However I only need to compute the expectation value of a radial quantity, which means I don't need to do that at all and $\chi$ is enough for me (since $\chi(r)^2=R(r)^2r^2$ which is the radial density).

The red part indicates that for a given $n$ there indeed should be $n$ numbers of $l$ which most likely also start from $l=0$.
But what are $n$ and $l$ in this context? Does that mean that the solutions for $l=0$ would be the equivalent of my $s$ orbitals for $n=1,2,3...$, the solutions for $l=1$ the $p$ orbitals for $n=2,3,4...$ etc.?

blue_leaf77
Homework Helper
But what are $n$ and $l$ in this context? Does that mean that the solutions for $l=0$ would be the equivalent of my $s$ orbitals for $n=1,2,3...$, the solutions for $l=1$ the $p$ orbitals for $n=2,3,4...$ etc.?
$n$ and $l$ are the same quantum numbers as those in hydrogen atoms, they also correspond to the same operators only that now the degeneracy in $l$ is removed. The degeneracy in $m$ however is still intact since the system has no preferential orientation in space.

$n$ and $l$ are the same quantum numbers as those in hydrogen atoms, they also correspond to the same operators only that now the degeneracy in $l$ is removed. The degeneracy in $m$ however is still intact since the system has no preferential orientation in space.
Ok, but in the hydrogenic atom solution, $n$ is a parameter. However in this case there's no obvious $n$. There's $l$, appearing in the centrifugal term of the modified potential $V_{rot}$. Now when I solve for that I get a number of eigenstates (ideally infinite ones, in practice as many as the grid size I integrated the thing on). But the index of the eigenstate can't be $n$, or I get to the absurd notion of having states with $n=0$ for any $l$. There's no obvious reason to think that for example the first solution for $l=1$ should be discarded; in theory it's as good a solution as any other. It would be different of course if I were perturbing the hydrogen atom solution, as I could use $n$ and $l$ to index the solutions and then perturb those, but if I'm starting from scratch there's no obvious way to do it. Unless it has to do with the limit solutions for small $r$ (which ought to share the same properties since the centrifugal term dominates, at least in my use case in which there are no $r^{-3}$ or even bigger terms).

A. Neumaier
I also assume I can apply a cutoff on l (for example only compute solutions for l≤5) since high l also brings high energy and I only care for the low energy solutions.
This is a Sturm-Liouville problem; thus the energies of the $k$th eigenvalue are strictly increasing in $\ell$. This gives an upper bound on the number of eigenvalues you need to try, once you have enough eigenvalues for $\ell=0$ and the lowest ones for $\ell=1,2,3,\ldots$.
The most reliable numerical procedure is to use a geometrically spaces grid (fixed ratios), discretize using cubic splines multiplied by the correct power of $r$ to match the conditions at the origin, and solve the associated matrix eigenvalue banded problem using LAPACK.

Ok, so I think I figured this out. I tried my code on the hydrogen atom problem and managed to map the functions to the known orbitals, so that clarifies things a bit. Pic incoming:

So what's happening here is we have eigenvectors for the 5 lowest energy solutions I found plus eigenvalues for the first 20 on the right. I've checked against analytical formulations and the orbitals match decently (with the obvious issues of precision). Here's the general approach I used:

1. I computed the solution to the above mentioned equation for values of $l$ ranging from 0 to 5. Each channel gives me a number of solutions (in my case, 800, since that's the size of the grid I used); let's index them with $i$ starting from 1;
2. I sorted the solutions by energy;
3. In the plot on the left, the value of $l$ is of course the one employed for the given channel; the value of $n$ though is in fact $n=l+i$;
4. As it can be seen in the plot on the right the states are roughly degenerate; we first have one dot for 1s, then two for 2s, 2p, and so on; numerical precision issues ruin the pattern pretty quickly though.
Basically what I got is that when I solve for $l=0$ I do in fact get all s orbitals, when I solve for $l=1$ I get all p orbitals, and so on. This means that the first solution for $l=1$, namely $l=1, i=1$ is in fact the 2p orbital; and so on.

So I think applying the same approach everywhere else should be fine, though precision is tricky (I'm using an existing algorithm as I mentioned but perhaps I can find one with better basis sets as @A. Neumaier suggested). What happens in non-Coulomb cases, as long as there are bound states, is simply that the so called 'accidental degeneracy' is removed.

DrDu
If you only write down a radial equation, l does not have to be integer. For example, it could stem from the coulombic system restricted to a conical region of space.
You will get infinitely many solutions for a chosen value of l. For small values of r, you can show that the solutions have to behave as $r^s$. How is s related to l?
This $r^s$ behaviour makes it difficult to use the numerov algorithm directly, as for higher s, both the wavefunction and its derivative vanish at r=0.

If you only write down a radial equation, l does not have to be integer. For example, it could stem from the coulombic system restricted to a conical region of space.
You will get infinitely many solutions for a chosen value of l. For small values of r, you can show that the solutions have to behave as $r^s$. How is s related to l?
This $r^s$ behaviour makes it difficult to use the numerov algorithm directly, as for higher s, both the wavefunction and its derivative vanish at r=0.
To clarify, I'm writing only the radial equation because I'm interested only in the radial density part of the solution (since I need it to compute the expectation value of a quantity that we assume has spherical symmetry too), but the potential is spherical. Therefore the usual logic with $l$ integer and the angular part being defined by $Y_{l,m}$ applies.

And yes, I know I get infinitely many solutions for a given value of $l$ - see above, I mentioned that I did (though they were limited by my grid). The question of the derivatives might be troublesome. I'm using this variation of the Numerov algorithm:

https://www.physics.wisc.edu/~tgwalker/448-9Mathematica/448 Mathematica/MatrixNumerov.pdf

meaning I'm basically diagonalising a matrix, the basis set is made of $\delta(x)$ functions and the derivatives are computed numerically to $O(d^6)$. I don't know if this still suffers from the same problems as the traditional, iterative Numerov method.

A. Neumaier
I don't know if this still suffers from the same problems as the traditional, iterative Numerov method.
You can check the error by varying the grid size and comparing the results.

You can check the error by varying the grid size and comparing the results.
Tried that, they're pretty consistent. Though one thing that messes with the results is if I start the grid too close to zero (zero itself is, of course, excluded, as it would make the potential diverge). I suppose numerical precision screws it up.

A. Neumaier
Tried that, they're pretty consistent. Though one thing that messes with the results is if I start the grid too close to zero (zero itself is, of course, excluded, as it would make the potential diverge). I suppose numerical precision screws it up.
That's because of the singularity at zero. It is probably best to pick grid points in a geometric progression.

That's because of the singularity at zero. It is probably best to pick grid points in a geometric progression.
I will try that at some point I guess but I'll need to redefine the derivatives then. As they are now they won't work without a regularly-spaced grid.

A. Neumaier
As already mentioned, the right way is to use cubic splines multiplied by the power of $r$ needed to cancel the singularity. These are $C^2$, so they can be differentiated analytically.