Question

How do you use spherical Bessel functions of the first and second kind to solve the following ordinary differential equation?

`(d^2R(r))/(dr^2) + 2/r (dR(r))/(dr) + [(2mE)/(ℏ^2) - (l(l+1))/(r^2)]R(r) = 0`

where `R(r)` is to be a radial wave function that vanishes at the boundary `r = R`.

This is the differential equation I solved for when I was working with an infinite spherically symmetric well with a bound electron. By multiplying through by `r^2`, this is a spherical Bessel differential equation:

`r^2(d^2R(r))/(dr^2) + 2r (dR(r))/(dr) + [(2mE)/(ℏ^2)r^2 - l(l+1)]R(r) = 0`

with `k^2 = (2mE)/(ℏ^2)`.

Another reference is:
http://farside.ph.utexas.edu/teaching/qmech/Quantum/node81.html

Answer 1

Using the substitution that `k^2 = (2mE)/(ℏ^2)`, the common way to solve this seems to be to let `rho = kr`, so that:

• `(rho/k)^2(d^2R(r))/(d(rho//k)^2) = rho^2(d^2R(r))/(drho^2)`

• `2(rho/k) (dR(r))/(d(rho//k)) = 2rho (dR(r))/(drho)`

This would then lead to

`r^2 (d^2R(r))/(dr^2) + 2r(dR(r))/(dr) + [(2mE)/(ℏ^2)r^2 - l(l+1)]R(r) = 0`

becoming:

`rho^2(d^2R(rho))/(drho^2) + 2rho (dR(r))/(drho) + [rho^2 - l(l+1)]R(r) = 0`

or

`(d^2R(rho))/(drho^2) + 2/rho (dR(r))/(drho) + [1 - (l(l+1))/rho^2]R(r) = 0`

From the second reference, the solutions to this are some linear combination of spherical Bessel functions of the first and second kind, `j_l(rho)` and `y_l(rho)`, where `l` is the angular momentum of the state.

These are given by:

`j_l(rho) = rho^l (-1/rho d/(drho))^l ((sin rho)/rho)`

`y_l(rho) = -rho^l (-1/rho d/(drho))^l ((cos rho)/rho)`

However, `y_l(rho)` is not square-integrable at `rho = 0`:

graph{(cosx)/x [-10, 10, -5, 5]}

So, we can choose only `j_l(rho)`. In this case, I was asked to only look for `s` states, so `l = 0` and the unnormalized wave function is:

`j_0(rho) = rho^0 (-1/rho d/(drho))^0 ((sin rho)/rho)`

`= (sin rho)/rho`

Applying the boundary condition that `R(r) = 0` at `r = R` (the sphere radius),

`j_0(kR) = (sin (kR))/(kR) = 0`

and this is only true when `kR -= rho_(nl) = npi`, where `n = 1, 2, 3, . . .`.

So, `rho_(nl) = npi = sqrt((2mE)/ℏ^2)R`. This allows us to solve for the energies of the `1s, 2s, 3s, . . .` states in terms of the roots of `j_0(rho)`:

`color(blue)(E_(n0) = (ℏ^2rho_(nl)^2)/(2mR^2))`

This also leads to the unnormalized radial wave function:

`R(r) = (sin (npir//R))/(npir//R)`

Normalizing this in spherical coordinates, one finds:

`1 = N^2 int_(0)^(2pi) 1^2 d phi int_(0)^(pi) sintheta d theta int_(0)^(R) R^"*"(r) R(r) r^2 dr`

`1 = N^2 int_(0)^(2pi) d phi int_(0)^(pi) sintheta d theta int_(0)^(R) (sin^2 (npir//R))/(npir//R)^2r^2dr`

`= 4pi N^2 int_(0)^(R) (sin^2 (npir//R))/(n^2pi^2//R^2cancel(r^2))cancel(r^2)dr`

`= (4piR^2)/(n^2pi^2) cdot N^2 int_(0)^(R) sin^2 (npir//R)dr`

Using the identity `sin^2u = 1/2(1 - cos2u)`, one would get:

`= (4piR^2)/(2n^2pi^2) cdot N^2 (R - cancel(R/(2npi)sin((2npiR)/R))^(0))`

Therefore:

`N = sqrt((2n^2pi^2)/(4piR^3))`

`= 1/sqrt(4pi)sqrt(2/R)(npi)/R`

`= 1/sqrt(2piR)(npi)/R`

As a result, the normalized wave function would be:

`color(blue)(R(r)) = 1/sqrt(2piR) cancel((npi)/R) sin(npir//R)/(cancel(npi)r//cancel(R))`

`= color(blue)(1/sqrt(2piR) sin(npir//R)/r)`

And this looks something like this:

graph{sin(x)/x [-0.5, 20, -1, 1.2]}