© Kamal K. Barley, Andreas Ruffing, and Sergei K. Suslov
“Old Quantum Mechanics by Bohr and Sommerfeld from Modern Perspective”

Dedicated to the 100th anniversary of the birth of wave mechanics.

© К. Барли, А. Руффинг, С.К. Суслов
“Старая квантовая механика Бора и Зоммерфельда с современной точки зрения”

Успехи физических наук, Январь 2026 г., Том 196, № 1, с. 83-103

DOI: https://doi.org/10.3367/UFNr.2025.08.040014

Note: You will be able to run each section independently from the others.

(Last modified/executed on December 14, 2026;  12:25 AM Arizona time.)

Evaluation of the integral (3.15)

f[ψ_]:=(((ε^2) ((Sin[ψ])^2))/(((1+ε Cos[ψ]))^2))

f[ψ]

(* Integral (3.15): *)Integrate[f[ψ],{ψ,0,2*Pi},Assumptions->1>ε>0]

(* Antiderivative *)Integrate[f[ψ],ψ,Assumptions->1>ε>0]

Derivation of Equations (3.18)-(3.21): - See also Appendix D

This is about solution of (3.4). Start with (3.9): s=1/r

s[θ_]:=(1+ϵ*Cos[ω*θ])/(a*(1-(ϵ^2)))

s[θ]

(* The first derivative needed in (3.4): *) D[s[θ],θ]

(%)^2

(* in terms of cosine *)%/. ((Sin[θ ω])^2)->1-((Cos[θ ω])^2)

((s[θ])^2)

(* Verify *) Simplify[%/.Cos[θ ω]->((-1+a S-a S (ϵ^2))/ϵ)]

(* Indeed,from (3.9): *) Solve[((1+ϵ C)/(a (1-(ϵ^2))))-S==0,C]

Simplify[/.((Cos[θ ω])^2)->((((-1+a S (1-(ϵ^2))))^2)/(ϵ^2))]

Expand[%,S]

Simplify[%+(((S))^2)]

Expand[%,S]

As a result, this expression for (ds/dθ)^2 + s^2 is correct.

Working on the quadratic polynomial in (3.4) - l.h.s minus r.h.s:

((((En/(m (c^2)))+((Z (e^2) S )/(m (c^2)))))^2)-1-(((((ℏ (Subscript[n, θ]))/(m c))))^2)*(((ω^2)/((a^2) (-1+(ϵ^2))))-((2 S (ω^2))/(a (-1+(ϵ^2))))-(S^2) (-1+(ω^2)))

%/.(e^2)->α*ℏ*c

Expand[%,S]

(* Denote this quadratic polynomials as P: *)P=-1+((En^2)/((c^4) (m^2)))+((2 En S Z α ℏ)/((c^3) (m^2)))+(((S^2) (Z^2) (α^2) (ℏ^2))/((c^2) (m^2)))-(((ω^2) (ℏ^2) (Subscript[n, θ]^2))/((a^2) (c^2) (m^2) (-1+(ϵ^2))))+((2 S (ω^2) (ℏ^2) (Subscript[n, θ]^2))/(a (c^2) (m^2) (-1+(ϵ^2))))+(((S^2) (-1+(ω^2)) (ℏ^2) (Subscript[n, θ]^2))/((c^2) (m^2)))

Collect[P,S]

The last coefficient of  S^2 must be zero:

Factor[((((Z^2) (α^2) (ℏ^2))/((c^2) (m^2)))+(((-1+(ω^2)) (ℏ^2) (Subscript[n, θ]^2))/((c^2) (m^2))))]

(* Derivation of (3.19): *) Solve[(Z^2) (α^2)-(Subscript[n, θ]^2)+(ω^2) (Subscript[n, θ]^2)==0,ω]

The last but one coefficient of S must be zero. By (3.17):

Solve[((En Z α)/((Subscript[n, θ]^2)*c*ℏ*(ω^2)))==(1/(a (1-(ϵ^2)))),En]

(* Verification: *) ((2 En Z α ℏ)/((c^3) (m^2)))+((2 (ω^2) (ℏ^2) (Subscript[n, θ]^2))/(a (c^2) (m^2) (-1+(ϵ^2))))/.En->-((c (ω^2) ℏ (Subscript[n, θ]^2))/(a Z α (-1+(ϵ^2))))

The constant term must be zero:

-1+((En^2)/((c^4) (m^2)))-(((ω^2) (ℏ^2) (Subscript[n, θ]^2))/((a^2) (c^2) (m^2) (-1+(ϵ^2))))/.(1/((-1+(ϵ^2))))->-((((Subscript[n, r])+(Sqrt[-(Z^2) (α^2)+(Subscript[n, θ]^2)]))^2)/((ω^2) (Subscript[n, θ]^2)))

In view of (3.16)-(3.17):

-1+((En^2)/((c^4) (m^2)))+(((ℏ^2) ((((Subscript[n, r])+(Sqrt[-(Z^2) (α^2)+(Subscript[n, θ]^2)])))^2))/((a^2) (c^2) (m^2)))/.En->((c ℏ ((((Subscript[n, r])+(Sqrt[-(Z^2) (α^2)+(Subscript[n, θ]^2)])))^2))/(a Z α))

