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.

**Shooting Method**

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.

**Hydrogen**

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?

Just a quick comment about units if they’re throwing you off. All energies are in eV and all lengths are in nm. Here are the relevant constants:

hbar c = 197 eV-nm

mc^2=511000 eV

ke^2 = 1.44 eV-nm

This is really nice. I’d have to go do it myself to understand all the details, of course. Pedagogically, it’s nice because it provides a mechanism for doing inquiry around these mathematical objects, which typically we just lecture on.

I used to use mathematica all the time. As an undergrad, mathematica was infused into nearly all our courses for majors… beginning sophmore year with modern physics and classical mechanics. This paid huge dividends later in upper division courses, and also when I was grad student. It’s certainly been a generative tool. You are reminding me of its importance and utility.

The Manipulate command that was introduced in the few years is really a game changer. The students can rapidly explore the parameter space and can produce code quickly that they are proud to share with others.

I, too, first learned Mathematica in undergrad but then took 10 years off from it. Now it’s a necessary tool for my teaching. It just sucks that it’s so expensive. Of course, the student version is “only” ~$150 and everything I do uses the student version to be sure that my students can do it also.

I made it to the end! (Ahhh, quiet Sunday mornings!) There is one move that I’m unfamiliar with and that is the symbol just to the left of the 0 and 2 in the first snippet of code (where you set the potential well values). How do you do that in Mathematica? What’s it called? I don’t even know how to look it up!

I wish we had Mathematica back when I was in college and grad school. I would have loved using it. I got nothing but lecture until I escaped to the lab. I have one HS student who uses it regularly (she’s been using differential equations to model the changes in Earth’s ice sheet as a result of orbital dynamics since she was a sophomore). She’ll be in my second year physics class next year and I hope to get her using Mathematica on some projects with a like-minded cohort of students.

Ah quite Mother’s day mornings when she is out doing a walk for cancer treatment fundraising and the boys are all working on cleaning up the garden for her!

The symbol you’re referring to is the PieceWise function. Hit esc-pw-esc and you get the curly bracket. For mine, type zero, then control comma then -1<x<1 then control enter to get a new line. I love it, it's much easier to use than a bunch of nested if-then statements.

So why do you expect \psi(-3) = \psi'(-3) ? I don’t see why this should be the case for a particular eigenfunction. Since you expect symmetry, why not use \psi(0) = 1, \psi'(0) = 0 , (when expecting even symmetry) and \psi(0) = 0, \psi'(0)=1 (when expecting odd symmetry)?

How do you justify these values to the students? One is of course a normalization, but the other is unknown.

It’s good that you seem to be at least discussing the initial values — but I assume that you recognize that your curve _cannot_ be an eigenfunction (only one of the initial values can be chosen arbitrarily, the other is determined), even if you knew the eigenvalue exactly in advance.

Numerically, you seem to be observing some stability — I’ll have to think about this a little bit. Do your eigenvalues match an analytically calculated value?

Here’s a variation I would be more comfortable with: use the same potential, but adjust it to be \infty for |x|>3 (or some such value — think of this as a “cut-off” potential). You expect the eigenfunction for the original potential to decay exponentially anyway, so the solution should be essentially 0 for |x| > 3. Now, your solution must have \psi(\pm 3) =0 . Choose the derivative arbitrarily at -3, and then look for a solution that is 0 again at +3. Note that this is essentially what you are doing, but it seems to me to be easier to justify the initial values. And you don’t have to worry about the embarrassment that _all_ of your numerical solutions (to your original problem) “blow up” if you plot far enough to the right!

It also might be useful to do some initial analysis — given the potential, can you estimate the rate of the exponential decay? You can use this to choose appropriate values for the “cut-off” potential.

Good points, Art. (full disclosure: Art is the math professor who teaches the lecture portion of the course) When the students adjust the initial conditions, as long as they are deep enough into the wall, the general shape of the eigenfunctions and the values of the eigenenergies don’t change. Why? Because the wavefunction is experiencing quite dramatic exponential decays in the walls of the well and as long as you start deep enough, the value is effectively zero and so is the slope. Now we can’t start at both value and slope equal to zero because then the solution is zero everywhere. But any other values are effectively very small because the diff eq allows for scaling by arbitrary constants (that’s the normalization piece).

With the students I have them monkey around with the initial conditions, and it’s interesting when they’re not deep enough into the wall. I like Art’s perspective, though, that we could do some simple analysis to figure out how deep is deep enough. That’s equivalent to the notion of a skin depth.

Please ignore the following, I just wanted to check if latex works in wordpress comments:

Pingback: Shooting method for hydrogen | SuperFly Physics | Solve Math & Science Problems - Solveable.com

Sweet Andy. I’m going to try this out in QM in the fall when we do the H atom. I like to give them a better reason than “because it works” that l(l+1) is a reasonable constant to try. Trying to get them to see l(l+1) from the constants that work seems like a much better approach.

Hi Joss. I’ll try speaking for Art here and say that the best way to proceed would be to find the list numerically but then say something like “aha, there must something really cool going on analytically for the acceptable answers to be this particular list.” If you have time to show the analytical solution, fine, but if not you can at least plant a seed about the power of an analytical solution. I like to point out that Schroedinger didn’t have Mathematica but he did have a really good grasp of differential equations.

Pingback: Trying out a new type of simulation-based pre-class assignment « Science Learnification

Thanks for your post. Very helpful for teaching purposes indeed. I’ll give it a shot next time I teach Radiative Processes, in 2016.