Computing solutions to the radial Schroedinger equation?

In summary, The conversation discusses the computation of solutions to a general case for a Schroedinger equation with a radial potential. The change of variables to spherical coordinates is mentioned, as well as the expression for the radial part of the wavefunction. The use of a Numerov solver is suggested, as well as the computation of solutions for various values of ##l## with a cutoff. The issue of degeneracy and the role of ##l## is also discussed, with a suggestion to start from ##l=0## and calculate the energies and eigenfunctions for each ##l##. The cutoff for ##l## is mentioned as being the same as for the hydrogen atom, but with the possibility of non-degenerate states.
  • #1
Gan_HOPE326
66
7
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!
 
Physics news on Phys.org
  • #2
Gan_HOPE326 said:
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 ...
[edit] oversight on my part ? ##\chi = rR##
 
  • #3
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 ?
 
  • #4
BvU said:
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)##.
 
  • #6
BvU said:
Seems like a tall order to me. Maybe @blue_leaf77 has an idea ?
In the thread https://www.physicsforums.com/threa...inding-the-energy-levels.906937/#post-5712497
@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!
 
  • #7
Gan_HOPE326 said:
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 wavefunction). 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 accommodate 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.
 
  • Like
Likes BvU
  • #8
blue_leaf77 said:
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 accommodate 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).
 
  • #9
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##.
 
  • #10
I don't think the boundary condition χ(0) = 0 is valid. Consider hydrogenic s-orbitals.
 
  • #11
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).
 
  • #12
John Park said:
I don't think the boundary condition χ(0) = 0 is valid. Consider hydrogenic s-orbitals.
##\chi(r)## is the actual radial wavefunction times ##r##.
 
  • #13
John Park said:
I don't think the boundary condition χ(0) = 0 is valid. Consider hydrogenic s-orbitals.
blue_leaf77 said:
##\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).

About this:

blue_leaf77 said:
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.?
 
  • #14
Gan_HOPE326 said:
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.
 
  • #15
blue_leaf77 said:
##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).
 
  • #16
Gan_HOPE326 said:
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.
 
  • #17
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:

download.png

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.
 
  • #18
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.
 
  • #19
DrDu said:
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.
 
  • #20
Gan_HOPE326 said:
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.
 
  • #21
A. Neumaier said:
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.
 
  • #22
Gan_HOPE326 said:
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.
 
  • #23
A. Neumaier said:
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.
 
  • #24
Gan_HOPE326 said:
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.
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.
 

1. What is the radial Schroedinger equation?

The radial Schroedinger equation is a mathematical equation used in quantum mechanics to describe the behavior of a particle in a spherically symmetric potential. It is used to determine the wave function and energy levels of a particle in a given potential.

2. What are "computing solutions" to the radial Schroedinger equation?

Computing solutions to the radial Schroedinger equation refers to the process of using numerical methods or computer algorithms to solve the equation and obtain the wave function and energy levels of a particle in a given potential. This allows for a more accurate and efficient way of obtaining solutions compared to analytical methods.

3. Why is the radial Schroedinger equation important?

The radial Schroedinger equation is important because it allows us to understand and predict the behavior of particles in a spherically symmetric potential, which is crucial in fields such as quantum chemistry and atomic physics. It also provides a fundamental basis for understanding the principles of quantum mechanics.

4. What are some common techniques used to compute solutions to the radial Schroedinger equation?

Some common techniques used to compute solutions to the radial Schroedinger equation include the Numerov method, finite difference methods, and the shooting method. These methods involve breaking down the equation into smaller, solvable parts and using numerical algorithms to find the solutions.

5. Are there any limitations to computing solutions to the radial Schroedinger equation?

Yes, there are limitations to computing solutions to the radial Schroedinger equation. These include the use of approximations and assumptions in the numerical methods, which may introduce errors in the solutions. Additionally, the computational cost may increase significantly for more complex potentials or higher energy levels.

Similar threads

Replies
4
Views
1K
  • Advanced Physics Homework Help
Replies
1
Views
2K
Replies
2
Views
233
Replies
1
Views
668
Replies
3
Views
430
  • Advanced Physics Homework Help
Replies
1
Views
2K
Replies
2
Views
1K
  • Advanced Physics Homework Help
Replies
11
Views
973
Replies
16
Views
1K
  • Quantum Physics
Replies
27
Views
2K
Back
Top