Replace (a^2) by A by hand:

-1+(((ℏ^2) ((((Subscript[n, r])+(Sqrt[-(Z^2) (α^2)+(Subscript[n, θ]^2)])))^2))/(A (c^2) (m^2)))+(((ℏ^2) ((((Subscript[n, r])+(Sqrt[-(Z^2) (α^2)+(Subscript[n, θ]^2)])))^4))/(A (c^2) (m^2) (Z^2) (α^2)))==0

Solve[(((ℏ^2) ((((Subscript[n, r])+(Sqrt[-(Z^2) (α^2)+(Subscript[n, θ]^2)])))^2))/((c^2) (m^2)))+(((ℏ^2) ((((Subscript[n, r])+(Sqrt[-(Z^2) (α^2)+(Subscript[n, θ]^2)])))^4))/((c^2) (m^2) (Z^2) (α^2)))==A,A]

Factor[(((ℏ^2) ((((Subscript[n, r])+(Sqrt[-(Z^2) (α^2)+(Subscript[n, θ]^2)])))^2))/((c^2) (m^2)))+(((ℏ^2) ((((Subscript[n, r])+(Sqrt[-(Z^2) (α^2)+(Subscript[n, θ]^2)])))^4))/((c^2) (m^2) (Z^2) (α^2)))]

Expand[((((Subscript[n, r])+(Sqrt[-(Z^2) (α^2)+(Subscript[n, θ]^2)])))^2)+(Z^2) (α^2)]

Therefore we can write:

(((ℏ^2) ((((Subscript[n, r])+(Sqrt[-(Z^2) (α^2)+(Subscript[n, θ]^2)])))^2)*(((((Subscript[n, r])+(Sqrt[-(Z^2) (α^2)+(Subscript[n, θ]^2)])))^2)+(Z^2) (α^2)))/((c^2) (m^2) (Z^2) (α^2)))

(* Verify: *) Simplify[%-(((ℏ^2) ((((Subscript[n, r])+(Sqrt[-(Z^2) (α^2)+(Subscript[n, θ]^2)])))^2) ((Subscript[n, r]^2)+(Subscript[n, θ]^2)+2 (Subscript[n, r]) (Sqrt[-(Z^2) (α^2)+(Subscript[n, θ]^2)])))/((c^2) (m^2) (Z^2) (α^2)))]

For a=\sqrt[A]:

Sqrt[(((ℏ^2) ((((Subscript[n, r])+(Sqrt[-(Z^2) (α^2)+(Subscript[n, θ]^2)])))^2) ((Z^2) (α^2)+((((Subscript[n, r])+(Sqrt[-(Z^2) (α^2)+(Subscript[n, θ]^2)])))^2)))/((c^2) (m^2) (Z^2) (α^2)))]

As a result we obtain (3.21):

((ℏ )/(c m Z α))*((Subscript[n, r])+(Sqrt[-(Z^2) (α^2)+(Subscript[n, θ]^2)])) (Sqrt[(Z^2) (α^2)+((((Subscript[n, r])+(Sqrt[-(Z^2) (α^2)+(Subscript[n, θ]^2)])))^2)])

By (3.16)-(3.17) for the energy levels:

((c ℏ ((((Subscript[n, r])+(Sqrt[-(Z^2) (α^2)+(Subscript[n, θ]^2)])))^2))/(a Z α))/.a->((ℏ ((Subscript[n, r])+(Sqrt[-(Z^2) (α^2)+(Subscript[n, θ]^2)])) (Sqrt[(Z^2) (α^2)+((((Subscript[n, r])+(Sqrt[-(Z^2) (α^2)+(Subscript[n, θ]^2)])))^2)]))/(c m Z α))

This the Sommerfeld formula (3.18):

Simplify[((Z^2) (α^2)+((((Subscript[n, r])+(Sqrt[-(Z^2) (α^2)+(Subscript[n, θ]^2)])))^2))/((((Subscript[n, r])+(Sqrt[-(Z^2) (α^2)+(Subscript[n, θ]^2)])))^2)-(1+(((Z^2) (α^2))/((((Subscript[n, r])+(Sqrt[-(Z^2) (α^2)+(Subscript[n, θ]^2)])))^2)))]

(* Equation (3.18): *) m*c^2/((1+((Z^2) (α^2))/((Subscript[n, r])+(Sqrt[-(Z^2) (α^2)+(Subscript[n, θ]^2)]))^2))^(1/2)

Asymptotics (3.22)-(3.25):

(* Definition (3.19): *) ω[α_]:=(Sqrt[1-((Z^2) (α^2))/(Subscript[n, θ]^2)])

ω[α]

(* Equation (3.22): *)FullSimplify[PowerExpand[Series[ω[α],{α,0,5},Assumptions->{-(Z^2) (α^2)+(Subscript[n, θ]^2)>0,(Subscript[n, θ])>0}]]]

(* Definition (3.20): *) ϵ[α_]:=((Sqrt[(Subscript[n, r]) ((Subscript[n, r])+2 (Sqrt[-(Z^2) (α^2)+(Subscript[n, θ]^2)]))])/((Subscript[n, r])+(Sqrt[-(Z^2) (α^2)+(Subscript[n, θ]^2)])))

