This semester I’m teaching the lab portion of a course called Applied Math. My colleague in the math department teaches the lecture and I have a 1.5 hour lab once a week. One of the cool things about it is that all the students are currently taking Modern Physics as well. I wanted to put down some thoughts about how the students and I are using the shooting method to solve for the quantum wavefunctions of hydrogen.
The shooting method for solving for eigenfunctions of a 1D quantum well is described in most Modern Physics textbooks. Essentially you use a numeric differential equation integrator to start at one side of the well, integrate through the well and into the other side’s wall, and then look to see if the wavefunction runs away to . If it does, you can’t accept that solution because that would lead to an infinite probability that’s unnormalizable. So what do you do? Well, you have to change something. So you change the energy in the differential equation:
Here’s an example from a finite square well (2-nm wide, 2 eV deep) with three different energies: one that doesn’t run to infinity and two neighboring energies that do:
The energies for those curves are 0.64, 0.64154, and 0.65 eV. It’s clear that going away from the eigenenergy either way causes the wavefunction to run away to infinity.
I like teaching this 1D well shooting method because we can quickly adjust to non-flat wells and gain some insight into the expected wavefunctions. Just one example is the ramp well where the wavefunction has a higher amplitude where the potential is higher. We discuss how the probability is higher there because the particle’s kinetic energy is lower and so it moves slower and hence spends more time in that region.
Mathematica is a great tool for this. Here’s the code for generating the curves above. The arrow points to the part that has to be changed to do a more complicated well:
Simply changing that zero to something like Sin[x] immediately takes you to a problem that would be difficult to solve analytically.
What has been fun for the last few years in this course is to go beyond 1D wells and use the same shooting method to show how we can get the energy states of hydrogen. Students seem to be jazzed about this because I tell them that we don’t have to use any of the (wrong) assumptions built into the Bohr model, and we don’t have to do any of the (hard) series solutions to the differential equations involved (you know, the Laguerre and Legendre Polynomials, blah).
Ok, so first the hard part. Going from 1D wells to hydrogen is difficult because the two derivatives in the equation above become a 3D laplacian with derivatives with respect to r, theta, and phi (the spherical coordinates because the potential has spherical symmetry). Ok, just because it’s fun to use here on wordpress, I might as well show you the ugliness:
That’s not even the whole Schroedinger equation, just the second derivative. The big problem is the multi-variable nature of it. Really the extra crud in front of the various terms is not really an issue when you’re doing numeric work because Mathematica handles it all the same: tell me something about wavefunction at one point, I’ll take a small step and reevaluate. Repeat. That’s it!
So what do we do? The first thing is to assume that the three dimensions aren’t really mixed once you find a solution. That means we assume . The beauty of that assumption is, after doing a little algebra, you can get all the phi crud on one side and the theta and r stuff on the other. Ah, separation of variables, you’re my friend! Doing that gives you a separation constant. Voila! The m-quantum number is born.
That’s the easy one and students recognize the solutions (sines and cosines) right away. They also recognize the boundary conditions: when you go all the way around a circle, you better end up where you started. Next you can get all the theta crud on one side and the r crud on the other (with the m floating along, of course). HERE’S WHERE THIS METHOD STARTS TO TAKE A TURN. Sorry about the all caps but some people just scroll to these sorts of statements.
Normally, if you’re trying to show the students the analytical solution, you say (emphasis added) “ok, we’ve got theta and r separated. Let’s set both sides equal to a constant. How about l(l+1)?” Aha, the l (pronounced ‘ell’) quantum number is born!
That’s not what I do. I simply call it a constant and see what Mathematica can make of it. What I mean is, can students monkey with the value of that constant in Mathematica and find any values that keep the wavefunction from running away to infinity? If you’re getting tired of this post, I’ll give it away. They find that the following values of the separation constant give reasonable (non-infinite) curves for : 0, 2, 6, 12, 20, 30, . . . Yep, it’s l(l+1).
Here’s an example. This is a polar plot of the solution to the theta version of the differential equation after all the separations.
In this example the non-infinite solution is 6.000000 while the other two are 5.9 and 6.1.
Here’s the Mathematica code for that. It’s a little more sophisticated but really right on par with what my students can do in lab:
Note that I call the separation constant lambda and that the m-quantum number is present.
R equation (or, how you get to -13.6 eV)
Ok, now that we know what m’s and l’s are needed to keep the phi and theta parts from running to infinity, we’re ready for the R equation, which is the only one with the energy term in it.
Here it’s very similar to the 1D well approach except that there’s no hard wall at some known r>0 point. Instead the potential smoothly goes from negative infinity at r=0 to 0 at r=infinity (cool symmetry, huh?). But for my students, it really feels the same as the 1D well models. Here’s an image of the lowest energy state:
The image shows the best solution at -13.6 and the solutions at -12.6 and -14.6.
What’s really fun is having the students mess around with the l-quantum number. It completely changes the differential equation and yet the energies where you get non-infinite solutions remain at the same values. You can really only show why that is with the analytical approach but it’s fun to have them see it organically.
Just for completeness, here’s the Mathematica code for the R equation shooting method:
That’s it in a nutshell. Ok, that’s it in an incredibly long post that no one will actually make it to the bottom of. I really have fun with this stuff and I appreciate the power of Mathematica’s NDSolve and Manipulate commands. So how do you teach the hydrogen atom?