Computational Physics Exam

Instructions

This exam, worth 40% of the total mark for Computational Physics, is due in by 4 pm on Thursday May 13.  Each student is to do the exam by themselves and submit their solutions individually.

The exam consists of 6 questions of approximately equal difficulty and several hints as to how to proceed. You should attempt to answer  whole questions whose total value is 100 marks.  The marks for each question is clearly indicated in red brackets (i.e., [n]) and, obviously, the marks provide a good guide to the question's difficulty.  Bonus marks will be given for elegant or unconventional solutions.

. Convergence Acceleration [20]

It is very common in mathematical physics and engineering to expand functions or solve differential equations by expansion into a Fourier basis or into a system of orthogonal polynomials.  However, because these basis functions are well-behaved, convergence near singular points or discontinuities is slow — a famous example being the Gibbs' phenomena for Fourier series. If the form of the singularity is known then one can explicitly take this into account using series re-arrangement to dramatically improve the convergence of series expansions.

As a specific example, consider the generating function for the Legendre polynomials, (Subscript[P, n])(x),

(1/(Sqrt[1 - 2  x t + (t^2)]))= (∑_(n = 0)^∞)(Subscript[P, n])(x)(t^n), -1 < x < 1,|t|<1.

(1/(Sqrt[1 - 2  x t + (t^2)]))- (∑_(n = 0)^5)(Subscript[P, n])(x)(t^n)+O[t]^6

R = (Sqrt[1 - 2  x t + (t^2)]);

Factor[(∂_t) R]==((t - x)/R)

Simplify[(((t^2) - 2 x t + 1)/R)] == R

Hence we have that

(∂_t) R==(∂_t)(((R^2)/R))==(∂_t)((((t^2) - 2 x t + 1)/R))==((t - x)/R) ⟹(∂_t)(((t^2) - 2 x t + 1)(∑_(n = 0)^∞)(Subscript[P, n])(x)(t^n))== (t - x)(∑_(n = 0)^∞)(Subscript[P, n])(x)(t^n),

since

(1/R)== (∑_(n = 0)^∞)(Subscript[P, n])(x)(t^n).

Term-by-term,

(∂_t)(((t^2) - 2 x t + 1)(Subscript[P, n])(x)(t^n))==(t - x)(Subscript[P, n])(x)(t^n)

x (Subscript[P, n])(x)=((n (Subscript[P, n - 1])(x)+(n + 1) (Subscript[P, n + 1])(x))/(2 n + 1))·

Since (∂_t)((1 - 2  x t + (t^2))(Subscript[P, n])(x)(t^n))==(t - x)(Subscript[P, n])(x)(t^n),

(∂_t)((1 - 2  x t + (t^2))(Subscript[P, n])(x)(t^n))==(t - x)(Subscript[P, n])(x)(t^n)//ExpandAll

Using pattern-matching, to re-arrange the (implied) summations,

