Question

Solve Schrödinger equation for hydrogen atom?

Answer 1

Talk about a tall order. Obviously (?) this will be an obscenely long answer. More like a tome.

Note that each wave function component given as part of the overall wave function is unnormalized, except for when explicit expressions are presented.

SUMMARY

1. Separation of variables into `r` and `(theta,phi)`
2. Separation of variables of `(theta,phi)` into `theta` and `phi`
3. Solving the `phi` part
4. Solving the `theta` part
5. Putting together the `phi` and `theta` parts
6. "Solving" the `r` part
7. Putting together the overall wave function

I would hope you know that you only need to know how to use the wave functions (how to normalize them, how to show they are orthogonal, how to find the radial and angular nodes, etc), and maybe kind of know basic steps of the following derivation.

But see below anyway, because it took me 4+ hours...

---

The following derivation was adapted from here and from Physical Chemistry: A Molecular Approach by McQuarrie & Simon.

0. INITIAL DEFINITIONS

We begin from the time-independent Schrodinger equation (SE)

`hatHpsi = Epsi`,

which for hydrogen atom, has the Hamiltonian `hatH` defined in spherical coordinates to be:

`hatH = -ℏ^2/(2mu) nabla^2 - e^2/(4piepsilon_0r)`,

where:

• `ℏ = h//2pi` is the reduced Planck's constant.
• `mu = (m_p m_e)/(m_p + m_e)` is the reduced mass of the electron and proton.
• `nabla^2 = 1/r^2 (del)/(delr) (r^2 (del)/(delr)) + 1/(r^2 sintheta) (del)/(del theta)(sintheta (del)/(del theta)) + 1/(r^2sin^2theta) (del^2)/(delphi^2)` is the Laplacian operator in spherical coordinates.
• `e` is the elementary charge, `1.602 xx 10^(-19)` `"C"`, e.g. positive for a proton, negative for an electron.
• `epsilon_0 = 8.854187817 xx 10^(-12) "F"cdot"m"^(-1)` is the vacuum permittivity.
• `r` is the radial position. Under the Born-Oppenheimer approximation, we assume the nucleus is fixed so that `r` becomes the radial distance of the electron from the nucleus.

1. SEPARATION OF VARIABLES

We first have to get the SE into standard 2nd order partial differential equation form. Plug in the Laplacian:

`-ℏ^2/(2mu)[1/r^2 (del)/(delr) (r^2 (del)/(delr)) + 1/(r^2 sintheta) (del)/(del theta)(sintheta (del)/(del theta)) + 1/(r^2sin^2theta) (del^2)/(delphi^2)]psi - e^2/(4piepsilon_0r) psi = Epsi`

Subtract over `Epsi`, multiply through by `-(2mu)/(ℏ^2)`, and distribute `psi`...

`1/r^2 (del)/(delr) (r^2 (delpsi)/(delr)) + 1/(r^2 sintheta) (del)/(del theta)(sintheta (delpsi)/(del theta)) + 1/(r^2sin^2theta) (del^2psi)/(delphi^2) + (2mu)/(ℏ^2)(E + e^2/(4piepsilon_0r))psi = 0`

We assume separability of the wave function such that `psi_(nlm_l)(r,theta,phi) = R_(nl)(r)Y_(l)^(m_l)(theta,phi)`, where `R` is the radial component and `Y` is a spherical harmonic.

Since `R = R(r)` and `Y = Y(theta,phi)`, the partial derivatives become regular derivatives and functions of `r` are floated out as constants (similar with functions of angles):

`Y/r^2 (del)/(delr) (r^2 (delR)/(delr)) + R/(r^2 sintheta) (del)/(del theta)(sintheta (delY)/(del theta)) + R/(r^2sin^2theta) (del^2Y)/(delphi^2) + (2mu)/(ℏ^2)(E + e^2/(4piepsilon_0r))RY = 0`

Multiply through by `r^2/(RY)` from the left to complete the first partial separation of variables:

`overbrace(1/R (d)/(dr) (r^2 (dR)/(dr)))^("Radial Component") + overbrace(1/Y 1/(sintheta) (del)/(del theta)(sintheta (delY)/(del theta)) + 1/Y 1/(sin^2theta) (del^2Y)/(delphi^2))^"Angular Component" + overbrace((2mu r^2)/(ℏ^2) (E + e^2/(4piepsilon_0r)))^"Radial Component" = 0`

By subtraction, these two functions are equal to each other:

`1/R (d)/(dr) (r^2 (dR)/(dr)) + (2mu r^2)/(ℏ^2) (E + e^2/(4piepsilon_0r)) = -[1/Y 1/(sintheta) (del)/(del theta)(sintheta (delY)/(del theta)) + 1/Y 1/(sin^2theta) (del^2Y)/(delphi^2)]`

But since they are functions of different variables, that can only work if they are equal to the same separation constant, which we arbitrarily designate `-lambda` (assign the dangling negative sign to the radial part):

`-[1/R (d)/(dr) (r^2 (dR)/(dr)) + (2mu r^2)/(ℏ^2) (E + e^2/(4piepsilon_0r))] = -lambda`