ϵ[α]

(* Equation (3.23): *) FullSimplify[PowerExpand[Series[ϵ[α],{α,0,5},Assumptions->{-(Z^2) (α^2)+(Subscript[n, θ]^2)>0,(Subscript[n, θ])>0,(Subscript[n, r])>0}]]]

(* Definition (3.21) in Bohr's atomic units: *) a[α_]:=((Subscript[n, r])+(Sqrt[-(Z^2) (α^2)+(Subscript[n, θ]^2)])) (Sqrt[(Z^2) (α^2)+((((Subscript[n, r])+(Sqrt[-(Z^2) (α^2)+(Subscript[n, θ]^2)])))^2)])/Z

a[α]

(* Formula (3.24)in Bohr's atomic units: *)FullSimplify[PowerExpand[Series[Z*a[α],{α,0,5},Assumptions->{-(Z^2) (α^2)+(Subscript[n, θ]^2)>0,(Subscript[n, θ])>0}]]]

(* Definition (3.18): *) En[α_]:=m*c^2/((1+((Z^2) (α^2))/((Subscript[n, r])+(Sqrt[-(Z^2) (α^2)+(Subscript[n, θ]^2)]))^2))^(1/2)

(* Formula (3.25): *)FullSimplify[PowerExpand[Series[En[α]/((c^2) m),{α,0,5},Assumptions->{-(Z^2) (α^2)+(Subscript[n, θ]^2)>0,(Subscript[n, θ])>0}]]]

Animation of the relativistic Kepler motion

(*---Constants and Setup---*)(*Pretty print subscript-style variables for display only*)MakeBoxes[nt,StandardForm]:=SubscriptBox["n","t"]
MakeBoxes[nr,StandardForm]:=SubscriptBox["n","r"]

(*Physical constants*)
Z=79; (*Atomic number of gold, Au*)
α=1/137.036; (*Fine structure constant*)

(*Quantum numbers*)
ntVal=1;
nrVal=1;

(*---Define symbolic formulas as functions---*)

(*Eq.(3.19)*)
ωFunc[nt_]:=Sqrt[nt^2-Z^2 α^2]/nt

(*Eq.(3.20)*)
ϵFunc[nt_,nr_]:=Sqrt[nr] Sqrt[(nr+2 Sqrt[nt^2-Z^2 α^2])]/(nr+Sqrt[nt^2-Z^2 α^2])

(*Eq.(3.21)*)
aFunc[nt_,nr_]:=((nr+Sqrt[nt^2-Z^2 α^2]) Sqrt[Z^2 α^2+(nr+Sqrt[nt^2-Z^2 α^2])^2])/Z
(*---Compute numeric values---*)

ωVal=ωFunc[ntVal]//N
ϵVal=ϵFunc[ntVal,nrVal]//N
aVal=aFunc[ntVal,nrVal]//N

(*---Define relativistic radial orbit r(t) from Eq.(3.9)---*)
r[t_]:=(aVal (1-ϵVal^2))/(1+ϵVal Cos[ωVal t])

(*---Polar plot the relativistic elliptical orbit---*)
PolarPlot[r[t],{t,0,2Pi},PlotLabel->"Relativistic Elliptical Orbit",AxesLabel->{"x","y"},PlotRange->All,PlotStyle->Thick]
(*---Polar plot the relativistic elliptical orbit---*)
PolarPlot[r[t],{t,0,4 Pi},PlotLabel->"Relativistic Elliptical Orbit",AxesLabel->{"x","y"},PlotRange->All,PlotStyle->Thick]

(*Animation in Cartesian Coordinates*)Manipulate[ParametricPlot[{r[t] Cos[t],r[t] Sin[t]},{t,0,tMax},PlotStyle->{Blue,Thick},AxesLabel->{"x","y"},PlotLabel->Style["Relativistic Elliptical Orbit (Animated, Cartesian)",14,Bold],PlotRange->{{-0.084,0.084},{-0.084,0.084}},(*Fixed range ensures full visibility from (Subscript[r, min])=0.004683 to (Subscript[r, max])=0.083014*)AspectRatio->1,PlotPoints->1000,PerformanceGoal->"Quality"],{tMax,0.1,23 Pi,Appearance->"Labeled"}]

Summary: 
In the above animation of the relativistic Kepler motion of an electron in hydrogen-like gold ion, Au, when Z=79, the nucleus is situated  in the fixed focus at the origin.  For this quantum state, we choose (Subscript[n, r])=nr=1 and (Subscript[n, θ])=nt=1 in the fine structure formula (3.18). By (3.19)-(3.21), one gets: ω=0.817104,  ϵ=0.893193, and a=0.0438488, in Bohr’s atomic units. 
The perihelion and aphelion move along two concentric circles around the nucleus with radii:  (Subscript[r, min])=a (1-ϵ)=0.004683 and (Subscript[r, max])=a (1+ϵ)=0.083014, respectively.  In the animation, the geometrical loci of the successive perihelia and aphelia, the outer and inner envelopes of the orbit, see Figure 3, are not shown for simplicity.
Rotation of the ellipse over one period is given by Δθ=2π(1/ω-1)=1.4064.

Fine Structure Formulas in Wave Mechanics

Relativistic Schrödinger equation

Here, from (3.42)-(3.43), we derive the expansion (3.44) for the relativistic Schrödinger equation (3.39).

(* Definition (3.44) *) EnergyRelSchr[μ_]:=(1/(Sqrt[1+((μ^2)/(((n+(1/2)+(Sqrt[((((1/2)+l))^2)-(μ^2)])))^2))]))

EnergyRelSchr[μ]

(* Asymptotic (3.44) *)Simplify[PowerExpand[Series[EnergyRelSchr[μ],{μ,0,5},Assumptions->2*l+1>0]]]

FullSimplify[%/.n->n-l-1]

This is the nonrelativistic expansion (3.44).

Section 6:  A “mistake” that Schrödinger never made (Figure 5):

Simplify[EnergyRelSchr[μ]/.n->n-l/2]

FullSimplify[PowerExpand[Series[%,{μ,0,5},Assumptions->l>0]]]

Thus, the non-relativistic term is wrongfully given by:

-((μ^2)/(2*(((1+l/2+n))^2)))

Dirac Equation

Here, from (3.47)-(3.48), we derive the expansion (3.49) for the Dirac equation (3.45).

(1/(Sqrt[1+(μ^2)/(((n+ν))^2)]))

(* Definition (3.47)-(3.48): *) EnergyDirac[μ_]=%/.ν->Sqrt[(j+1/2)^2-(μ^2)];

EnergyDirac[μ]

Simplify[PowerExpand[Series[EnergyDirac[μ],{μ,0,5},Assumptions->2*j+1>0]]]

FullSimplify[%/.n->n-j-1/2]

This is the expansion (3.49), indeed:

1-((μ^2)/(2 (n^2)))+((((3/4)-(n/((1/2)+j))) (μ^4))/(2 (n^4)))+((O[μ])^6)

(* Verification:*) Simplify[%%-%]

Next term in (3.51) for the Dirac equation.

EnergyDiracExp[μ_]=1-((2 (μ^2))/(((1+2 j+2 n))^2))-((2 (1+2 j+8 n) (μ^4))/((1+2 j) (((1+2 j+2 n))^4)))

EnergyDirac[μ]-EnergyDiracExp[μ]

Simplify[PowerExpand[Series[EnergyDirac[μ]-EnergyDiracExp[μ],{μ,0,7},Assumptions->2*j+1>0]]]

FullSimplify[%/.n->n-j-1/2]

The last term can be rewritten as follows:

-(( (μ^6))/(16  (n^3)))*(((5 )/( (n^3)))-((24 )/((1+2 j) (n^2)))+((24 )/((((1+2 j))^2) n))+(16/(((1+2 j))^3)))

(* Verification: *) Simplify[%%-%]

The next term can also be rewritten as in (3.51):

-(( (μ^6))/(4  (n^6)))*(((5 )/( 4))-((3*n )/(j+1/2))+((3*(n^2))/(2*(((j+1/2))^2)))+((n^3)/(2*(((j+1/2))^3))))

(* Verification: *) Simplify[%%%-%]

Next terms as in (3.51)-(3.52) both together:

Simplify[PowerExpand[Series[EnergyDirac[μ],{μ,0,9},Assumptions->2*j+1>0]]]

FullSimplify[%/.n->n-j-1/2]

EnergyDiracExpMore[μ_]:=1-((μ^2)/(2 (n^2)))+((((3/4)-(n/((1/2)+j))) (μ^4))/(2 (n^4)))-((((5/4)-((3 n)/((1/2)+j))+((3 (n^2))/(2 ((((1/2)+j))^2)))+((n^3)/(2 ((((1/2)+j))^3)))) (μ^6))/(4 (n^6)))-(((-(35/8)+((15 n)/((1/2)+j))-((15 (n^2))/((((1/2)+j))^2))+((n^3)/((((1/2)+j))^3))+((3 (n^4))/((((1/2)+j))^4))+((n^5)/((((1/2)+j))^5))) (μ^8))/(16 (n^8)))

EnergyDiracExpMore[μ]

EnergyDiracExpMoreJ[μ_]:=1-((μ^2)/(2 (n^2)))+(((3-((8 n)/(1+2 j))) (μ^4))/(8 (n^4)))-(((5 (((1+2 j))^3)-24 (((1+2 j))^2) n+24 (1+2 j) (n^2)+16 (n^3)) (μ^6))/(16 ((((1+2 j))^3) (n^6))))-(1/(128 ((((1+2 j))^5) (n^8))))(-35 (((1+2 j))^5)+240 (((1+2 j))^4) n-480 (((1+2 j))^3) (n^2)+64 (((1+2 j))^2) (n^3)+384 (1+2 j) (n^4)+256 (n^5)) (μ^8)

EnergyDiracExpMoreJ[μ]

(* Verification: *) Simplify[EnergyDiracExpMore[μ]-EnergyDiracExpMoreJ[μ]]

Appendix F

Introduction

In this Introductory Section, the general formulas are provided. They are stored in global variables that will be used in all the subsequent sections. For this purpose, allow Mathematica to evaluate all initialization cells. After that you will be able to run each case independently from the others.
Use the following command, if you need to save/export graphs created by Mathematica into the current notebook directory:

SetDirectory[NotebookDirectory[]];

(*---Constants and Setup---*)(*Pretty print subscript-style variables for display only*)MakeBoxes[nt,StandardForm]:=SubscriptBox["n","t"]
MakeBoxes[nr,StandardForm]:=SubscriptBox["n","r"]
(*Physical constants*)
α=1/137.036; (*Fine structure constant*)

(*---Define symbolic formulas as functions---*)

(*Eq.(3.19)/(1.3)*)
ωFunc[nt_]:=Sqrt[nt^2-Z^2 α^2]/nt

(*Eq.(3.20)/(1.4)*)
ϵFunc[nt_,nr_]:=Sqrt[nr] Sqrt[(nr+2 Sqrt[nt^2-Z^2 α^2])]/(nr+Sqrt[nt^2-Z^2 α^2])

(*Eq.(3.21)/(1.5)*)
aFunc[nt_,nr_]:=((nr+Sqrt[nt^2-Z^2 α^2]) Sqrt[Z^2 α^2+(nr+Sqrt[nt^2-Z^2 α^2])^2])/Z

(*Eq.(3.18)/(1.1) in mc^2 units*)
eFunc[nt_,nr_]:=(1+(Z^2 α^2)/(nr+Sqrt[nt^2-Z^2 α^2])^2)^(-1/2)

Uranium, U, Z=92

Clear[ntVal,nrVal,ωVal,aVal,rMin,rMax,Z,DeltaTheta, PeriodTheta, eVal]

(*Physical constants*)
Z=92; (*Atomic number of Uranium, U*)

(*Quantum numbers*)
ntVal=1;
nrVal=1;

(*---Compute numeric values---*)
ωVal=ωFunc[ntVal]//N
ϵVal=ϵFunc[ntVal,nrVal]//N
aVal=aFunc[ntVal,nrVal]//N
{rMin=aVal*(1-ϵVal) //N, rMax=aVal*(1+ϵVal) //N}
DeltaTheta=2*Pi*((1/ωVal) -1) //N
PeriodTheta=2*Pi*(1/ωVal) //N
eVal=eFunc[ntVal,nrVal]//N
{Z,Z/137 //N}
wNum=2*((1/ωVal) -1) //N

(*---Define relativistic radial orbit r(t) from Eq.(3.9)/(1.2)---*)
r[t_]:=(aVal (1-ϵVal^2))/(1+ϵVal Cos[ωVal t])

(*---Parametric Plots in Polar Coordinates---From the first (Subscript[r, min]) to conseqitive (Subscript[r, min]) ones: n=1,1.5,1,2,3,4---*)
PolarPlot[r[t],{t,0,PeriodTheta},PlotLabel->"Relativistic Elliptical Orbit: one period, (Subscript[r, min]) to (Subscript[r, min])",AxesLabel->{"x","y"},PlotRange->All,PlotStyle->{Blue, Thick}]
 PolarPlot[{rMin,r[t],rMax},{t,0,PeriodTheta},PlotLabel->"Relativistic Elliptical Orbit: one period, (Subscript[r, min]) to (Subscript[r, min]) with (Subscript[r, min]) and (Subscript[r, max])",AxesLabel->{"x","y"},PlotRange->All,PlotStyle->{Blue, Thick}]
PolarPlot[r[t],{t,0,1.5*PeriodTheta},PlotLabel->"Relativistic Elliptical Orbit: one and a half periods, (Subscript[r, min]) to (Subscript[r, max])",AxesLabel->{"x","y"},PlotRange->All,PlotStyle->{Blue, Thick}]
PolarPlot[r[t],{t,0.5*PeriodTheta,1.5*PeriodTheta},PlotLabel->"Relativistic Elliptical Orbit: one periond, (Subscript[r, max]) to (Subscript[r, max])",AxesLabel->{"x","y"},PlotRange->All,PlotStyle->{Blue, Thick}]
PolarPlot[{rMin,r[t],rMax},{t,0.5*PeriodTheta,1.5*PeriodTheta},PlotLabel->"Relativistic Elliptical Orbit: one periond, (Subscript[r, max]) to (Subscript[r, max]) with (Subscript[r, min]) and (Subscript[r, max])",AxesLabel->{"x","y"},PlotRange->All,PlotStyle->{Blue, Thick}]
PolarPlot[r[t],{t,0,2*PeriodTheta},PlotLabel->"Relativistic Elliptical Orbit: 2 rotations",AxesLabel->{"x","y"},PlotRange->All,PlotStyle->{Blue,Thick}]
PolarPlot[r[t],{t,PeriodTheta,2*PeriodTheta},PlotLabel->"Relativistic Elliptical Orbit: the second period",AxesLabel->{"x","y"},PlotRange->All,PlotStyle->{Blue,Thick}]
PolarPlot[r[t],{t,0,3*PeriodTheta},PlotLabel->"Relativistic Elliptical Orbit: 3 rotations",AxesLabel->{"x","y"},PlotRange->All,PlotStyle->{Blue,Thick}]
PolarPlot[r[t],{t,0,4*PeriodTheta},PlotLabel->"Relativistic Elliptical Orbit; 4 rotations",AxesLabel->{"x","y"},PlotRange->All,PlotStyle->{Blue,Thick}]

Export["UraniumIon91.pdf",%%]

r[t]

(*Animation in Cartesian Coordinates*)Manipulate[ParametricPlot[{r[t] Cos[t],r[t] Sin[t]},{t,0,tMax},PlotStyle->{Blue,Thick},AxesLabel->{"x","y"},PlotLabel->Style["Relativistic Elliptical Orbit(Animated, Cartesian)",14,Bold],PlotRange->{{-rMax,rMax},{-rMax,rMax}},(*Fixed range ensures full visibility from (Subscript[r, min]) to (Subscript[r, max])*)AspectRatio->1,PlotPoints->1000,PerformanceGoal->"Quality"],{tMax,0.01,2*4*PeriodTheta,Appearance->"Labeled"}]

(*Animation in Rotating Cartesian Coordinates -- only radial motion*)Manipulate[ParametricPlot[{r[t] ,0},{t,0,tMax},PlotStyle->{Blue,Thick},AxesLabel->{"x","y"},PlotLabel->Style["Relativistic Elliptical Orbit(Animated, Cartesian)",14,Bold],PlotRange->{{-rMax,rMax},{-rMin,rMin}},(*Fixed range ensures full visibility from (Subscript[r, min]) to (Subscript[r, max])*)AspectRatio->1,PlotPoints->1000,PerformanceGoal->"Quality"],{tMax,0.01,PeriodTheta/2,Appearance->"Labeled"}]

Summary: 
In the above animation of the relativistic Kepler motion of an electron in hydrogen-like Uranium, U, when Z=92, the nucleus is situated  in the fixed focus at the origin.  For this quantum state, we choose (Subscript[n, r])=nr=1 and (Subscript[n, θ])=nt=1 in the fine structure formula (3.18). By (3.19)-(3.21), one gets: ω=0.741135,  ϵ=0.904882, and a=0.0353163, in Bohr’s atomic units. 
The perihelion and aphelion move along two concentric circles around the nucleus with radii:  (Subscript[r, min])=a (1-ϵ)=0.003359209 and (Subscript[r, max])=a (1+ϵ)=0.067273472, respectively.  In the animation, the geometrical loci of the successive perihelia and aphelia, the outer and inner envelopes of the orbit, see Figure 3, are not shown for simplicity.
CCW Rotation of the ellipse over one period is given by Δθ=2π(1/ω-1)=2.19461.  ‘Winding’ number is one.

Verification:

{ωVal, ϵVal, aVal}

{rMin, rMax}

{DeltaTheta,PeriodTheta}

{eVal,Z,ntVal,nrVal}

r[t]

(*End of Uranium, Z=92, section*)

Oganesson, Og, Z=118

Clear[ntVal,nrVal,ωVal,aVal,rMin,rMax,Z,DeltaTheta, PeriodTheta, eVal,wNum]

(*Physical constants*)
Z=118; (*Atomic number of Oganesson,Og*)

(*Quantum numbers*)
ntVal=1;
nrVal=1;

(*---Compute numeric values---*)
ωVal=ωFunc[ntVal]//N
ϵVal=ϵFunc[ntVal,nrVal]//N
aVal=aFunc[ntVal,nrVal]//N
{rMin=aVal*(1-ϵVal) //N, rMax=aVal*(1+ϵVal) //N}
DeltaTheta=2*Pi*((1/ωVal) -1) //N
PeriodTheta=2*Pi*(1/ωVal) //N
eVal=eFunc[ntVal,nrVal]//N
{Z,Z/137 //N}
wNum=2*((1/ωVal) -1) //N

(*---Define relativistic radial orbit r(t) from Eq.(3.9)/(1.2)---*)
r[t_]:=(aVal (1-ϵVal^2))/(1+ϵVal Cos[ωVal t])

(*---Parametric Plots in Polar Coordinates---From the first (Subscript[r, min]) to conseqitive (Subscript[r, min]) ones: n=1,1.5,2,3,4---*)
 PolarPlot[r[t],{t,0,PeriodTheta},PlotLabel->"Relativistic Elliptical Orbit: 1 rotation per period",AxesLabel->{"x","y"},PlotRange->All,PlotStyle->{Blue, Thick}]
PolarPlot[r[t],{t,PeriodTheta,2*PeriodTheta},PlotLabel->"Relativistic Elliptical Orbit: the second period",AxesLabel->{"x","y"},PlotRange->All,PlotStyle->{Blue,Thick}]
PolarPlot[r[t],{t,0.5*PeriodTheta,1.5*PeriodTheta},PlotLabel->"Relativistic Elliptical Orbit: 1 rotation, max->max",AxesLabel->{"x","y"},PlotRange->All,PlotStyle->{Blue, Thick}]
PolarPlot[r[t],{t,0,1.5*PeriodTheta},PlotLabel->"Relativistic Elliptical Orbit: 1.5 rotations",AxesLabel->{"x","y"},PlotRange->All,PlotStyle->{Blue, Thick}]
PolarPlot[r[t],{t,0,2*PeriodTheta},PlotLabel->"Relativistic Elliptical Orbit: 2 rotations",AxesLabel->{"x","y"},PlotRange->All,PlotStyle->{Blue,Thick}]
PolarPlot[r[t],{t,0,3*PeriodTheta},PlotLabel->"Relativistic Elliptical Orbit: 3 rotation",AxesLabel->{"x","y"},PlotRange->All,PlotStyle->{Blue,Thick}]
PolarPlot[r[t],{t,0,4*PeriodTheta},PlotLabel->"Relativistic Elliptical Orbit; 4 rotation",AxesLabel->{"x","y"},PlotRange->All,PlotStyle->{Blue,Thick}]

Export["OganessonIon117.pdf",%%%]

r[t]

(*Animation in Cartesian Coordinates*)Manipulate[ParametricPlot[{r[t] Cos[t],r[t] Sin[t]},{t,0,tMax},PlotStyle->{Blue,Thick},AxesLabel->{"x","y"},PlotLabel->Style["Relativistic Elliptical Orbit(Animated, Cartesian)",14,Bold],PlotRange->{{-rMax,rMax},{-rMax,rMax}},(*Fixed range ensures full visibility from (Subscript[r, min]) to (Subscript[r, max])*)AspectRatio->1,PlotPoints->1000,PerformanceGoal->"Quality"],{tMax,0.01,2*4*PeriodTheta,Appearance->"Labeled"}]

Summary: 
In the above animation of the relativistic Kepler motion of an electron in hydrogen-like Oganesson, Og, when Z=118, the nucleus is situated  in the fixed focus at the origin.  For this quantum state, we choose (Subscript[n, r])=nr=1 and (Subscript[n, θ])=nt=1 in the fine structure formula (3.18). By (3.19)-(3.21), one gets: ω= 0.508457,  ϵ=0.941479, and a=0.0222041, in Bohr’s atomic units. 
The perihelion and aphelion move along two concentric circles around the nucleus with radii:   (Subscript[r, min])=a (1-ϵ)=0.001299396360213745` and (Subscript[r, max])=a (1+ϵ)=0.0431087, respectively.  In the animation, the geometrical loci of the successive perihelia and aphelia, the outer and inner envelopes of the orbit, see Figure 3, are not shown for simplicity. ‘Winding’ number:  1.93347, about two. CW Rotation of the ellipse over one period is given by Δθ=2π(1/ω-1)=6.07418.  
At the point of self-intersection one gets, approximately,  r[2.3799243523877007`]=r[8.52445631673413`]=0.00281928 with  r[2.3799243523877007`]-r[8.52445631673413`]=-1.734723475976807`*^-18.

Verification:

{ωVal, ϵVal, aVal}

{rMin, rMax}

{DeltaTheta,PeriodTheta,wNum}

{eVal,Z,ntVal,nrVal}

(* points of self-intersection *)r[9.283329709624294`]-r[3.074037420536865`]

{r[9.283329709624294`],r[3.074037420536865`]}

ωVal*(9.283329709624294`+3.074037420536865`)/2//N

3.141592653589793`-Pi//N

r0[t_]:=0.004308954533019283`

r[t_]:=((0.0025227511535937256`)/(1+0.9414793136551063` Cos[0.5084566348962754` t]))

{r0[t], r[t]}

PeriodTheta

PolarPlot[{r0[t], r[t]},{t,0,1.05*PeriodTheta},PlotLabel->"Relativistic 'Elliptical' Orbit: comparison with the ground state",AxesLabel->{"x","y"},PlotRange->All,PlotStyle->{Blue, Thick}]

Export["OganessonIon117SP.pdf",%]

(*End of Oganesson, Z=118, section*)

Unbibium, Ubb, Z=122

Clear[ntVal,nrVal,ωVal,aVal,rMin,rMax,Z,DeltaTheta, PeriodTheta, eVal,wNum]

(*Physical constants*)
Z=122; (*Atomic number of a hypotetical element Unbibium, Ubb*)

(*Quantum numbers*)
ntVal=1;
nrVal=1;

(*---Compute numeric values---*)

ωVal=ωFunc[ntVal]//N
ϵVal=ϵFunc[ntVal,nrVal]//N
aVal=aFunc[ntVal,nrVal]//N
{rMin=aVal*(1-ϵVal) //N, rMax=aVal*(1+ϵVal) //N}
DeltaTheta=2*Pi*((1/ωVal) -1) //N
PeriodTheta=2*Pi*(1/ωVal)  //N
eVal=eFunc[ntVal,nrVal]//N
{Z,Z/137 //N}
wNum=2*((1/ωVal) -1) //N

(*---Define relativistic radial orbit r(t) from Eq.(3.9)/(1.2)---*)
r[t_]:=(aVal (1-ϵVal^2))/(1+ϵVal Cos[ωVal t])

 PolarPlot[r[t],{t,0,0.68*PeriodTheta},PlotLabel->"Relativistic 'Elliptical' Orbit: 1.5 rotations",AxesLabel->{"x","y"},PlotRange->All,PlotStyle->{Blue,Thick}]
PolarPlot[r[t],{t,0.5*PeriodTheta,1.5*PeriodTheta},PlotLabel->"Relativistic Elliptical Orbit: 1 rotation, max->max",AxesLabel->{"x","y"},PlotRange->All,PlotStyle->{Blue, Thick}]
 PolarPlot[r[t],{t,0,0.889*PeriodTheta},PlotLabel->"Relativistic 'Elliptical' Orbit: 2 rotations",AxesLabel->{"x","y"},PlotRange->All,PlotStyle->{Blue,Thick}]
PolarPlot[r[t],{t,0,1.15*PeriodTheta},PlotLabel->"Relativistic 'Elliptical' Orbit: 2.5 rotations",AxesLabel->{"x","y"},PlotRange->All,PlotStyle->{Blue,Thick}]
PolarPlot[r[t],{t,0,1.369*PeriodTheta},PlotLabel->"Relativistic 'Elliptical' Orbit 3 rotations",AxesLabel->{"x","y"},PlotRange->All,PlotStyle->{Blue,Thick}]
PolarPlot[r[t],{t,0,1.6*PeriodTheta},PlotLabel->"Relativistic 'Elliptical' Orbit: 3.5 rotations",AxesLabel->{"x","y"},PlotRange->All,PlotStyle->{Blue,Thick}]

Export["UnbibiumIon121.pdf",%%]

(*Animation in Cartesian Coordinates*)Manipulate[ParametricPlot[{r[t] Cos[t],r[t] Sin[t]},{t,0,tMax},PlotStyle->{Blue,Thick},AxesLabel->{"x","y"},PlotLabel->Style["Relativistic 'Elliptical' Orbit(Animated, Cartesian)",14,Bold],PlotRange->{{-rMax,rMax},{-rMax,rMax}},(*Fixed range ensures full visibility from (Subscript[r, min]) to (Subscript[r, max])*)AspectRatio->1,PlotPoints->1000,PerformanceGoal->"Quality"],{tMax,0.01,4*0.725*PeriodTheta,Appearance->"Labeled"}]

Summary: 
In the above animation of the relativistic Kepler motion of an electron in hydrogen-like Unbibium, Ubb, when Z=122, the nucleus is situated  in the fixed focus at the origin.  For this quantum state, we choose (Subscript[n, r])=nr=1 and (Subscript[n, θ])=nt=1 in the fine structure formula (3.18). By (3.19)-(3.21), one gets: ω= 0.455419,  ϵ=0.949782, and a=0.0203534, in Bohr’s atomic units.  The ‘winding’ number:  2.39156, more than 2.
The perihelion and aphelion move along two concentric circles around the nucleus with radii:   (Subscript[r, min])=a (1-ϵ)=0.00102211 and (Subscript[r, max])=a (1+ϵ)=0.0396847, respectively.  In the animation, the geometrical loci of the successive perihelia and aphelia, the outer and inner envelopes of the orbit, see Figure 3, are not shown for simplicity. ‘Winding’ number is three.
CCW Rotation of the ellipse over one period is given by Δθ=2π(1/ω-1)=7.5133.  
At the point of the first self-intersection one gets, approximately,  r(10.002448573515926)=r(3.7940322175405243)=0.00234067  with  r(10.002448573515926)-r(3.7940322175405243)=0. At the second self-intersection, approximately, r[(6.9329615821129)=r(10.66)=0.0017562 with

Verification:

{ωVal, ϵVal, aVal}

{rMin, rMax}

{DeltaTheta,PeriodTheta,wNum}

{eVal,Z,ntVal,nrVal}

(* Start of evaluating of self-intersections: *)0.725*PeriodTheta

r[10.002448573515926`]

Solve[ωVal*(10.002448573515926`+t)/2-Pi==0,t]

r[3.794032217540524`]

10.002448573515926`/3.794032217540524`

FindRoot[r[2.636363636363636`*t]-r[t],{t,3.7}]

2.636363636363636`*3.7940322175405243`

r[10.002448573515926`]-r[3.7940322175405243`]

{r[10.002448573515926`],r[3.7940322175405243`]}

PolarPlot[r[t],{t,3.7940322175405243`,10.002448573515926`},PlotLabel->"Relativistic 'Elliptical' Orbit",AxesLabel->{"x","y"},PlotRange->All,PlotStyle->{Blue,Thick}]

Solve[ωVal*(10.66`+t)/2-2*Pi==0,t]

PolarPlot[r[t],{t,10.66`, 16.9329615821129`},PlotLabel->"Relativistic 'Elliptical' Orbit",AxesLabel->{"x","y"},PlotRange->All,PlotStyle->{Blue,Thick}]

r[16.932961582112902`]-r[10.660000000000002`]

16.932961582112902`/10.660000000000002`

FindRoot[r[1.588457934532167`*t]-r[t]==0,{t,10}]

1.588457934532167`*10.66`

r[16.9329615821129`]-r[10.660000000000002`]

16.9329615821129/10.660000000000002

FindRoot[r[1.5884579345321665`*t]-r[t]==0,{t,10.6}]

1.5884579345321665`*10.660000000000002`

r[16.9329615821129`]-r[10.660000000000002`]

{r[16.9329615821129`],r[10.660000000000002`]}

(*End of Unbibium, Ubb,  Z=122, section*)