\(TraditionalForm
    \`% /. a_\ t\^\(n + k_. \) :> 
        \((a /. n -> n - k)\) t\^n\)

we find that x (Subscript[P, n])(x) is equal to

( x (Subscript[P, n])(x)/.First[Solve[%, (Subscript[P, n])(x)]]//Simplify)

Hence, the recurrence relation can be implemented as a general replacement rule:

ℛ=(x^m)_.  (Subscript[P, n_])(x)->Expand[(x^(m - 1)) %]

Table[ExpandAll[x (Subscript[P, n])(x)== ((n (Subscript[P, n - 1])(x)+(n + 1) (Subscript[P, n + 1])(x))/(2 n + 1))],{n, 0, 5}]

Collect[(x^3) (Subscript[P, n])(x)//.ℛ,{(Subscript[P, n_])(x)},Factor]

(1 - x)(∑_(n = 0)^∞)(Subscript[a, n]) (Subscript[P, n])(x)≡  (∑_(n = 0)^∞)(-((n (Subscript[a, n - 1]))/(2 n - 1)) + (Subscript[a, n]) - 
                (((n + 1) (Subscript[a, n + 1]))/(2 n + 3))) (Subscript[P, n])(x).

Expand[(Subscript[a, n]) (1 - x)(Subscript[P, n])(x)]/.ℛ

Re-indexing the summation using pattern-matching we obtain

Collect[Expand[%]/.anything_ (Subscript[P, n + k_. ])(x):>(anything /. n -> n - k)(Subscript[P, n])(x) ,{(Subscript[P, n])(x),\(a\__\)},Factor]

(1/(Sqrt[1 - x]))=(Sqrt[2]) (∑_(n = 0)^∞)(Subscript[P, n])(x).

pl=Plot[Evaluate({(1/(Sqrt[1 - x])),(Sqrt[2]) (∑_(n = 0)^40)(Subscript[P, n])(x)}),{x, -1, 1},PlotStyle -> {Hue[1], Hue[(1/3)]},PlotRange -> {0, 4}];

The sum is slowly convergent.

Multiplying the numerator and denominator of the right hand side of equation (.) by (1 - x) (to remove the (Sqrt[1 - x]) singularity),

(1/(Sqrt[1 - x]))=(1/(1 - x))((1 - x)(Sqrt[2]) (∑_(n = 0)^∞)(Subscript[P, n])(x)),

and, from equation (.), we compute

\(TraditionalForm
    \`\(-\(\(n\ a\_\(n - 1\)\)\/\(2\ n - 1\)\)\) + a\_n - 
          \(\((n + 1)\)\ a\_\(n + 1\)\)\/\(2\ n + 3\) /. a\_n_ -> \@2 // 
      Factor\)

to obtain

(1/(Sqrt[1 - x]))=((-2 (Sqrt[2]))/(1 - x)) (∑_(n = 0)^∞)(((Subscript[P, n])(x))/((2 n + 3) (2 n - 1)))·

Here are the comparative plots:

Plot[-((2 (Sqrt[2]))/(1 - x)) (∑_(n = 0)^4)(((Subscript[P, n])(x))/((2 n + 3) (2 n - 1))),{x, -1, 1},PlotStyle -> Hue[(2/3)],Epilog -> pl[[1]],PlotRange -> PlotRange[pl]];

The convergence acceleration is quite dramatic. We have used only 4 terms in the accelerated sum compared to 40 in the original expression!  One reason for the acceleration is obvious: the coefficients of the (Subscript[P, n])(x) in the accelerated sum go asymptotically like 1/(n^2). The technique outlined here is easily generalized to other eigenfunction expansions.

. Maxwell Construction [20]

The equation of state of an ideal gas of 𝒩 noninteracting classical particles reads

p V == 𝒩 (Subscript[k, B]) T,

where (Subscript[k, B]) is the Boltzmann constant. A simple approximation, taking into account the interaction between particles, is the van der Waals equation of state,

\(TraditionalForm
    \`\(\(\((p + a\/V\^2)\) \((V - b)\) == 
      𝒩\ \(k\_B\) T\)\(,\)\)\)

with parameters a and b.  As you will see, there exists a critical temperature (Subscript[T, c]) below which loops appear in the p - V diagram.  These loops are forbidden for thermodynamical reasons. For T < (Subscript[T, c]), a phase transition between gas and liquid takes place; at T = (Subscript[T, c]) the difference between the liquid and gaseous phases vanishes.

(Subscript[T, c]) and (Subscript[V, c]) are determined by the condition that the first and second derivatives of p(V, T) with respect to V vanish, i.e., ((∂p(V, T))/(∂V))==0 and (((∂^2) p(V, T))/(∂(V^2)))==0.

Define p(V, T),

\(TraditionalForm
    \`p(V_, T_) = 
      \(p /. Solve[
              \((p + a\/V\^2)\) \((V - b)\) == 
                𝒩\ \(k\_B\) T, p] // First\) // Simplify\)

and solve for the critical values,

{(Subscript[V, c]), (Subscript[T, c])}={V, T}/.First[Solve[{((∂p(V, T))/(∂V)),(((∂^2) p(V, T))/(∂(V^2)))}=={0, 0},{V, T}]]

The critical pressure is then,

(Subscript[p, c]) = p((Subscript[V, c]), (Subscript[T, c]))

After re-scaling,

\(TraditionalForm
    \`\(p\&~\)(v_, t_) = Apart[\(p(v\ V\_c, t\ T\_c)\)\/p\_c, v]\)

we obtain the dimensionless form of the van der Waals equation:

\(TraditionalForm
    \`\(Solve[p\&~ == \(p\&~\)(V\&~, T\&~), T\&~] // First\) // Factor
      \)

\(TraditionalForm
    \`\(vanderWaals = 
        Plot[Evaluate[Table[\(p\&~\)(v, t), {t, 0.85, 1.05, 0.05}]], 
          {v, 0.5, 5}];\)\)

\(TraditionalForm
    \`∫\_v\_1\%v\_3\(\( p\&~\)(v)\) dv = 
      \((v\_3 - v\_1)\) \(\(\(p\_t\)\&~\)\(·\)\)\)

({{\(\(p\&~\) \((v\_1, t)\)\), ==, \(\(\(p\&~\)(v\_3, t)\)\(,
            \)\)}, {\(∫\_v\_1\%v\_3\(\( p\&~\)(v, t)\) dv\), ==, \(\((v\_3 - v\_1)\) \(\(\(p\&~\)(v\_1, t)\)\(.\)\)\)}})

The two equations are

\(TraditionalForm
    \`\(ℰ\_1\)(t_) = 
      \(p\&~\)(v\_1, t) == \(p\&~\)(v\_3, t)\)

and

\(TraditionalForm
    \`\(ℰ\_2\)(t_) = 
      \((\(∫\_a\%b\(\( p\&~\)(v, t)\) \[DifferentialD]v /. 
                {a -> v\_1, b -> v\_3} // Expand\) // Simplify)\) 
        == \(\(p\&~\)(v\_1, t)\) \((v\_3 - v\_1)\)\)

An automatic way of finding initial guesses is as follows:

initial[t_ /; t < 1]:=Module[{v, v1, v3},{v1, v3}=v/.\(Solve[∂\_v\(\( p\&~\)(v, t)\) == 0, v]
                        \)//Rest;\(\((v /. Solve[
                      \(p\&~\)(v, t) == \(p\&~\)(\(v1 + v3\)\/2, t), 
                      v])\)[[{1, 3}]]\)]

These initial guesses allow the volumes to be determined using FindRoot:

\(TraditionalForm
    \`volumes[t_ /; t < 1] := 
      Module[{v1, v3}, {v1, v3} = initial[t]; 
        FindRoot(
          Evaluate(
            {\(ℰ\_1\)(t), \(ℰ\_2\)(t)}, 
            {v\_1, v1}, {v\_3, v3}))]\)

Here is (Subscript[v, 1]) and (Subscript[v, 3]) for each temperature, \(T\&~\)=0.85`,0.9`,0.95:

Table(volumes[t], {t, 0.85, 0.95, 0.05})

Defining

\(TraditionalForm\`\(p\&~\_Maxwell\)(v_, t_) := \(p\&~\)(v, t)\)

and

\(TraditionalForm
    \`\(p\&~\_Maxwell\)(V_, t_ /; t < 1) := 
      With[{vols = volumes(t)}, 
        If[V < v\_1 ∨ V > v\_3 /. vols, \(p\&~\)(V, t), 
          \(p\&~\)(v\_1 /. vols, t)]]\)

here are plots of the Maxwell construction:

\(TraditionalForm
    \`\(Plot(Evaluate(Table(\(p\&~\_Maxwell\)(v, t), {t, 0.85, 1.05, 0.05})), 
        {v, 0.5, 5}, PlotRange -> {0, 2}, 
        AxesLabel -> {v\&~, p\&~}, 
        PlotStyle -> Table(Hue(i\/5), {i, 5}), 
        Epilog -> vanderWaals[[1]]);
      \)\)

The three lowest curves show the separation into a liquid (small volume) and gaseous phase. Along each two-phase line, liquid and gas co-exist in thermal equilibrium.

. Normal Modes of a Circular Membrane [20]

If there is no external force, the (wave) equation for free vibrations of a uniform isotropic membrane reads

(((∂^2) u)/(∂(t^2)))==(c^2) Δ u,

where Δ is the two-dimensional Laplacian. For a circular membrane of radius 1 fixed at its edge, polar coordinates (r, θ) are convenient; the boundary condition is u(1, θ, t) == 0 and the Laplacian reads

\(TraditionalForm
    \`Δ[f_] := 
      \(1\/r\) ∂\((r ∂f\/∂r)\)\/∂r + 
        \(1\/r\^2\) ∂\^2 f\/∂θ\^2\)

Substituting the separation ansätz,

\(TraditionalForm
    \`\(u(r_, θ_, t_) := 
        \(R(r)\)\ \(Θ(
            θ)\)\ E\^\(I\ ω\ t\);
      \)\)

into equation (.) we obtain

(separation=((r^2)/((c^2) u(r, θ, t)))((((∂^2) u(r, θ, t))/(∂(t^2)))-(c^2) Δ(u(r, θ, t)))//Expand)==0

From this we can conclude that

-(((ω^2) (r^2))/(c^2))-((R''(r) (r^2))/(R(r)))-((R'(r) r)/(R(r)))==-(n^2),

and

-((Θ''(θ))/(Θ(θ)))==(n^2),

where we have chosen the separation constant to be (n^2).

Using Mathematica, selecting the term independent of r, gives the equation for Θ:

\(TraditionalForm
    \`Θeqn = 
      Select[separation, FreeQ[#, r] &] == n\^2\)

and, by subtraction, the differential equation for R:

\(TraditionalForm
    \`Reqn = \((separation - First[Θeqn])\) == 
        \(-n\^2\)\)

The solution to the differential equation for Θ is immediate:

\(TraditionalForm
    \`DSolve(Θeqn, Θ(θ), θ)\)

We obviously require n to be an integer, i.e., n ∈ ℤ, for Θ(θ) to be single-valued.

The differential equation for R reads

Reqn=Collect[-((R(r))/(r^2)) (Subtract @@ Reqn),{R(r),R'(r),R''(r),r}]==0/.((ω^2)/(c^2)) -> (k^2)

DSolve[Reqn,R(r),r]

Since we require n ∈ ℤ, the solution obtained in [c] simplifies as follows.

Table[Simplify[%], {n, 0, 10}]

We see that each result only depends on one function, i.e., (Subscript[J, -n])(k r) and (Subscript[J, n])(k r) are dependent.

\(TraditionalForm
    \`\(Plot[Evaluate[Table[\(J\_n\)(z), {n, 0, 4}]], {z, 0, 10}, 
        PlotStyle -> Table[Hue[i\/5], {i, 5}]];\)\)

\(TraditionalForm
    \`\(Plot[Evaluate[Table[\(Y\_n\)(z), {n, 0, 4}]], {z, 0, 10}, 
        PlotStyle -> Table[Hue[i\/5], {i, 5}], 
        PlotRange -> PlotRange[%]];\)\)

A physical solution must have R finite at r = 0, so the coefficient of (Subscript[Y, n])(k r) must be 0 since (Subscript[Y, n])(0) is infinite (see the plots in [e]). The general solution reads

(Subscript[J, n])(k r)((Subscript[c, 2]) cos(n θ)+(Subscript[c, 1]) sin(n θ))(e^(i ω t)).

The boundary condition u(1, θ, t) = 0 implies that (Subscript[J, n])(k) = 0.  Hence, the allowed values of k are k = (Subscript[j, n, m]), where the (m^th) root of (Subscript[J, n])(z) is denoted (Subscript[j, n, m]).

<< NumericalMath`BesselZeros`

zeros = Table(BesselJZeros(n, 10), {n, 0, 10})

ordered = Sort(Flatten(zeros));

topten = Take(ordered, 10)

Determine the indices using Position:

\(TraditionalForm
    \`orderedpairs = 
      First /@ \((\(\((Position[zeros, \ #])\) &\) /@ topten)\)\)

Here are the corresponding radial functions, (Subscript[R, n, m])(r).

k /: (Subscript[k, n_, m_]) := zeros[[n + 1, m]]

\(TraditionalForm
    \`R /: R\_\(n_, m_\) := Function[r, \(J\_n\)(k\_\(n, m\)\ r)]\)

\(TraditionalForm
    \`Map[R\_\(First[#] - 1, Last[#]\)[r] &, orderedpairs]\)

\(TraditionalForm
    \`\(Plot[Evaluate[%], {r, 0, 1}, 
        PlotStyle -> Table(Hue(i\/Length[%]), {i, Length[%]})];\)\)

\(TraditionalForm
    \`ω /: ω\_\(n_, m_, c_:  1\) := c\ k\_\(n, m\)\)

\(TraditionalForm
    \`u /: u\_\(n_, m_, c_:  1\) := 
      Function[{r, ϕ, t}, 
        R\_\(n, m\)[r]\ \(cos(n\ ϕ)\)\ \(cos(ω\_\(n, m, c\)\ t)\)]
        \)

<< Graphics`ParametricPlot3D`

Table(CylindricalPlot3D(Evaluate[(Subscript[u, 2, 1])(r, ϕ, t)],{r, 0, 1, 0.1},{ϕ, 0, 2 π, (π/15)},Boxed -> False,BoxRatios -> {1, 1, 1},Axes -> None,PlotRange->(({{-1, 1}, {-1, 1}, {-1, 1}}))),{t, 0, ((2 π)/(Subscript[ω, 2, 1])), 
              ((2 π)/(8 (Subscript[ω, 2, 1])))});

. Ground State Energy of Helium [20]

In atomic units, the Hamiltonian for a two-electron atom with nuclear charge (SuperPlus[𝒵]) reads

\(TraditionalForm
    \`\(\(ℋ = 
      \(-\(1\/2\)\) \((Δ\_1 + Δ\_2)\) - 
        𝒵\/r\_1 - 𝒵\/r\_2 + 1\/r\_12\)\(,
    \)\)\)

where the interparticle coordinates —  (Subscript[r, 1]), (Subscript[r, 2]), (Subscript[r, 12])— are defined in Figure 1.

$FormatType = TraditionalForm;

$TextStyle = {FontSize -> 10};

\(TraditionalForm
    \`\(Show[Graphics[
          {Arrow[{0, 0}, \ {1, 0.1}]\ , Arrow[{0, 0}, \ {0.7, 1}], 
            Arrow[{0.7, 1}, {1, 0.1}], AbsolutePointSize[40], \ 
            Point[{0, 0}], AbsolutePointSize[20], Point[{1, 0.1}], 
            Point[{0.7, 1}], Text[r\_1, \ {0.3, 0.55}], 
            Text[r\_2, \ {0.5, 0}], Text[r\_12, \ {0.9, 0.55}], GrayLevel[1], 
            Text[1, \ {0.7, 1}], Text[2, \ {1, 0.1}], 
            Text[\(𝒵\^+\), \ {0, 0}]}], 
        AspectRatio -> Automatic, 
        PlotRange -> {{\(-0.1\), 1.1}, {\(-0.1\), 1.1}}];\)\)

In interparticle coordinates, the matrix element of a general operator, ℒ, is defined by the following 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 assignment 6 you showed that for

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

the matrix element of an arbitrary monomial, (Subscript[r, 1]^i) (Subscript[r, 2]^j) (Subscript[r, 12]^k) with i, j, k ≥  -1, can be computed by differentiating a particular matrix element, viz,

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

where

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

The importance of this result 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]).

In interparticle coordinates the (s-state) kinetic energy operator, \(TraditionalForm
      \`𝒯 = 
        \(-\(1\/2\)\) \((Δ\_1 + Δ\_2)\)\), reads

𝒯(ψ_):=-(1/2) ((((∂^2) ψ)/(∂((Subscript[r, 1])^2)))+(((∂^2) ψ)/(∂((Subscript[r, 2])^2)))+2(((∂^2) ψ)/(∂((Subscript[r, 12])^2)))+(2/(Subscript[r, 1]))((∂ψ)/(∂(Subscript[r, 1])))+(2/(Subscript[r, 2]))((∂ψ)/(∂(Subscript[r, 2])))+(4/(Subscript[r, 12]))((∂ψ)/(∂(Subscript[r, 12])))+((((Subscript[r, 1]^2) - (Subscript[r, 2]^2) + (Subscript[r, 12]^2)))/((Subscript[r, 1]) (Subscript[r, 12])))(((∂^2) ψ)/(∂(Subscript[r, 1]) ∂(Subscript[r, 12])))+((((Subscript[r, 2]^2) - (Subscript[r, 1]^2) + (Subscript[r, 12]^2)))/((Subscript[r, 2]) (Subscript[r, 12])))(((∂^2) ψ)/(∂(Subscript[r, 2]) ∂(Subscript[r, 12]))))

the potential energy operator, 𝒱, is

\(TraditionalForm
    \`\(𝒱 = 
        \(-\(𝒵\/r\_1\)\) - 𝒵\/r\_2 + 1\/r\_12;
      \)\)

and the Hamiltonian, ℋ, is

\(TraditionalForm
    \`ℋ(ψ_) := 
      𝒯(ψ) + 𝒱\ ψ\)

In Mathematica notation, the matrix element \(TraditionalForm
      \`⟨i, j, k⟩ ≡ 
        ⟨ϕ 
            |\(r\_1\%i\) \(r\_2\%j\) 
                r\_12\%k| ϕ⟩\) can be computed and saved using the following code:

\(TraditionalForm
    \`⟨i_, j_, k_⟩ := 
      \(⟨i, j, k⟩ = 
        \(1\/\((\(-2\))\)\^\(i + j + k + 3\)\) 
          ∂\_\({λ\_1, i + 1}, {λ\_2, j + 1}, {λ\_12, k + 1}\)\(\(2  π\^2\)\/\(\((λ\_1 + λ\_2)\) 
                  \((λ\_1 + λ\_12)\) 
                  \((λ\_2 + λ\_12)\)\)\)\)\)

Using this result we can compute matrix elements of any zeroth order approximation consisting of a linear combination of terms of the form  \(TraditionalForm
      \`\(P(r\_1, r\_2, r\_12)\) 
        e\^\(\(-λ\_1\) r\_1 - \(λ\_2\) r\_2 - \(λ\_12\) r\_12\)\), where P((Subscript[r, 1]), (Subscript[r, 2]), (Subscript[r, 12])) is a polynomial in the interparticle coordinates.

Consider the trial wavefunction \(TraditionalForm
      \`ψ\_0 = 
        e\^\(\(-𝒵\)\ r\_1 - 𝒵\ r\_2 + r\_12/2\)\).

\(TraditionalForm
    \`\(ψ\_0 = 
        E\^\(\(-𝒵\)\ r\_1 - 𝒵\ r\_2 + r\_12/2\)\ ;\)\)

\(TraditionalForm
    \`\(Λ = 
        {λ\_1 -> 𝒵, 
          λ\_2 -> 𝒵, 
          λ\_12 -> \(-\(1\/2\)\)};\)\)

\(TraditionalForm
    \`Simplify[
      ⟨0, 0, 0⟩ 
        /. Λ]\)

\(TraditionalForm
    \`\(p = 1 + 2\/\(\(a\)\(\ \)\) + 3\/b + 4  a + 5  a\ b;\)\)

Expand[(a^i) (b^j) p]

\(TraditionalForm
    \`% /. \(a\^\(i + m_. \)\) b\^\(j + n_. \) -> 
        ⟨m, n⟩\)

Here is the local energy:

\(TraditionalForm
    \`Expand[\(ℋ(ψ\_0)\)\/ψ\_0]\)

Using the suggested trick,

\(TraditionalForm
    \`ℰ(𝒵_) = 
      Factor[\(Expand[r\_1\%i\ r\_2\%j\ r\_12\%k\ %] /. r\_1\%\(i + l_. \)\ r\_2\%\(j + m_. \)\ r\_12\%\(k + n_. \) -> ⟨l, m, n⟩\)\/⟨0, 0, 0⟩ /. Λ]\)

and for helium we obtain

ℰ(2. )

The non-relativistic energy is

(Subscript[ℰ, NR]) = -2.903724377034;

so this answer is too high by (((Subscript[ℰ, NR]) - 
                ℰ(2))/(Subscript[ℰ, NR]))=1.66%.

ℰ(1. )

Since the ionization energy of hydrogen is ((-1)/2) and this energy is an upper bound for the \(TraditionalForm\`\(H\^-\)\) ground state energy, we cannot determine from this whether \(TraditionalForm\`\(H\^-\)\) is bound.

\(TraditionalForm
    \`\(Plot[ℰ(𝒵), {𝒵, 1, 5}, 
        Epilog -> 
          {Hue[1], PointSize[0.03], 
            Point /@ Table[
                {𝒵, ℰ(𝒵)}, 
                {𝒵, 1, 5}]}];\)\)

\(TraditionalForm
    \`\(ψ\_0 = 
        E\^\(\(-μ\)\ r\_1 - μ\ r\_2 - λ\ r\_12\);\)\)

\(TraditionalForm
    \`\(Λ = 
        {λ\_1 -> μ, λ\_2 -> μ, 
          λ\_12 -> λ};\)\)

\(TraditionalForm
    \`\(ℰ\_\(λ_, μ_\)\)(𝒵_) = 
      Simplify[\(Expand[\(r\_1\%i\ r\_2\%j\ r\_12\%k\ \(ℋ(ψ\_0)\)\)\/ψ\_0] /. r\_1\%\(i + l_. \)\ r\_2\%\(j + m_. \)\ r\_12\%\(k + n_. \) -> ⟨l, m, n⟩\)\/⟨0, 0, 0⟩ 
          /. Λ]\)

\(TraditionalForm
    \`FindMinimum[\(ℰ\_\(λ, μ\)\)(2), {μ, 2}, 
      {λ, \(-\(1\/2\)\)}]\)

This answer is, as expected, lower than that obtained in [b] because \(TraditionalForm
      \`ψ\_0 = 
        e\^\(\(-μ\)\ r\_1 - μ\ r\_2 - λ\ r\_12\)\) includes the previous example as a special case and we now have two variational parameters.  The minimum energy is quite close to the "exact" value: \(TraditionalForm
      \`\(ℰ\_NR - \(ℰ\_\(λ, μ\)\)(2)\)\/ℰ\_NR = 0.486  %\).

\(TraditionalForm
    \`FindMinimum[\(ℰ\_\(λ, μ\)\)(1), {μ, 1}, 
      {λ, \(-\(1\/2\)\)}]\)

Since the variational energy is less than the ionization energy of hydrogen we can conclude that \(TraditionalForm\`\(H\^-\)\) is (very weakly) bound.

. Pauli spin matrices and Dirac gamma matrices [20]

The Pauli spin matrices,

σ/:(Subscript[σ, 1])=(({{0, 1}, {1, 0}}));

σ/:(Subscript[σ, 2])=(({{0, -I}, {I, 0}}));

σ/:(Subscript[σ, 3])=(({{1, 0}, {0, -1}}));

along with the 2×2 identity matrix:

(Subscript[ℐ, n])_ := IdentityMatrix[n]

(Subscript[ℐ, 2])

are used in a matrix representation of spin operators in quantum mechanics.

\(TraditionalForm
    \`Table[σ\_i . σ\_i == ℐ\_2, {i, 3}]
      \)

(Subscript[ϵ, i])__ := Signature[{i}]

\(TraditionalForm
    \`test[{i_, j_, k_}] := 
      σ\_i . σ\_j - σ\_j . σ\_i == 
        2  I\ ϵ\_\(i, j, k\)\ σ\_k\)

test /@ Permutations[{1, 2, 3}]

Solve[(({{(Subscript[x, 1]) + I (Subscript[y, 1]), (Subscript[x, 2]) + I (Subscript[y, 2])}, {(Subscript[x, 3]) + I (Subscript[y, 3]), (Subscript[x, 4]) + I (Subscript[y, 4])}}))==(Subscript[c, 0]) (Subscript[ℐ, 2]) + 
                (∑_(i = 1)^3) (Subscript[c, i]) (Subscript[σ, i]),Table[(Subscript[c, i]), {i, 0, 3}]]//Simplify

((Subscript[a, 1]) (Subscript[σ, 1]) + (Subscript[a, 2]) (Subscript[σ, 2]) + 
              (Subscript[a, 3]) (Subscript[σ, 3])) . 
          ((Subscript[b, 1]) (Subscript[σ, 1]) + (Subscript[b, 2]) (Subscript[σ, 2]) + 
              (Subscript[b, 3]) (Subscript[σ, 3]))=={(Subscript[a, 1]), (Subscript[a, 2]), (Subscript[a, 3])} . {(Subscript[b, 1]), (Subscript[b, 2]), (Subscript[b, 3])} (Subscript[ℐ, 2])+i ({(Subscript[a, 1]), (Subscript[a, 2]), (Subscript[a, 3])}×{(Subscript[b, 1]), (Subscript[b, 2]), (Subscript[b, 3])}).{(Subscript[σ, 1]), (Subscript[σ, 2]), (Subscript[σ, 3])}. [3]

We can show this directly:

\(TraditionalForm
    \`\((\(a\_1\) σ\_1 + \(a\_2\) σ\_2 + 
              \(a\_3\) σ\_3)\) . 
          \((\(b\_1\) σ\_1 + \(b\_2\) σ\_2 + 
              \(b\_3\) σ\_3)\) == 
        {a\_1, a\_2, a\_3} . {b\_1, b\_2, b\_3} ℐ\_2 + 
          I\ \(({a\_1, a\_2, a\_3}\[Cross]{b\_1, b\_2, b\_3})\) . 
              {σ\_1, σ\_2, σ\_3} // ExpandAll\)

or indirectly:

𝒜 = {(Subscript[a, 1]), (Subscript[a, 2]), (Subscript[a, 3])}; 
    ℬ = {(Subscript[b, 1]), (Subscript[b, 2]), (Subscript[b, 3])}; 
    Σ = 
      {HoldForm[(Subscript[σ, 1])], HoldForm[(Subscript[σ, 2])], HoldForm[(Subscript[σ, 3])]}

\(TraditionalForm
    \`\(\((Σ . 𝒜)\) . 
          \((Σ . ℬ)\) // ReleaseHold\) // Expand
      \)

\(TraditionalForm
    \`𝒜 . ℬ\ ℐ\_2 + \ 
        \((I\ Σ . 
                \((𝒜\[Cross]ℬ)\) // 
            ReleaseHold)\) // Expand\)

% == %%

(Subscript[γ, k])=(({{𝒪, -i (Subscript[σ, k])}, {i (Subscript[σ, k]), 𝒪}})), k = 1,2,3,
(Subscript[γ, 4])=(({{(Subscript[ℐ, 2]), 𝒪}, {𝒪, -(Subscript[ℐ, 2])}})),
(Subscript[γ, 5]) =  
          (Subscript[γ, 1]) . (Subscript[γ, 2]) . (Subscript[γ, 3]) . (Subscript[γ, 4]),

<< LinearAlgebra`MatrixManipulation`

𝒪 = ZeroMatrix[2]

γ/:(Subscript[γ, k : (1 | 2 | 3)]):=(Subscript[γ, k])=BlockMatrix[(({{𝒪, -I (Subscript[σ, k])}, {I (Subscript[σ, k]), 𝒪}}))]

γ/:(Subscript[γ, 4])=BlockMatrix[(({{(Subscript[ℐ, 2]), 𝒪}, {𝒪, -(Subscript[ℐ, 2])}}))]

\(TraditionalForm
    \`γ /: γ\_5 = \ 
      γ\_1 . γ\_2 . γ\_3 . γ\_4\)

\(TraditionalForm
    \`δ\_\(μ_, ν_\) := If[μ == \ ν, 1, 0]\)

\(TraditionalForm
    \`Table[γ\_μ . γ\_ν + 
          γ\_ν . γ\_μ == 
        2  δ\_\(μ, ν\)\ ℐ\_4, {μ, 1, 4}, 
      {ν, 1, 4}]\)

\(TraditionalForm
    \`\[LeftDoubleBracketingBar]𝒜_, 
        ℬ_\[RightDoubleBracketingBar] := 
      𝒜 . ℬ - 
        ℬ . 𝒜\)

\(TraditionalForm
    \`σ /: σ\_\(μ_, ν_\) := 
      \(σ\_\(μ, ν\) = 
        \[LeftDoubleBracketingBar]γ\_μ, γ\_ν\[RightDoubleBracketingBar]\/\(2\ I\)\)\)

\(TraditionalForm
    \`Table[σ\_\(μ, ν\) == \ 
        \(-I\)\ γ\_μ . γ\_ν, {μ, 1, 4}, 
      {ν, 1, 4}]\)

\(TraditionalForm
    \`Γ = 
      {HoldForm[γ\_1], HoldForm[γ\_2], HoldForm[γ\_3], 
        HoldForm[γ\_4]}\)

𝒜 = {(Subscript[a, 1]), (Subscript[a, 2]), (Subscript[a, 3]), (Subscript[a, 4])};

ℬ = {(Subscript[b, 1]), (Subscript[b, 2]), (Subscript[b, 3]), (Subscript[b, 4])};

𝒞 = {(Subscript[c, 1]), (Subscript[c, 2]), (Subscript[c, 3]), (Subscript[c, 4])};

\(TraditionalForm
    \`\(∑\+\(μ = 1\)\%4 γ\_μ . 
              \((Γ . 𝒜)\) . 
              \((Γ . ℬ)\) . 
              \((Γ . 𝒞)\) . γ\_μ 
          == 
          \(-2\) \((Γ . 𝒞)\) . 
              \((Γ . ℬ)\) . 
              \((Γ . 𝒜)\) // ReleaseHold\) // 
      ExpandAll\)

\(TraditionalForm
    \`Tr[γ\_5 . \((Γ . 𝒜)\) . 
          \((Γ . ℬ)\) // ReleaseHold]\)

. Series [20]

.. Discrete nonlinear transform [12]

Consider the formal Taylor series

\(TraditionalForm
    \`\(\(∑\+\(n = 0\)\%∞\(\((\(-1\))\)\^n\) \(a\_n\) x\^n\)\(,
    \)\)\)

where (Subscript[a, 0]) = 1and the continued fraction

\(TraditionalForm
    \`\(\(1\/\(1 + 
            \(\(b\_1\) x\)\/\(1 + \(\(b\_2\) x\)\/\(1 + \(\(b\_3\) x\)\/\(1 + …\)\)\)\)\)\(·\)\)\)

m = 5;

\(TraditionalForm
    \`as = Solve[
        ∑\+\(n = 0\)\%m\(\((\(-1\))\)\^n\) \(a\_n\) x\^n - 
            Fold[1\/\(1 + #1\ #2\) &, 1, Table[x\ b\_n, {n, m}] // Reverse] + 
            O[x]\^\(m + 1\) == 0, Table[a\_i, {i, 0, m}]]\)

\(TraditionalForm
    \`bs = Solve[
          ∑\+\(n = 0\)\%m\(\((\(-1\))\)\^n\) \(a\_n\) x\^n - 
              Fold[1\/\(1 + #1\ #2\) &, 1, Table[x\ b\_n, {n, m}] // Reverse] 
              + O[x]\^\(m + 1\) == 0, Table[b\_i, {i, m}]] // Simplify
        \)

Self-consistency is easy to verify:

First[bs] /. First[as] // Simplify

Table[a/:(Subscript[a, n])=(((-1))^n) (Subscript[E, 2 n]),{n, 5}]

First[bs] // Reverse

\(TraditionalForm
    \`Fold[1\/\(1 + #1\ #2\) &, 1, Table[x\ n\^2, {n, m}] // Reverse]\)

sec(x) + O[x]^16

Table[(((I^n) (Subscript[E, n]))/(n!)),{n, 0, 14}]

Hence sec(x)=(∑_(n = 1)^∞)(i^n)(Subscript[E, n])((x^n)/(n!)).

.. Kepler Equation [8]

The Kepler equation, s = u + e sin s, arises in celestial mechanics. Defining s = u + a reduces the equation to a = e sin(u + a) which can be formally solved by introducing the recurrence relation (Subscript[a, k]) = e sin(u + (Subscript[a, k - 1])). Since e is a small parameter, as k -> ∞, (Subscript[a, k - 1]) -> (Subscript[a, k]) -> a, and hence the Kepler equation is formally satisfied.

Clear[a]

a /: (Subscript[a, 0]) = 0;

\(TraditionalForm
    \`a /: a\_k_ := \(a\_k = e\ Sin[u + a\_\(k - 1\)] + O[e]\^\(k + 1\)\)\)

TrigReduce /@ (Subscript[a, 5])

Collect[%, sin(_)]

(Subscript[a, 5]) - e Sin[u + (Subscript[a, 5])] // Simplify