`1/Y 1/(sintheta) (del)/(del theta)(sintheta (delY)/(del theta)) + 1/Y 1/(sin^2theta) (del^2Y)/(delphi^2) = -lambda`

Now we rearrange these ordinary differential equations separately. Take the radial component first. Switch the signs and then multiply through by `R` from the left. Subtract over `lambdaR` to then get:

`ul((d)/(dr) (r^2 (dR)/(dr)) + (2mu r^2)/(ℏ^2) (E + e^2/(4piepsilon_0r))R - lambdaR = 0)`

An analogous process follows for the angular equation and we obtain:

`ul(1/(sintheta) (del)/(del theta)(sintheta (delY)/(del theta)) + 1/(sin^2theta) (del^2Y)/(delphi^2) + lambdaY = 0)`

2. FURTHER SEPARATING THE ANGULAR PART

The angular part is still a function of `theta` and `phi`, so we need to separate this further. Assume that `Y_(l)^(m_l) (theta,phi) = Theta(theta)Phi(phi)`. Repeat the same process as seen in the first separation of variables.

`Phi/(sintheta) (del)/(del theta)(sintheta (delTheta)/(del theta)) + Theta/(sin^2theta) (del^2Phi)/(delphi^2) + lambdaThetaPhi = 0`

Multiply through by `(sin^2theta)/(ThetaPhi)` to complete the separation of variables:

`overbrace((sintheta)/(Theta) (d)/(d theta)(sintheta (dTheta)/(d theta)) + lambdasin^2theta)^(Theta) + overbrace(1/Phi (d^2Phi)/(dphi^2))^(Phi) = 0`

As before, these are equal to a separation constant. Call it `-B`. With some minor manipulation, we obtain as before (assign the dangling negative sign to the `theta` part):

`ul((sintheta)/(Theta) (d)/(d theta)(sintheta (dTheta)/(d theta)) + lambdasin^2theta - B = 0)`

`ul(1/Phi (d^2Phi)/(dphi^2) + B = 0)`

3. SOLVING THE PHI PART

Obviously, the `Phi` part is the simplest to solve to obtain `B`.

`(d^2Phi)/(dphi^2) + BPhi = 0`

We assume a solution to this ordinary differential equation with constant coefficients by choosing `phi` on the `xy` plane:

`Phi = e^(im_lphi)`

This gives the auxiliary equation:

`-m_l^2 cancel(e^(im_lphi))^(ne 0) + Bcancel(e^(im_lphi))^(ne 0) = 0`

`=> sqrtB = pm |m_l|`

The solution is then a linear combination of `Phi`'s,

`Phi(phi) = c_1e^(im_lphi) + c_2e^(-im_lphi)`

For one full rotation on the `xy` plane in one direction (say, clockwise), the wave function must coincide with itself (be cyclic), so `Phi(0) = Phi(2pi)`.

To eliminate variables, currently, `m_l >= 0` by assumption. So, let `c_1 + c_2 = A` and allow `m_l` to be negative, so that:

`Phi(phi) = (c_1 + c_2)e^(im_lphi) = Ae^(im_lphi)`

From Euler's relation,

`e^(i cdot m_l cdot 0) = cos(0m_l) + isin(0m_l) = e^(i cdot m_l cdot 2pi) = cos(2pim_l) + isin(2pim_l)`

This only holds true for integer values of `m_l`, which we define to be the magnetic quantum number:

`m_l = 0, pm1, pm2, . . .`

The unnormalized `phi` component to the angular wave function is then:

`color(green)(barul(|" "stackrel(" ")(Phi(phi) prop e^(im_lphi))" "|))`

Note that we currently do not know that `|m_l| <= l`, but it can be proven using ladder operators as shown here.

(That fact would be determined the traditional way at the end of solving the `Theta` part of the SE.)

4. SOLVING THE THETA PART

As the separation constant is consistent throughout the process, `B = m_l^2` here too. Multiply through from the left by `(Theta)/(sin^2theta)`:

`1/(sintheta) (d)/(d theta)(sintheta (dTheta)/(d theta)) + (lambda - m_l^2/(sin^2theta))Theta = 0`

The "trick" for this is to use a function of `costheta`. Let:

`P(costheta) := Theta(theta)`
`x := costheta`

This transformation changes the differential operator by the chain rule from `d/(d theta)` to `(dx)/(d theta) d/(dx) = -sintheta d/(dx)` so that we get:

`d/(dx)(sin^2theta (dP)/(dx)) + (lambda - m_l^2/(sin^2theta))P = 0`

Since `sin^2theta = 1 - cos^2theta = 1 - x^2`:

`d/(dx)((1 - x^2) (dP)/(dx)) + (lambda - m_l^2/(1 - x^2))P = 0`

And after the product rule and chain rule to rewrite in terms of just `x`, we get the ordinary differential equation as a function of `costheta`:

`ul((1 - x^2) (d^2P)/(dx^2) - 2x(dP)/(dx) + (lambda - m_l^2/(1 - x^2))P = 0)`

