. Schrödinger Equation

Paul Abbott

Pre-requisites

To improve the appearance of graphics, TraditionalForm will be used as the default format type for text that appears in graphics,

$FormatType = TraditionalForm;

and the default style for text in graphics is set as follows.

$TextStyle = {FontSize -> 10};

Several Mathematica graphics packages are used in this work. Entering

<< Graphics`

and for “pretty” plots:

SetOptions[ParametricPlot3D, 
        PlotPoints -> {50, 20}, PlotRange -> All, 
        Axes -> True, AxesLabel -> {x, y, z}, Ticks -> None, 
        Boxed -> True, BoxRatios -> {1, 1, 1}];

We define a routine for plotting angular momentum diagrams:

vecs(l_):=Show[Graphics[Table[{Text[(Sqrt[j (j + 1)]),  
                        0.85 {(Sqrt[l (l + 1)]), (Sqrt[l (l + 1)])}],GrayLevel[0.5],Line[(({{0, m}, {(Sqrt[l (l + 1) - (m^2)]), m}}))],GrayLevel[0],Arrow({0, 0}, {(Sqrt[l (l + 1) - (m^2)]), 
                          m}),Hue[1],Circle[{0, 
                          0},  (Sqrt[l (l + 1)]), {((-π)/2), (π/2)}]},{m, -l, l}],AspectRatio -> Automatic,Axes -> True,AxesLabel -> {"", m},Ticks -> {Table[i, {i, 0, l}], 
                  Table[i, {i, -l, l}]}]]

Complex conjugation is easily achieved using a replacement rule:

\(TraditionalForm\`\(exp_\^*\) := 
      exp /. Complex(a_, b_) -> Complex(a, \(-b\))\)

We use the Notation package,

Utilities`Notation`AutoLoadNotationPalette  =  False;

<< Utilities`Notation`

SetOptions[Notation, 
            WorkingForm  ->  TraditionalForm];

to define custom input and output notation for the 3-j and Clebsch-Gordan coefficients.

Notation[(Subscript[(({{j1_, j2_, j3_}, {m1_, m2_, m3_}}) ), 3  j]) ⟺ ThreeJSymbol[{j1_,  m1_},  {j2_,  m2_},  {j3_,  
                m3_}]]

Notation[⟨j1_,,j2_,,m1_,,m2_ \[VerticalSeparator] j_,,m_⟩ ⟺ ClebschGordan[{j1_,  m1_},  {j2_,  m2_},  {j_,  
                  m_}]]

Atomic orbitals are obtained by solving the Schrödinger equation,

H Ψ(r)=E Ψ(r), H=T+V(r),

where T and V are the kinetic and potential energy operators, respectively.

For a one-electron atom the potential is spherically symmetric and the solutions can be separated into radial and angular components via spherical polar coordinates r=(r, θ, ϕ). The complete solution, Ψ(r), is simply the product of these functions and depends upon three quantum numbers, n, l, and m, i.e.,

(Subscript[Ψ, n, l, m])(r)=(Subscript[R, n, l])(
            r) (Subscript[Y, l]^m)(θ, ϕ).

To assist you in understanding the concepts presented in this assignment, you may find the information at the Atomic Orbitals and Computational Chemistry hyperlinks useful.

... Laplacian operator

For a one-electron atom the kinetic energy operator (in atomic units) is T = -(1/2) Δ where Δ is the three-dimensional Laplacian operator,

\(TraditionalForm\`Δ = ∂\^2\/∂x\^2 + ∂\^2\/∂y\^2 + \(\(∂\^2\/∂z\^2\)\(·\)\)\)

In spherical polar coordinates, r=(r, θ, ϕ), this becomes

Δ=(1/(r^2))((∂((r^2) (∂/(∂r))))/(∂r))-(1/(r^2)) (l^2)(θ, ϕ),

where the angular part, (l^2)(θ, ϕ), is the square of the (total) angular momentum operator,

(l^2)(θ, ϕ)=-((1/(sin(θ)))((∂(sin(θ) (∂/(∂θ))))/(∂θ))+(1/((sin^2)(θ)))(((∂^2))/(∂(ϕ^2)))).

... Schrödinger equation for one-electron atoms

For a spherically symmetric potential, V(r)=V(r), using separation of variables the solution can be written as Ψ(r)=R(r) Y(θ, ϕ). Substituting this expression into the Schrödinger equation and choosing the separation constant to be l (l + 1) leads to a pair of equations: the angular equation,

\(TraditionalForm\`\(\(\(\(l\^2\)(θ, ϕ)\) \(Y(θ, ϕ)\) = 
      l\ \((l + 1)\) \(Y(θ, ϕ)\)\)\(,\)\)\)

and the radial equation,

(Subscript[R, n, l]^″)(
          r) + (2/r) (Subscript[R, n, l]^′)(
            r) + (2  E - 
              2  V(r) - ((l (l + 1))/(r^2))) (Subscript[R, n, l])(
            r) = 0.

... Angular solutions

Table of values

The angular eigenfunctions are the spherical harmonics, (Subscript[Y, l]^m)(θ, ϕ), which are built-in.

Eigenfunction equation

In Mathematica notation, the (squared) total angular momentum operator, ℒ, is

ℒ[f_]:=-((1/(sin(θ)))((∂(sin(θ) ((∂f)/(∂θ))))/(∂θ))+(1/((sin^2)(θ)))(((∂^2) f)/(∂(ϕ^2))))

and the z-component operator is

(Subscript[ℒ, z])[f_]:=-I ((∂f)/(∂ϕ))

The spherical harmonics satisfy the eigenfunction relations:

\(TraditionalForm\`\(\(\(\(l\^2\)(θ, ϕ)\) \(\(Y\_l\%m\)(θ, ϕ)\) = 
      l\ \((l + 
            1)\) \(\(Y\_l\%m\)(θ, ϕ)\)\)\(,\)\)\)

\(TraditionalForm\`\(\(l\_z\)(θ, ϕ)\) \(\(Y\_l\%m\)(θ, ϕ)\) = 
      m\ \(\(\(Y\_l\%m\)(θ, ϕ)\)\(.\)\)\)

Orthonormality

The usefulness of the spherical harmonics, (Subscript[Y, l]^m)(θ, ϕ), is that they are orthogonal (actually orthonormal) and form a complete basis set for the expansion of any angular function in 3-space (including the atomic orbitals of any atom). The orthonormality is defined by

(∫_0^(2  π))(∫_0^π)\(\(\(\(Y\_l\%m\)(θ, ϕ)\)\^*\) 
                    \(\(Y\_L\%M\)(θ, ϕ)\)\ \(sin(θ)\)\)dθdϕ=(Subscript[δ, l, L]) (Subscript[δ, m, M]),

where \(TraditionalForm\`\(\(\(Y\_l\%m\)(θ, ϕ)\)\^*\)\) is the complex conjugate of (Subscript[Y, l]^m)(θ, ϕ) (implemented as a Mathematica rule in the pre-requisites section above), and the sin(θ) term comes from the angular part of the integration volume element, dr=(r^2)drsin(θ)dθ\(dϕ\).

The general inner product can be defined using bra-ket notation as

⟨{l_, m_}, {L_, 
          M_}⟩ := (∫_0^(2 π))(∫_0^π) 
                  (Subscript[Y, l]^m)(θ, -ϕ) (Subscript[Y, L]^M)(θ, ϕ) sin(θ) \[DifferentialD]θ \[DifferentialD]ϕ

The spherical harmonics are eigenfunctions of the (squared) total angular momentum operator and therefore are a representation of a state of angular momentum, l, and projection quantum number, m.

Visualization

The (Subscript[Y, l]^m)(θ, ϕ) are, in general, complex. We therefore can either plot the real or imaginary parts or the absolute value. The angular charge density distribution

\(TraditionalForm\`\(\(\(Ω\_\(l, 
            m\)\)(θ) = |\(Y\_l\%m\)(θ, ϕ)|\^2\)\(,\)\)\)

is axially symmetric with respect to the azimuthal angle ϕ. After defining the angular distribution function:

\(TraditionalForm\`\(Ω\_\(l_, 
            m_\)\)(θ_) := \(\(Y\_l\%m\)(θ, ϕ)\)\ \(\(\(Y\_l\%m\)(θ, ϕ)\)\^*\)\)

individual charge densities can be immediately examined, e.g.,

(Subscript[Ω, 3, 1])(θ)

You can see that this result is independent of ϕ. In such cases, PolarPlot can provide a complete view of the atomic orbitals.

... Radial solutions

Assuming nuclear charge Z, the (Coulomb) potential is V = -(Z/r) and the radial equation reads

(Subscript[R, n, l]^″)(
          r) + (2/r) (Subscript[R, n, l]^′)(
            r) + (2  
                E + ((2  Z)/r) - ((l (l + 1))/(r^2))) (Subscript[R, n, 
                l])(r) = 0.

It turns out that this equation can be solved directly in terms of confluent hypergeometric functions.

Simplify[DSolve[R''(r)+((2 R'(r))/r)+(2  ℰ + ((2 Z)/r) - ((l (l + 1))/(r^2))) R(r)==0,R(r),r],l ≥ 0]

To examine bound states (i.e. ℰ < 0) we put ℰ -> -((Z^2)/(2 
                 (n^2))) which simplifies things slightly (note that we are not yet assuming n to be an integer).

Simplify[% /. ℰ -> -((Z^2)/(2  (n^2))),n > 0 ∧ Z > 0]

However, you might find it difficult to reconcile this result with that presented in your textbook so we now proceed to examine the solution using a general method for analysing differential equations.

General considerations

Consider the differential equation y''(x)+p(x) y'(x)+q(x) y(x)=0.

Re-scaled equation

To solve the radial differential equation for a spherically symmetric potential, let's follow the standard solution method. We would like to find a suitable re-scaling of the radial variable to simplify the radial equation which can be entered as

ℛ(R_)=(((\[DifferentialD]^2) R)/(\[DifferentialD](r^2)))+(2/r) ((\[DifferentialD]R)/(\[DifferentialD]r))+(2 ℰ + ((2 Z)/r) - ((l (l + 1))/(r^2))) R;

Total derivatives (e.g., ((dR)/(dr))) are used so that we can easily change variables.

If you have followed all the steps correctly, the re-scaling of the radial variable has changed the radial equation to the dimensionless form

(Subscript[R, n, l]^(′′))(ρ)+((2 (Subscript[R, n, l]^′)(ρ))/ρ)-(((l (l + 1))/(ρ^2)) - (n/ρ) + 
                  (1/4)) (Subscript[R, n, l])(ρ)=0,

which only has the two parameters (quantum numbers) n and l.

ℛ(R_)=(((\[DifferentialD]^2) 
                R)/(\[DifferentialD](ρ^2)))+(2/ρ) ((\[DifferentialD]R)/(\[DifferentialD]ρ))-(((l (l + 1))/(ρ^2)) - (n/ρ) + (1/4)) 
              R;

Behaviour as ρ -> ∞

In the limit ρ -> ∞, equation (..) becomes

R''(ρ)-(1/4) R(ρ)==0,

where the dependence on n and l has been suppressed.

Behaviour as ρ -> 0

(∫_0^∞)((((1/r) (E^(-r))))^2) (r^2) \[DifferentialD]r

For a solution to be acceptable, expectation values of the kinetic energy and potential energy operators must exist. Since the Coulomb potential energy operator is proportional to 1/r, it should be clear that a 1/r singularity in the wave function means that the potential expectation value would not be finite. A similar argument can be applied to the kinetic energy operator. (Consider what happens when Δ is applied to a 1/r singularity and then attempt to compute the expectation value of the kinetic energy operator.)

Factoring off the asymptotic behaviour

The physical solution, (Subscript[R, n, l])(ρ), to the original differential equation must have the correct behaviour for both ρ -> ∞ and ρ -> 0. One way to ensure this is to construct the solution in the form of a product \(TraditionalForm\`\(R\_\(n, 
              l\)\)(ρ) == \(𝒜(ρ)\)\ \(𝒫(ρ)\)\ \(\(L\_\(n . l\)\)(ρ)\)\) where (Subscript[L, n . l])(ρ) is to be determined. Note that there is no power series solution of equation (..) because (truncated) power series cannot display an exponential decrease (or increase) for large values of the expansion variable.  In fact, series always diverge to ±∞!  Hence you must factor off the exact asymptotic behaviour to be able to obtain a power series solution valid over a wide domain. After factoring off the large and small ρ behaviour, the resulting differential equation admits a simple series solution (which, for hydrogen, turns out to be a polynomial).

You should obtain (n - l - 1) L(ρ)+(2 l - ρ + 2) L'(ρ)+ρ L''(ρ). Mathematica can solve this differential equation — again in terms of hypergeometric functions:

DSolve[(n - l - 1)  L(ρ)+(2 l - ρ + 2) L'(ρ)+ρ L''(ρ)==0,L(ρ),ρ]

Instead, we will solve this differential equation using power series.

Power Series Solution

To solve differential equations such as x y''(x)+y'(x)+y(x) ==0 using series we substitute in the generic term of the series solution y(
          x) = (∑_(i = 0)^∞) (Subscript[a, i]) (x^i), i.e.,

x y''(x)+y'(x)+y(x) ==0/.y -> Function[x, (Subscript[a, i]) (x^i)]//ExpandAll

Next we use pattern matching to find the recurrence relation that must be satisfied by the (Subscript[a, i]).

% /.  
        c_ (x^(i + k_. )) :> (c /. i -> i - k - 1) // 
      Simplify

Now we can solve the recurrence relation to find (Subscript[a, i]) in terms of (Subscript[a, i - 1]).

Solve[%, (Subscript[a, i])]

For finiteness, this series must terminate and this requires n to be an integer (i.e., the states are quantized). If the series did not terminate, it can be shown that the recurrence relation implies that the solution would blow up exponentially for large ρ.

It turns out that the recurrence produces a special function known as the associated Laguerre polynomials (which are special cases of Hypergeometric1F1 functions), (Subscript[L, n - l - 1]^(2 l + 1))(ρ).

General formula

Combining the above results, since ρ = ((2  Zr)/n) you should be able to confirm that the (un-normalised) radial orbitals are given by (for n - l - 1 ≥ 0)

(Subscript[R, n_, l_])(r_) = 
        (r^l) (E^(-((Z r)/n))) (Subscript[L, n - l - 
                  1]^(2 l + 1))(((2 Z r)/n));

From the quantization condition, the energy of each state is (Subscript[ℰ, n]) = ((-(Z^2))/(2  (n^2))), and the radial equation (..) becomes

(Subscript[R, n, l]^(′′))(r)+(2/r) (Subscript[R, n, l]^′)(r)+(-((l(l + 1))/(r^2)) - 
                ((Z^2)/(n^2)) + ((2  Z)/r)) (Subscript[R, n, l])(r)=0.

In Mathematica notation one can write

(Subscript[𝒟, n_, l_])[R_]=(((\[DifferentialD]^2) R)/(\[DifferentialD](r^2)))+(2/r) ((\[DifferentialD]R)/(\[DifferentialD]r))+(-((l (l + 1))/(r^2)) - 
                  ((Z^2)/(n^2)) + ((2  Z)/r)) R;

Normalization

The radial wavefunctions given above are not normalised. To achieve this, one requires that

(∫_0^∞)(((Subscript[R, n, l])(r))^2)(r^2)dr==1,

where the (r^2) dr term comes from the integration volume element.

A closed form expression can be found for the normalization integral. It turns out to be \(TraditionalForm\`\(\((\(2\ Z\)\/n)\)\^\(l + 
              3/2\)\) \@\(\(\((n - l - 1)\)!\)\/\(2  n\ \(\((n + l)\)!\)\)\)\).

You are now in a position to define the normalized radial wavefunctions for n > l ≥ 0.

Energy orthogonality

Consider (∫_0^∞)(Subscript[R, n, l])(r)(Subscript[R, p, l])(r)(r^2)dr which involves two states with different principal quantum numbers (i.e., n ≠ p which implies the states have different energy).

Other integrals

Consider (∫_0^∞)(Subscript[R, n, l])(r)(Subscript[R, n, L])(r)(r^2)dr which  involves two radial states of different orbital angular momentum (i.e., l ≠ L).

Visualization

It is physically more relevant to examine the radial probability density which is given by

P(r)=(r^2)(((Subscript[R, n, l])(r))^2).

After normalisation, this will integrate to unity as required for a probability density.

The nodes of this function are places with zero probability of finding an electron. Since the radial functions depend only on the radius, r, this leads to a number of concentric electron "shells".

.. Stark effect for hydrogen

Consider a hydrogen atom in a weak uniform external electric field (Subscript[E, 0])(not related to the ground-state energy) which is in the z direction, i.e., the Hamiltonian perturbation is given by (H^′) = (Subscript[E, 0]) z. Applying first-order perturbation theory, the associated energy shift between two states with the same principal quantum number, n, is given by

\(TraditionalForm\`\(\(Δ\ E = \(\(E\_0\) ⟨n\ L\ M |z| 
              n\ l\ m⟩ = \(E\_0\) ⟨n\ L\ M |r\ \(cos(θ)\)| 
              n\ l\ m⟩\)\)\(,\)\)\)

where

| 
              n l m⟩≡(Subscript[R, n, l])(
                r) (Subscript[Y, l]^m)(θ, ϕ),

and the (Subscript[R, n, l])(r) are normalised. Explicitly, this requires the computation of

\(⟨n\ L\ M |r\ \(cos(θ)\)| 
                n\ l\ m⟩ ≡ ⟨n\ L |r| 
                  n\ l⟩ ⟨L\ M |cos(θ)| 
                  l\ m⟩\)=(∫_0^∞)(Subscript[R, n, L])(r) (Subscript[R, n, l])(r)(r^3)dr×(∫_0^(2  π))(∫_0^π)\(\(\(\(Y\_L\%M\)(θ, ϕ)\)\^*\)\ \(\(Y\_l\%m\)(θ, ϕ)\)\ \(sin(θ)\)\ \(cos(θ)\)\)dθ\(dϕ\),

where the separability of the radial and angular integrals is made evident. In this section you will use the expressions for the hydrogen wavefunction developed above to compute these integrals and find the associated energy shifts.

... Radial Integral

... Angular Integrals

⟨{L_, M_}, {l_, m_}⟩:=⟨{L, M}, {l, m}⟩=(∫_0^(2  π))(∫_0^π)(((-1))^M) (Subscript[Y, L]^(-M))(θ, ϕ) (Subscript[Y, l]^m)(θ, ϕ) sin(θ) cos(θ)\[DifferentialD]θ\[DifferentialD]ϕ

Table[⟨{L, M}, {l, 
          m}⟩, {L, 0, 2}, {l, 0, 2}, {M, -L, 
        L}, {m, -l, l}]

... An aside — Clebsch-Gordan coefficients

Introducing the Clebsch-Gordan coefficients, ⟨(Subscript[l, 1]) (Subscript[l, 2]) (Subscript[m, 1]) 
            (Subscript[m, 2]) | L M⟩ (there are many notations for these coefficients — see, e.g., Appendix 4 of Bransden and Joachain, Physics of Atoms and Molecules (Longman, 1983)), and the closely related (more symmetric) 3-j symbols, (({{(Subscript[j, 1]), (Subscript[j, 2]), (Subscript[j, 3])}, {(Subscript[m, 1]), (Subscript[m, 2]), (Subscript[m, 3])}})),

⟨(Subscript[l, 1]) (Subscript[l, 2]) (Subscript[m, 1]) 
                (Subscript[m, 2]) | L M⟩=(((-1))^(M + (Subscript[l, 1]) - (Subscript[l, 2]))) (Sqrt[2 L + 1]) (({{(Subscript[l, 1]), (Subscript[l, 2]), L}, {(Subscript[m, 1]), (Subscript[m, 2]), -M}})),

which vanish unless (Subscript[m, 1]) + (Subscript[m, 2]) ==  M (physical condition for addition of projection quantum number) and |(Subscript[l, 1]) - 
            (Subscript[l, 2])| ≤ L ≤ 
        (Subscript[l, 1]) + (Subscript[l, 2]) (triagularity condition, i.e., imagine the angular momenta (Subscript[l, 1]), (Subscript[l, 2]), and L as the (integer) side-lengths of a triangle), the angular integrals can be computed directly using the following general identity:

(∫_0^(2  π))(∫_0^π)\(\(\(Y\_\(l\_1\)\%\(m\_1\)\)(θ, ϕ)\)\ \(\(Y\_\(l\_2\)\%\(m\_2\)\)(θ, ϕ)\)\ \(\(Y\_\(l\_3\)\%\(m\_3\)\)(θ, ϕ)\)\ \(sin(θ)\)\)dθ\(dϕ\)=
	(Sqrt[(((2  (Subscript[l, 1]) + 1) (2  (Subscript[l, 2]) + 1) (2  (Subscript[l, 3]) + 
                        1))/(4  π))])(({{(Subscript[l, 1]), (Subscript[l, 2]), (Subscript[l, 3])}, {0, 0, 0}}))(({{(Subscript[l, 1]), (Subscript[l, 2]), (Subscript[l, 3])}, {(Subscript[m, 1]), (Subscript[m, 2]), (Subscript[m, 3])}})).

Since (Subscript[Y, 1]^0)(θ, ϕ)=(1/2) (Sqrt[(3/π)]) cos(θ), and \(\(\(Y\_l\%m\)(θ, ϕ)\)\^*\)=(Subscript[Y, l]^m)(θ, -ϕ)=\(\((\(-1\))\)\^m\ \(\(Y\_l\%\(-m\)\)(θ, ϕ)\)\), we find that

(∫_0^(2  π))(∫_0^π)\(\(\(\(Y\_L\%M\)(θ, ϕ)\)\^*\)\ \(\(Y\_l\%m\)(θ, ϕ)\)\ \(sin(θ)\)\ \(cos(θ)\)\)dθ\(dϕ\)=(((-1))^m)(Subscript[δ, M, 
              m])(Sqrt[(2  L + 1) (2  l + 1)])(({{L, l, 1}, {0, 0, 0}}))(({{L, l, 1}, {-M, m, 0}})).

Since the 3-j symbol (({{(Subscript[j, 1]), (Subscript[j, 2]), (Subscript[j, 3])}, {(Subscript[m, 1]), (Subscript[m, 2]), (Subscript[m, 3])}})) is built-in to Mathematica as ThreeJSymbol, computing these integrals can be sped up.

Here is a palette for entering this special notation.

(Subscript[(({{\[Placeholder], \[Placeholder], \[Placeholder]}, {\[Placeholder], \[Placeholder], \[Placeholder]}}) ), 3  j])

... Total Matrix Element

n = 2

Fix attention to the case n = 2. The possible {n, l, m} triples (kets) are

kets=(({{2, 0, 0}, {2, 1, -1}, {2, 1, 0}, {2, 1, 1}}));

and the perturbation matrix — computed from the matrix element of each pair of kets using Outer — is simply

Outer[AngleBracket, kets, kets, 1]

The (matrix) eigenvalues of the perturbation matrix give the energy shifts.

{λ, u} = Eigensystem[%]

λ

The (matrix) eigenvectors give the appropriate linear combinations of the initial states that correspond to these energy shifts.

u . Apply[Ket, kets, 1] // TableForm

n = 3

.. Helium perturbation calculation

The Schrödinger equation (..) for a two-electron atom (with nuclear charge Z) reads

ℋ Ψ(r)=(Underscript[\(\(-\(1\/2\)\) \((Δ\_1 + Δ\_2)\)\), (Underscript[⏟, 𝒯])])Ψ(r)+(Underscript[(-(Z/(Subscript[r, 1])) - (Z/(Subscript[r, 2])) + (1/(Subscript[r, 12]))), (Underscript[⏟, 𝒱])])Ψ(r)=ℰ Ψ(r),

where ℋ is the Hamiltonian, 𝒯 is the kinetic energy operator, and 𝒱 is the potential energy. The interparticle coordinates — (Subscript[r, 1]), (Subscript[r, 2]), and (Subscript[r, 12]) — are defined below:

Show[
        Graphics[{
Arrow({0, 0}, {1, 0}), 
            Arrow({1, 0}, {1.3, 1}), 
            Arrow({0, 0}, {1.3, 1}), 
Text[
              (Subscript[r, 1]),  {1/2, -0.1}], Text[(Subscript[r, 12]),  {1.2, 0.5}], 
            Text[(Subscript[r, 2]),  {0.7, 0.63}]}]];

(Subscript[Δ, i]) is the three-dimensional Laplacian operator for electron i,

(Subscript[Δ, i])=\(1\/r\_i\%2\)((∂((Subscript[r, i]^2) (∂/(∂(Subscript[r, i])))))/(∂(Subscript[r, i])))-\(\(1\/r\_i\%2\) \(\(\(l\_i\%2\)(θ\_i, ϕ\_i)\)\(.\)\)\)

Considering the inter-electron potential, 1/(Subscript[r, 12]), as a perturbation, the Hamiltonian  can be written as

ℋ=(Underscript[\(\(-\(1\/2\)\) \((Δ\_1 + Δ\_2)\) - Z\/r\_1 - Z\/r\_2\), (Underscript[⏟, (ℋ^((0)))])])+(Underscript[(1/(Subscript[r, 12])), (Underscript[⏟, ℋ'])]),

where

\(TraditionalForm\`\(\(ℋ\^\((0)\) = \(-\(1\/2\)\) \((Δ\_1 + Δ\_2)\) - Z\/r\_1 - 
          Z\/r\_2 ≡ \(\(-\(1\/2\)\) Δ\_1 - Z\/r\_1\)\+\(⏟\+\(ℋ\_1\)\) - \(1\/2\) Δ\_2 - \(Z\/r\_2\)\+\(⏟\+\(ℋ\_2\)\)\)\(,\)\)\)

and

\(TraditionalForm\`ℋ\^′ = \(\(1\/r\_12\)\(·\)\)\)

To zeroth order, the Hamiltonian is simply separated into the sum of two single-particle (hydrogenic) operators, \(TraditionalForm\`ℋ\^\((0)\) = ℋ\_1 + ℋ\_2\). Hence the (unnormalised) ground state of the a two-electron atom can be approximated by a product of hydrogenic functions (with nuclear charge Z),

\(TraditionalForm\`\(\(ψ\^\((0)\) = e\^\(\(-Z\)\ \((r\_1 + r\_\(\(2\)\(\ \)\))\)\)\)\(,\)\)\)

and the energy is the sum of the two hydrogenic (ground-state) energies, i.e.,

(ℋ^((0))) (ψ^((0)))=(ℰ^((0)))(ψ^((0))) where (ℰ^((0)))=-((Z^2)/2) - ((Z^2)/2) = -(Z^2).

The first-order perturbation energy, (ℰ^((1))), is \(TraditionalForm\`ℰ\^\((1)\) = ℰ\^\((0)\) + Δ ℰ\), where the energy shift is given by

\(TraditionalForm\`Δ ℰ = \(\(⟨\(ψ\^\((0)\)\) |1\/r\_12| ψ\^\((0)\)⟩\/⟨ψ\^\((0)\) | ψ\^\((0)\)⟩\)\(.\)\)\)

In interparticle coordinates  the matrix element of a general operator, 𝔒, can be computed from the integral

\(TraditionalForm\`⟨ψ |𝔒| ψ⟩ = 
      8 \( π\^2\) \(∫\_0\%∞\( 
                r\_1\) dr\_1 \(∫\_0\%∞\( 
                    r\_2\) dr\_2 \(∫\_\(|r\_1 - r\_2|\)\%\(r\_1 + r\_2\)\(ψ\^*\) 𝔒\ ψ\ \(r\_12\) \(\(dr\_12\)\(.\)\)\)\)\)\)

In the following sections we will compute the first-order perturbation energy from equations (..) and (..).

... A generating integral

In this section you will verify that

\(TraditionalForm\`\(\(⟨ϕ |1\/\(\(r\_1\) \(r\_2\) r\_12\)| ϕ⟩ = \(2  π\^2\)\/\(\((λ\_1 + λ\_2)\) \((λ\_1 + λ\_12)\) \((λ\_2 + λ\_12)\)\)\)\(,\)\)\)

where

ϕ(
        r) = (e^(-(Subscript[λ, 1]) (Subscript[r, 1]) - (Subscript[λ, 2]) 
                (Subscript[r, 2]) - (Subscript[λ, 12]) 
                (Subscript[r, 12]))).

The importance of this integral is that matrix elements of any combination of integral powers (≥-1) of the interparticle coordinates can be obtained by differentiation with respect to the parameters, (Subscript[λ, i]).  Moreover, this general integral permits the examination of a whole class of zeroth order approximations (including exponential terms of the above form multiplied by polynomials in the interparticle coordinates).

(Subscript[r, 1]) >  (Subscript[r, 2])

The integral for the case (Subscript[r, 1]) > (Subscript[r, 2]) can be computed through the steps

SetOptions[Integrate, 
      GenerateConditions -> False]

8 (π^2) (∫_((Subscript[r, 1]) - (Subscript[r, 2]))^((Subscript[r, 1]) + (Subscript[r, 2])))(E^(-2 (Subscript[r, 1]) (Subscript[λ, 1]) - 
                2 (Subscript[λ, 2]) (Subscript[r, 2]) - 
                2 (Subscript[λ, 12]) (Subscript[r, 12]))) \[DifferentialD](Subscript[r, 12])

(∫_0^((Subscript[r, 1])))% \[DifferentialD](Subscript[r, 2]) // 
      Simplify

and

(∫_0^∞) % \[DifferentialD](Subscript[r, 1]) // 
      Simplify

Under the conditions on λ we have

(Subscript[ℐ, (Subscript[r, 1]) > (Subscript[r, 2])]) = ((2 (π^2))/(((Subscript[λ, 1]) + (Subscript[λ, 2])) ((Subscript[λ, 1]) + (Subscript[λ, 12])) ((Subscript[λ, 1]) + (Subscript[λ, 2]) + 2 (Subscript[λ, 12]))));

(Subscript[r, 1]) <  (Subscript[r, 2])

Total integral

General matrix elements

Since \(TraditionalForm\`ϕ(
          r) = e\^\(\(-λ\_1\) r\_1 - \(λ\_2\) r\_2 - \(λ\_12\) r\_12\)\),

\(TraditionalForm\`∂\_\(λ\_i\)⟨ϕ |1\/\(\(r\_1\) \(r\_2\) r\_12\)| ϕ⟩ = \(-2\) ⟨ϕ |r\_i\/\(\(r\_1\) \(r\_2\) r\_12\)| ϕ⟩, \ i = 1, 2, 12. \)

General relation

From this it can be shown that

\(TraditionalForm\`\(∂\^\(i + j + m + 3\)\/\(∂λ\_1\%\(i + 1\) ∂λ\_2\%\(j + 1\) ∂λ\_12\%\(m + 1\)\)\) ⟨ϕ |1\/\(\(r\_1\) \(r\_2\) r\_12\)| ϕ⟩ = \(\((\(-2\))\)\^\(i + j + m + 
              3\)\) \(\(⟨ϕ |\(r\_1\%i\) \(r\_2\%j\) 
                  r\_12\%m| ϕ⟩\)\(.\)\(\ \)\)\)

Define this general matrix element (after moving the power of 2 to the left-hand side of the equation) as

⟨i_, j_, 
        m_⟩ := ⟨i, j, 
          m⟩ = 
        Simplify[(1/(((-2))^(i + j + m + 3))) (∂_({(Subscript[λ, 1]), i + 1}, {(Subscript[λ, 2]), j + 1}, {(Subscript[λ, 12]), m + 1}))ℐ]

... Normalisation

... Electron-electron interaction

... Total energy