The properties of a chemical element are
determined by the available energy levels of an electron's possible state in
orbit around

the nucleus. (The element's chemical properties depend on how
tightly it holds onto its outer electrons and how readily extra electrons can
hitch a ride;
its absorption and emission
spectra depend on the energy levels of all of the electrons.) To determine
these, we must solve Schrödinger's
equation. (In writing the following I leaned heavily on Schiff's chapter
14.)

The electrostatic field of the nucleus is spherically symmetric, and we are concerned with solutions of definite energy, so we can skip directly to the solution of the spherically symmetric time-independent case; which reduces to a simple form.

We have an energy E and some natural b, with a total of 2.b+1 distinct solutions for the directional dependence yielding a typical solution, in terms of time t, radius r and direction v (encodable, for example, via a unit vector or a pair of angular co-ordinates),

- ψ(t,r,v) = exp(−i.E.t/ℏ).R(r).Y
_{b a}(v)

to our field equation, with R satisfying

- 2.m.r.r.(V(r) −E)/ℏ
^{2}+b.(1+b) = d(R'(r).r.r)/dr/R(r)

Since R's equation doesn't depend on v, each solution, for given E and
b, will yield 1+2.b distinct states of the object described by our wave
function, one for each possible value for a. Note that the whole derivation up
to this point only relied on V depending only on r. We now substitute in the
Coulomb potential V(r) = −Q/r, with Q some integer multiple of
e.e/4/π/ε_{0} = α.ℏ.c. Suppose E to be negative
(for a bound state), introduce a dimensionless radial co-ordinate ρ = r.K
with K = √(−8.m.E)/ℏ and radial wave function factor F defined
by F(ρ).exp(−ρ/2) = R(r); then d/dr = K.d/dρ, R'(r) =
(F'(ρ)−F(ρ)/2).exp(−ρ/2).K so
that

- d(R'(r).r.r)/dr
- = d((F'(ρ)−F(ρ)/2).exp(−ρ/2).ρ.ρ)/dρ
- and our differential equation becomes:
- d((F'(ρ)−F(ρ)/2).exp(−ρ/2).ρ.ρ)/dρ/F(ρ)/exp(−ρ/2)
- = ρ.ρ.(V(r) −E)/(−4.E) +b.(1+b)
- = b.(1+b) +ρ.ρ.(1/4 −Q.K/ρ/(−4.E))

so introduce λ = Q.K/(−4.E) = (Q/ℏ).√(m/2/(−E)) to express this as

- (b.(1+b) +ρ.ρ/4 −ρ.λ).F(ρ)
- = d((F'(ρ)−F(ρ)/2).exp(−ρ/2).ρ.ρ)/dρ/exp(−ρ/2)
- = −(F'(ρ)−F(ρ)/2).ρ.ρ/2 +(F''(ρ)−F'(ρ)/2).ρ.ρ +2.(F'(ρ)−F(ρ)/2).ρ
- = F''(ρ).ρ.ρ +(2−ρ).F'(ρ).ρ +F(ρ).(ρ.ρ/4 −ρ)

whence

- F''(ρ).ρ.ρ +(2−ρ).F'(ρ).ρ = F(ρ).(b.(1+b) −ρ.λ +ρ)

We now (just as when solving for the spherical harmonics) look for a polynomial solution, F = sum(: f(j).power(j) ←j :) for some ({scalars}:j:{integers}). This yields

- 0
- = F''(ρ).ρ.ρ +F'(ρ).(2−ρ).ρ −F(ρ).(b.(1+b) +ρ.(1−λ))
- = sum(: j.(j−1).f(j).power(j−2,ρ) ←j :).ρ.ρ +sum(: j.f(j).power(j−1,ρ) ←j :).(2−ρ).ρ −sum(: f(j).power(j,ρ) ←j :).(b.(1+b) +ρ.(1−λ))
- = sum(: j.(j−1).f(j).power(j,ρ) ←j :) +sum(: 2.j.f(j).power(j,ρ) ←j :) −sum(: j.f(j).power(j+1,ρ) ←j :) −sum(: b.(1+b).f(j).power(j,ρ) ←j :) +sum(: (λ−1).f(j).power(j+1,ρ) ←j :)
- = sum(: ((j.(j−1) +2.j −b.(1+b)).f(j) +(λ −1 −(j−1)).f(j−1)).power(j,ρ) ←j :)

Each coefficient of a power of ρ must be zero, so this implies, for each j,

- (j.(j+1) −(1+b).b).f(j) = (j −λ).f(j−1)

The left-hand side of this is zero at j = b and at j = −(1+b), implying f(j−1) = 0 in each case unless j = λ. We know λ > 0, so j−λ is non-zero for all j ≤ −(1+b) and we can infer that f(j) = 0 for all j < −(1+b). For large j, greater than either b or λ, if f(j−1) is non-zero we have f(j) asymptotically equal to f(j−1)/j making F equal to exp(ρ) plus some function whose power series has high-order coefficients tending to zero. This would make R behave as exp(ρ/2) for large ρ, which would not be normalizable. We can thus infer that f(j) must be zero for all sufficiently large j; implying that λ is an integer, with f(j) = 0 for j ≥ λ. We thus find possibilities;

- b ≥ λ &implies; f(j) non-zero for λ > j ≥ −(b+1), making F(ρ) equal some polynomial in ρ divided by power(1+b, ρ);
- λ > b &implies; f(j) non-zero for λ > j ≥ b, making F(ρ) equal some polynomial in ρ times power(b, ρ).

Now, when we come to normalize ψ we'll be integrating ψ.*ψ over volume, which will imply integrating R(r).R(r).r.r with respect to r; this will yield an infinite result if, near r = 0, it is 1/r or worse; so R(r) can afford to be a polynomial times 1/r but not times 1/r/r. Hence F(ρ) can afford to be a polynomial divided by ρ, but if it is divided by any higher power we shan't be able to normalize ψ. This limits the first case above to b = 0; but this case requires b ≥ λ and we know λ > 0, so the first case can't arise. Only the second possibility above remains.

Thus λ is a positive integer > b; we can write F(ρ) = power(b, ρ).sum(: g(j).power(j, ρ) ←j :) with (:g|λ−b) = (: f(j+b) ←j :) so

- (j+1+b −λ).g(j)/g(j+1)
- = (j+1+b).(j+b+2) −b.(1+b)
- = j.(j+b+2) +(1+b).(j+2)
- = j.j +j.b +2.j +j +b.j +2 +2.b
- = j.j +2.j.b +3.j +2 +2.b
- = (j+1).(j+2+2.b)

or

- (2.(1+b)+j).(j+1).g(j+1) = (1+b+j−λ).g(j)

The resulting polynomial F has order λ but its behaviour near the origin will be dominated by its lowest-order term, of order b; consequently the solutions for F with b small are predominantly near the origin, while those with b close to λ are predominantly far from the origin. Since b is non-negative – and R incorporates a factor exp(−ρ/2) – we are assured that our over-all wave function, ψ, will be normalizable.

Our resulting R(r) = F(ρ).exp(−ρ/2), with ρ = K.r, is R(r) = F(K.r).exp(−K.r/2) with K = √(−8.m.E)/ℏ and when we normalize our wave function we'll integrate our spherical harmonic's squared modulus over the whole sphere and integrate R(r).R(r).r.r over all positive r. Squaring R yields the square of F scaled by exp(−K.r); each term is then a scalar times a power of K.r times this exponential. Now,

- integral(: power(j, K.r).r.r.exp(−K.r) ←r :{positives})
- = integral(: power(j+2, s).exp(−s) ←s :{positives})/K/K/K
- = Γ(j+3)/K/K/K
- = (j+2)!/K/K/K

so integration of the square of (: R(r).r ←r :{positives}) reduces to computing the polynomial square of F, multiplying each coefficient of it by the factorial of two plus its order, summing (i.e. evaluating the thus scaled polynomial at 1) and dividing the answer by the cube of K.

For each positive integer λ we then have (re-arranging the definition of λ) solutions with energy

- E = −(m/2).(Q/ℏ/λ)
^{2}

For each natural b < λ and each integer a with −b ≤ a ≤ b, we get a separate solution; there are 2.b+1 such a; so this energy level's total degeneracy is

- sum(: 2.b+1 ←b |λ)
- = λ +2.sum(λ)
- = λ +λ.(λ−1)
- = λ
^{2}

For each atomic number Z and each λ, b and a as above, let
Ψ(Z,λ,b,a) be a suitably normalized solution to Schrödinger's
equation for this system. I'll refer to the outputs of any given
Ψ(Z,λ,b) collectively as an orbital

; and to the
union of orbitals associated with outputs of any given Ψ(Z,λ) as
a shell

. Although each shell and orbital formally depends
on Z, I'll elide mention of this dependence when discussing the properties that
hold true for each Z; likewise, each orbital formally depends on λ as
well but I'll elide mention of this when the value of λ is clear from
context or when discussing properties that hold true for each λ (and
Z).

Since the electrons occupying these states have intrinsic spin, with two
possible states, they contribute a further factor of two to the degeneracy; and
no two electrons can be in the same state, so no more than
2.λ^{2} electrons can fit into the λ-shell. When such a
shell is filled, the total energy of the electrons in it would (if we could
ignore the mutual repulsion of the electrons) thus be

- 2.λ
^{2}.E = −m.(Q/ℏ)^{2}

independent of λ; each filled shell contributes (roughly) the same amount to the total binding energy of the atom. Note that the length scale, r/ρ, of a state varies in inverse proportion to √(−E) and hence is proportional to λ/Z, where Z is the atomic number of the nucleus, since Q is proportional to Z. Indeed, Q/ℏ is Z.α.c, where α is the fine structure constant, so we can write the above, and our length scale, as:

- E = −(m/2).(α.c.Z/λ)
^{2} - 2/K = ℏ/√(−2.m.E) = λ/Z/(m.α.c/ℏ)

in which −E's co-efficient of the square of Z/λ is
(m/2).(α.c)^{2} = 13.6 eV, the Rydberg constant
(c.h.R_{∞}, where R_{∞} is the inverse wavelength of
the photon emitted when a previously free electron falls into the Z=λ
orbital), and 2/K's coefficient of λ/Z is ℏ/(m.α.c) = 52.92
pico metres, the Bohr radius. Similarly, 2.λ^{2}.E's
co-efficient of the square of Z is just twice the Rydberg energy, known as the
Hartree energy.

The above strictly only describes the energy levels of a nucleus accompanied
by a single electron (and, formally, m should then be a reduced mass

obtained by summing the reciprocals of the masses of electron and nucleus, then
taking the reciprocal of the sum; though this will differ from the electron's
mass by only one part in some thousands). However, we can reasonably interpret
various aspects of the more complex case by consideration of the above. In
particular, it is reasonable to expect electron states kindred to the above,
each characterized by positive λ, natural b < λ and integer a
between −b and b; albeit the energy of such a state may be perturbed from
the above and may even depend on b and a as well as λ, though it is
reasonable to suppose that the dependence on λ will remain
dominant. Indeed, we may also reasonably expect b to be more significant than a
(since b plays a rôle in determining the radial distribution); this is why
I follow orthodoxy in grouping states into orbitals and shells.

When we consider a real atom, we can imagine filling up the available states
for electrons near it. Each added electron sees a partially filled set of
states, possibly breaking the spherical symmetry and, at least, reducing the
effective charge of the apparent nucleus. If, for every λ and b that
have any states filled, we've filled all states with that λ and b, we can
expect that the state will have spherical symmetry. In such a state, say with N
electrons orbiting a nucleus of atomic number Z, the potential's variation with
radius will be more complex than the simple Coulomb potential: at large radius,
the potential.4.π.ε_{0}/e/e will effectively be (Z−N)/r
but, further in, the effective value of N will be reduced, yielding a bigger
potential. This will tend to give states predominantly close to the origin
lower (i.e. more negative) energy than ones further from the origin; thus
orbitals with low b will be filled in preference to those with high b.

Indeed, in the actual periodic table of
the elements, the states in each orbital have roughly the same energy as one
another; within a given shell, the orbitals' energies increase with b; indeed,
they do so swiftly enough that each Ψ(Z,λ,1+n) has similar energy to
Ψ(Z,λ+n,0) for n = 1, 2 at least; it seems reasonable to suspect
similar will apply for larger n, but we only know about the first hundred or so
elements and this is nowhere near enough to see what happens for larger n. No
Ψ(Z,λ,4) can arise until λ = 5; assume that Ψ(Z,5,4) has
distinctly higher energy than Ψ(Z,5,3), which quite plainly does have
comparable energy to Ψ(Z,6,2) and Ψ(Z,7,0); these last orbitals comprise
2.(7+5+1) = 26 states and only start filling at _{87}Fr, so won't be
full until element 112, so this is the earliest Ψ(Z,5,4), the 5g orbital,
could possibly start filling; and, in all likelihood, Ψ(Z,7,1), the 7p
orbital, has lower energy than it, so shall make it wait until after element 118
to tell us its story.

However, it is not impossible to get some hint at how things may work out. If we compute the radial dependencies of the solutions above, we can use them to determine electron densities u.*u and hence, for each orbital:

- its contribution, when full, to reducing the potential and – for any given subset of the other orbitals, using the thus-adjusted V(p) after their reductions are taken into account –
- its integral of u(p).*u(p).V(p), which should be a good estimator of its actual energy.

By choosing, iteratively, the orbital with lowest actual energy per
electron in the orbital, then adding it to the set of orbitals being used to
adjust V, this approach should reveal the correct order of filling of orbitals.
By looking at the energy of individual states within an orbital, it may also be
possible to discern which orbitals will *begin* filling first, even where
other orbitals may finish before them – as, for example, in the sequence
noted above: each Ψ(Z,λ,0) fills (groups I and II of the periodic
table) before Ψ(Z,λ−1,2) begins filling; but no sooner has each
Ψ(Z,λ−i,1+i) started filling than
Ψ(Z,λ−i−1,2+i) does so as well; the ones with higher i
fill before the ones with lower i, and borrow electrons off them as they go; as
indeed Ψ(Z,λ−1,2) borrows off Ψ(Z,λ,0) when finally
it comes to fill. Only once all Ψ(Z,λ−i,1+i) with 0 < i and
λ > 1+2.i have filled does Ψ(Z,λ,0) finally refill, after
which Ψ(Z,λ,1) fills, completing a cycle of the table, as groups III
through VIII.