This is something called an associated Legendre-type differential equation. This is a known equation type with solutions called the Associated Legendre Polynomials `P_(l)^(|m_l|)(costheta)`.

I'm not going to bother with the general formula for the solution, but the first few are:

`ul(l = 0: )` `" "" "" "" "" "" "ul(l = 1: )`
`P_(0)^(0) = 1` `" "" "" "" "" "P_(1)^(0) = costheta`
`" "" "" "" "" "" "" "" "P_(1)^(1) = sintheta`

`ul(l = 2: )`
`P_(2)^(0) = 1/2 (3cos^2theta - 1)`
`P_(2)^(1) = 3sinthetacostheta`
`P_(2)^(2) = 3 - 3cos^2theta`

What we further get from this is that `ul(lambda = l(l+1))` is our separation constant from the very first separation of variables, and that `m_l = {-l, -l+1, . . . , l-1, l}`, but we knew the latter from general chemistry.

5. PUTTING TOGETHER THE ANGULAR PARTS

Now that we have two solutions `P_(l)^(|m_l|)(costheta)` and `Phi(phi)`, we construct the form of the unnormalized spherical harmonics `Y_(l)^(m_l)(theta,phi)`:

`color(green)(barul(|Y_(l)^(m_l) prop P_(l)^(|m_l|)(costheta) e^(im_lphi)|))`

You can look at the first few normalized spherical harmonics here:

`Y_(0)^(0) = 1/(4pi)^(1//2)" "" "" "" "" "Y_(1)^(0) = (3/(4pi))^(1//2) costheta`
`Y_(1)^(1) = (3/(8pi))^(1//2) sintheta e^(iphi)" "" "Y_(1)^(-1) = (3/(8pi))^(1//2) sinthetae^(-iphi)`

and you can find the rest in your textbook.

6. "SOLVING" THE RADIAL PART

Now that we finally have `lambda`:

`(d)/(dr) (r^2 (dR)/(dr)) + (2mu r^2)/(ℏ^2) (E + e^2/(4piepsilon_0r))R - l(l+1)R = 0`

Apply the product rule on the first half, factor out `R`, then divide by `r^2` from the left to get a standard form:

`ul((d^2R)/(dr^2) + 2/r (dR)/(dr) + [(2mu)/(ℏ^2)(E + e^2/(4piepsilon_0r)) - (l(l+1))/r^2]R = 0)`

The solution we care about here is near the nucleus. Most textbooks don't show any more detail than this...

Here, `R(r)` takes the following unnormalized form, which utilizes the associated Lageurre Polynomials `L_(n+1)^(2l+1)((2r)/(na_0))` as a function of the principal quantum number `n` and radial position `r`:

`color(green)(barul(|R_(nl)(r) prop r^l e^(-r//na_0)L_(n+1)^(2l+1)((2r)/(na_0))|))`

where `a_0 = "52.9177 pm"` is the bohr radius.

The first few associated Laguerre polynomials are shown here:

`n = 1, " "l = 0, " "L_(1)^(1)(x) = -1`
`n = 2, " "l = 0, " "L_(1)^(2)(x) = -2!(2-x)`
`" "" "" "" "l = 1, " "L_(3)^(3)(x) = -3!`

and the rest should be in your textbook.

7. PUTTING TOGETHER THE OVERALL WAVE FUNCTION

Thus (thus?), from the first substitution,

`color(blue)(barul(|psi_(nlm_l)(r,theta,phi) = R_(nl)(r)Y_(l)^(m_l)(theta,phi)|))`

we know the form of the hydrogen atom wave functions. (I hope you recognize that none of the above green rectangled equations are normalized.)

I will list the normalized wave functions `n = 1` and `n = 2`, but your textbook surely holds many more... Let `sigma = r//a_0`.

`(n,l,m_l) = (1,0,0): " "" "psi_(100) = 1/sqrtpi (1/a_0)^(3//2) e^(-sigma)`

`(n,l,m_l) = (2,0,0): " "" "psi_(200) = 1/sqrt(32pi)(1/a_0)^(3//2) (2-sigma)e^(-sigma//2)`

`(n,l,m_l) = (2,1,0): " "" "psi_(210) = 1/sqrt(32pi)(1/a_0)^(3//2) sigma e^(-sigma//2) costheta`

`(n,l,m_l) = (2,1,pm1): " "psi_(21pm1) = 1/sqrt(64pi)(1/a_0)^(3//2) sigmae^(-sigma//2) sintheta e^(pmiphi)`

And lastly, if you wish to find the energy based on a particular atomic orbital wave function, evaluate the triple integral in allspace:

`E = int_(0)^(2pi) int_(0)^(pi) int_(0)^(oo) psi_(nlm_l)^"*"(r,theta,phi) hatH psi_(nlm_l)(r,theta,phi)r^2dr sinthetad thetadphi`

In atomic units, the ground-state energy for the hydrogen atom is `-"0.5000 Hartrees"`, or `-"13.61 eV"`, the potential energy of its `1s` orbital.