Stark-Zeeman linear perturbation theory
A supplementary file for Alex Thorman’s 2018 ANU PhD thesis demonstrating the calculations contained in Chapter 2
The notebook should run in under 1 hour

First do Stark-Zeeman calculations without fine structure to derive analytic results

Load |n,l,m> hydrogen wavefunctions and possible combinations for n,l,m for the Balmer-α transition
Note ψ used for wavefunction and θ for polar angle. ψ is used for polar angle and θ for polarisation orientation in the thesis.

ψ[{n_,l_,m_}]=Sqrt[(2/(n a))^3 (n-l-1)!/(2n((n+l)!)^3)]Exp[-r/(n a)](2r/(n a))^l LaguerreL[n-l-1,2l+1,2r/(n a)](n+l)!SphericalHarmonicY[l,m,θ,ϕ];(*Hyrdogen |n,l,m> wavefunctions without spin Eq. 2.9*)
x=r Sin[θ]Cos[ϕ];
y=r Sin[θ]Sin[ϕ];
z=r Cos[θ];
nlm2={{2,0,0},{2,1,0},{2,1,1},{2,1,-1}};
nlm3={{3,0,0},{3,1,0},{3,2,0},{3,2,2},{3,2,-2},{3,1,1},{3,2,1},{3,1,-1},{3,2,-1}};

Evaluate individual Stark and Zeeman contributions to the (Subscript[H, n,αβ]) matrices
Could use some selection rules to simplify but brute force is fast enough, Lz and Lx operators applied in polar coordinates

starkz2=Table[Integrate[ψ[nlm2[[i]]]\[Conjugate]ψ[nlm2[[j]]]ϵ z/(3a)r^2 Sin[θ],{r,0,∞},{θ,0,π},{ϕ,0,2π},Assumptions->a>0],{i,1,4},{j,1,4}];(*Stark effect contribution to (Subscript[H, 2,αβ])*)
zeemanz2=Table[-I γ Integrate[ψ[nlm2[[i]]]\[Conjugate](D[ψ[nlm2[[j]]],ϕ])r^2 Sin[θ],{r,0,∞},{θ,0,π},{ϕ,0,2π},Assumptions->a>0],{i,1,4},{j,1,4}];(*Zeeman effect orbital angular momentum (Subscript[L, z ])contribution to (Subscript[H, 2,αβ])*)
zeemanx2=Table[I γ Integrate[ψ[nlm2[[i]]]\[Conjugate](Sin[ϕ]D[ψ[nlm2[[j]]],θ]+Cot[θ]Cos[ϕ]D[ψ[nlm2[[j]]],ϕ])r^2 Sin[θ],{r,0,∞},{θ,0,π},{ϕ,0,2π},Assumptions->a>0],{i,1,4},{j,1,4}];(*Zeeman effect orbital angular momentum (Subscript[L, x ])contribution to (Subscript[H, 2,αβ])*)
starkz3=Table[Integrate[ψ[nlm3[[i]]]\[Conjugate]ψ[nlm3[[j]]]ϵ z/(3a)r^2 Sin[θ],{r,0,∞},{θ,0,π},{ϕ,0,2π},Assumptions->a>0],{i,1,9},{j,1,9}];(*Stark effect contribution to (Subscript[H, 3,αβ])*)
zeemanz3=Table[-I γ Integrate[ψ[nlm3[[i]]]\[Conjugate](D[ψ[nlm3[[j]]],ϕ])r^2 Sin[θ],{r,0,∞},{θ,0,π},{ϕ,0,2π},Assumptions->a>0],{i,1,9},{j,1,9}];(*Zeeman effect orbital angular momentum (Subscript[L, z ])contribution to (Subscript[H, 3,αβ])*)
zeemanx3=Table[I γ Integrate[ψ[nlm3[[i]]]\[Conjugate](Sin[ϕ]D[ψ[nlm3[[j]]],θ]+Cot[θ]Cos[ϕ]D[ψ[nlm3[[j]]],ϕ])r^2 Sin[θ],{r,0,∞},{θ,0,π},{ϕ,0,2π},Assumptions->a>0],{i,1,9},{j,1,9}];(*Zeeman effect orbital angular momentum (Subscript[L, x ])contribution to (Subscript[H, 3,αβ])*)

Calculate eigenvalues (energies En,i in Eq. 2.12) and corresponding eigenvectors (mixing coefficients bn,iη in Eq. 2.11) for the n=2 matrix and then n=3
i.e. reproduce Tables 2.5 and 2.6 for |n,k,m>C states
Take limit that E and B are perpendicular (τ=0) to get |n,k,m>C orthogonalisation of degenerate Stark-Zeeman states. Also need to normalise. Because FullSimplify struggles to give nice result some substitutions are used to coax the solution into nice form. This was generated with Mathematica 11.0.1.0 and also works with 10.3.1.0, it is conceivable but unlikely that other versions may order the eigenvalues differently or apply the simplification differently in which case some adjustments would be needed to get the same formatting of the results

{energy2,vector2}=Eigensystem[starkz2 + (zeemanx2 Cos[τ]+zeemanz2 Sin[τ])];(*Combine all contributions to (Subscript[H, 2,αβ]) before finding eigenvalues and eigenvectors*)
order2={4,3,2,1};(*Required to get same ordering of Stark-Zeeman states in thesis. This may be version dependent and require adjustemnt.*)
energy2=Limit[energy2[[order2]],τ->0,Assumptions->ϵ>0&&γ>0](*n=2 Energy splitting*)
rescale2={-1,-1,-1,1};(*Scaling to get sign to agree with 2.5, a sign change is inconsequential*)
vector2=Table[FullSimplify[Simplify[Apart[Limit[rescale2[[i]] vector2[[order2[[i]]]]/Sqrt[vector2[[order2[[i]]]].vector2[[order2[[i]]]]],τ->0,Assumptions->ϵ>0&&γ>0]],Assumptions->ϵ>0&&γ>0]/.{((γ^2)+2 ϵ (ϵ+(Sqrt[(γ^2)+(ϵ^2)])))->(ϵ+(Sqrt[(γ^2)+(ϵ^2)]))^2,((γ^2)+2 ϵ (ϵ-(Sqrt[(γ^2)+(ϵ^2)])))->(ϵ-(Sqrt[(γ^2)+(ϵ^2)]))^2},Assumptions->ϵ>0&&γ>0],{i,1,4}](*Normalised n=2 mixing coefficients i.e. Table 2.5*)

{energy3,vector3}=Eigensystem[starkz3 + (zeemanx3 Cos[τ]+zeemanz3 Sin[τ])];(*Combine all contributions to (Subscript[H, 3,αβ]) before finding eigenvalues and eigenvectors*)
order3={5,4,1,3,2,9,7,6,8};(*Required to get same ordering of Stark-Zeeman states in thesis. This may be version dependent and require adjustemnt.*)
energy3=Limit[energy3[[order3]],τ->0,Assumptions->ϵ>0&&γ>0](*n=3 Energy splitting*)
rescale3={1,-1,Sin[τ],Sin[τ],Sin[τ],1,-1,-1,1};(*Rescale some eigenvectors to by factors of τ to avoid getting zero in the limit or infinity.*)
vector3=Table[FullSimplify[Limit[rescale3[[i]]vector3[[order3[[i]]]],τ->0,Assumptions->ϵ>0&&γ>0],Assumptions->ϵ>0&&γ>0],{i,1,9}];(*Not yet normalised eigenvectors, note this is the slowest calculation in the notebook*)
norm=FullSimplify[FullSimplify[Table[Sqrt[vector3[[i]].vector3[[i]]],{i,1,9}],Assumptions->ϵ>0&&γ>0]/.{(2 (γ^2)+3 ϵ (3 ϵ+(Sqrt[4 (γ^2)+9 (ϵ^2)])))->1/2(3 ϵ+(Sqrt[4 (γ^2)+9 (ϵ^2)]))^2,(2 (γ^2)-3 ϵ (-3 ϵ+(Sqrt[4 (γ^2)+9 (ϵ^2)])))->1/2((Sqrt[4 (γ^2)+9 (ϵ^2)])-3ϵ)^2,(2 (γ^2)+9 (ϵ^2)-3 ϵ (Sqrt[4 (γ^2)+9 (ϵ^2)]))->1/2((Sqrt[4 (γ^2)+9 (ϵ^2)])-3ϵ)^2},Assumptions->ϵ>0&&γ>0];
vector3=Table[FullSimplify[Apart[vector3[[i]]/norm[[i]]],Assumptions->ϵ>0&&γ>0]/.(Sqrt[6+((27 (ϵ^2))/(2 (γ^2)))])->(Sqrt[(3/(2 ))])(Sqrt[4 (γ^2)+9 (ϵ^2)])/γ,{i,1,9}](*Normalised n=3 mixing coefficients i.e. Table 2.6*)

Construct |nkm>L states (Eqs. 2.44 -2.47)
It is shown later the fine-structure Stark-Zeeman states converge to these states

vector2L=vector2;
vector2L[[{3,4}]]={(vector2[[3]]-vector2[[4]])/Sqrt[2],(vector2[[3]]+vector2[[4]])/Sqrt[2]};
vector3L=vector3;
vector3L[[{4,5}]]={(vector3[[4]]+vector3[[5]])/Sqrt[2],(vector3[[4]]-vector3[[5]])/Sqrt[2]};
vector3L[[{6,7}]]={(vector3[[6]]-vector3[[7]])/Sqrt[2],(vector3[[6]]+vector3[[7]])/Sqrt[2]};
vector3L[[{8,9}]]={(vector3[[8]]-vector3[[9]])/Sqrt[2],(vector3[[8]]+vector3[[9]])/Sqrt[2]};

Evaluate dipole vector for |nlm> states (BraKet in far right of Eq. 2.13)
Gives same result as in Howard Yuh 2005 MIT PhD thesis
For simplicity scale out Bohr Radius, can add back in later if needing an absolute intensity calibration
Again could use some selection rules to simplify but brute force is fast enough

rx=Table[Integrate[ψ[nlm2[[i]]]\[Conjugate]ψ[nlm3[[j]]]x r^2 Sin[θ],{r,0,∞},{θ,0,π},{ϕ,0,2π},Assumptions->a>0],{i,1,4},{j,1,9}]/a;
ry=Table[Integrate[ψ[nlm2[[i]]]\[Conjugate]ψ[nlm3[[j]]]y r^2 Sin[θ],{r,0,∞},{θ,0,π},{ϕ,0,2π},Assumptions->a>0],{i,1,4},{j,1,9}]/a;
rz=Table[Integrate[ψ[nlm2[[i]]]\[Conjugate]ψ[nlm3[[j]]]z r^2 Sin[θ],{r,0,∞},{θ,0,π},{ϕ,0,2π},Assumptions->a>0],{i,1,4},{j,1,9}]/a;

Evaluate dipole vector for |nkm(Subscript[>, C]) states (perform sum in Eq. 2.13)

Dx=FullSimplify[vector2\[Conjugate].rx.Transpose[vector3],Assumptions->ϵ>0&&γ>0];
Dy=FullSimplify[vector2\[Conjugate].ry.Transpose[vector3],Assumptions->ϵ>0&&γ>0];
Dz=FullSimplify[vector2\[Conjugate].rz.Transpose[vector3],Assumptions->ϵ>0&&γ>0];

Evaluate Stokes Vector (Eqs. 2.18-2.20)
Remember ψ is used for wavefunction so using θ for polar angle now

Ep1=Simplify[{-Sin[ϕ],Cos[ϕ],0}.{Dx,Dy,Dz}]; (*Eq 2.18*)
Ep2=Simplify[{-Cos[ϕ]Cos[θ],-Sin[ϕ]Cos[θ],Sin[θ]}.{Dx,Dy,Dz}];(*Eq 2.19*)
s0=Simplify[Apart[ComplexExpand[Abs[Ep1]^2+Abs[Ep2]^2]],Assumptions->ϵ>0&&γ>0&&0<θ<π&&-π<ϕ<π];(*Stokes vector components*)
s1=Simplify[Apart[ComplexExpand[Abs[Ep1]^2-Abs[Ep2]^2]],Assumptions->ϵ>0&&γ>0&&0<θ<π&&-π<ϕ<π];
s2=Simplify[Apart[ComplexExpand[2 Re[Ep1 Ep2\[Conjugate]]]],Assumptions->ϵ>0&&γ>0&&0<θ<π&&-π<ϕ<π];
s3=Simplify[Apart[ComplexExpand[2 Im[Ep1 Ep2\[Conjugate]]]],Assumptions->ϵ>0&&γ>0&&0<θ<π&&-π<ϕ<π];

Take series approximations to give some basic results

MatrixForm[Transpose[Simplify[Table[2(energy3[[i]]-energy2[[j]])/ϵ,{j,1,4},{i,1,9}]/.γ->0,Assumptions->ϵ>0]]](*Reproduce pure Stark energy splitting in units of ϵ/2*)

MatrixForm[Transpose[Simplify[5^14/(2^14*3^6)Limit[s0,γ->0,Assumptions->ϵ>0],Assumptions->ϵ>0]]](*Reproduce pure Stark relative intensities*)

MatrixForm[Transpose[FullSimplify[Normal[Series[ArcTan[s1,s2]/2,{γ,0,2}]],Assumptions->0<θ<π&&-π<ϕ<π&&ϵ>0]]]//Quiet(*Second order series for linear polarisation orientation, first order effects cancel for (SuperPlus[σ]) (SuperMinus[σ]) pairings*)

MatrixForm[Transpose[FullSimplify[Normal[Series[s3/s0,{γ,0,1}]],Assumptions->0<θ<π&&-π<ϕ<π&&ϵ>0]]](*First order series for circular polarisation fraction, zeroth order effects cancel for (SuperPlus[σ]) (SuperMinus[σ]) pairings*)

Now include spin and fine structure to get numeric results

Load |n,l,(Subscript[m, l]), (Subscript[m, s]) > hydrogen wavefunctions, possible quantum number combinations and fine structure variables
Now using s in name of each variable to indicate spin now included

ψs[{n_,l_,m_,ms_}]=Sqrt[(2/(n a))^3 (n-l-1)!/(2n((n+l)!)^3)]Exp[-r/(n a)](2r/(n a))^l LaguerreL[n-l-1,2l+1,2r/(n a)](n+l)!SphericalHarmonicY[l,m,θ,ϕ]{1/2+ms,1/2-ms};(*Hyrdogen |n,l,m> wavefunctions, now with spin as per Eq. 2.9*)
nlms2={{2,0,0,1/2},{2,0,0,-1/2},{2,1,0,1/2},{2,1,0,-1/2},{2,1,1,1/2},{2,1,1,-1/2},{2,1,-1,1/2},{2,1,-1,-1/2}};
nlms3={{3,0,0,1/2},{3,0,0,-1/2},{3,1,0,1/2},{3,1,0,-1/2},{3,2,0,1/2},{3,2,0,-1/2},{3,2,2,1/2},{3,2,2,-1/2},{3,2,-2,1/2},{3,2,-2,-1/2},{3,1,1,1/2},{3,1,1,-1/2},{3,2,1,1/2},{3,2,1,-1/2},{3,1,-1,1/2},{3,1,-1,-1/2},{3,2,-1,1/2},{3,2,-1,-1/2}};
En[2,0,1/2]=E2s1;(*Fine stucture variables*) 
En[2,1,1/2]=E2p1;
En[2,1,3/2]=E2p3;
En[3,0,1/2]=E3s1;
En[3,1,1/2]=E3p1;
En[3,1,3/2]=E3p3;
En[3,2,3/2]=E3d3;
En[3,2,5/2]=E3d5;

Evaluate individual Stark, Zeeman and Fine Structure contributions to the (Subscript[H, n,αβ]) matrices
Again could use some selection rules but brute force is fast enough, Lz and Lx operators applied in polar coordinates, Pauli Spin matrices used for Sz and Sx. Clebsch-Gordon coefficients used to get fine structure contribution as per Eq 2.10.

starkz2s=Table[Integrate[ψs[nlms2[[i]]]\[Conjugate].ψs[nlms2[[j]]]ϵ z/(3a)r^2 Sin[θ],{r,0,∞},{θ,0,π},{ϕ,0,2π},Assumptions->a>0],{i,1,8},{j,1,8}];(*Stark effect contribution to (Subscript[H, 2,αβ])*)
zeemanLz2s=Table[-I γ Integrate[ψs[nlms2[[i]]]\[Conjugate].(D[ψs[nlms2[[j]]],ϕ])r^2 Sin[θ],{r,0,∞},{θ,0,π},{ϕ,0,2π},Assumptions->a>0],{i,1,8},{j,1,8}];(*Zeeman effect orbital angular momentum (Subscript[L, z ])contribution to (Subscript[H, 2,αβ])*)
zeemanSz2s=Table[γ Integrate[ψs[nlms2[[i]]]\[Conjugate].{{1/2,0},{0,-1/2}}.ψs[nlms2[[j]]]r^2 Sin[θ],{r,0,∞},{θ,0,π},{ϕ,0,2π},Assumptions->a>0],{i,1,8},{j,1,8}];(*Zeeman effect spin angular momentum (Subscript[S, z ])contribution to (Subscript[H, 2,αβ])*)
zeemanLx2s=Table[I γ Integrate[ψs[nlms2[[i]]]\[Conjugate].(Sin[ϕ]D[ψs[nlms2[[j]]],θ]+Cot[θ]Cos[ϕ]D[ψs[nlms2[[j]]],ϕ])r^2 Sin[θ],{r,0,∞},{θ,0,π},{ϕ,0,2π},Assumptions->a>0],{i,1,8},{j,1,8}];
(*Zeeman effect orbital angular momentum (Subscript[L, x ])contribution to (Subscript[H, 2,αβ])*)
zeemanSx2s=Table[γ Integrate[ψs[nlms2[[i]]]\[Conjugate].{{0,1/2},{1/2,0}}.ψs[nlms2[[j]]]r^2 Sin[θ],{r,0,∞},{θ,0,π},{ϕ,0,2π},Assumptions->a>0],{i,1,8},{j,1,8}];
(*Zeeman effect spin angular momentum (Subscript[S, x ])contribution to (Subscript[H, 2,αβ])*)
fs2s=Table[Sum[ClebschGordan[{nlms2[[i,2]],nlms2[[i,3]]},{1/2,nlms2[[i,4]]},{j1,nlms2[[i,3]]+nlms2[[i,4]]}]ClebschGordan[{nlms2[[j,2]],nlms2[[j,3]]},{1/2,nlms2[[j,4]]},{j2,nlms2[[j,3]]+nlms2[[j,4]]}]KroneckerDelta[j1,j2]KroneckerDelta[nlms2[[i,2]],nlms2[[j,2]]]KroneckerDelta[nlms2[[i,3]]+nlms2[[i,4]],nlms2[[j,3]]+nlms2[[j,4]]]En[2,nlms2[[i,2]],j1],{j1,Max[Abs[nlms2[[i,3]]+nlms2[[i,4]]],nlms2[[i,2]]-1/2],nlms2[[i,2]]+1/2},{j2,Max[Abs[nlms2[[j,3]]+nlms2[[j,4]]],nlms2[[j,2]]-1/2],nlms2[[j,2]]+1/2}],{i,1,8},{j,1,8}];
(*Fine structure contribution to (Subscript[H, 2,αβ])*)
starkz3s=Table[Integrate[ψs[nlms3[[i]]]\[Conjugate].ψs[nlms3[[j]]]ϵ z/(3a)r^2 Sin[θ],{r,0,∞},{θ,0,π},{ϕ,0,2π},Assumptions->a>0],{i,1,18},{j,1,18}];(*Stark effect contribution to (Subscript[H, 3,αβ])*)
zeemanLz3s=Table[-I γ Integrate[ψs[nlms3[[i]]]\[Conjugate].(D[ψs[nlms3[[j]]],ϕ])r^2 Sin[θ],{r,0,∞},{θ,0,π},{ϕ,0,2π},Assumptions->a>0],{i,1,18},{j,1,18}];(*Zeeman effect orbital angular momentum (Subscript[L, z ])contribution to (Subscript[H, 3,αβ])*)
zeemanSz3s=Table[γ Integrate[ψs[nlms3[[i]]]\[Conjugate].{{1/2,0},{0,-1/2}}.ψs[nlms3[[j]]]r^2 Sin[θ],{r,0,∞},{θ,0,π},{ϕ,0,2π},Assumptions->a>0],{i,1,18},{j,1,18}];(*Zeeman effect spin angular momentum (Subscript[S, z ])contribution to (Subscript[H, 3,αβ])*)
zeemanLx3s=Table[I γ Integrate[ψs[nlms3[[i]]]\[Conjugate].(Sin[ϕ]D[ψs[nlms3[[j]]],θ]+Cot[θ]Cos[ϕ]D[ψs[nlms3[[j]]],ϕ])r^2 Sin[θ],{r,0,∞},{θ,0,π},{ϕ,0,2π},Assumptions->a>0],{i,1,18},{j,1,18}];
(*Zeeman effect orbital angular momentum (Subscript[L, x ])contribution to (Subscript[H, 3,αβ])*)
zeemanSx3s=Table[γ Integrate[ψs[nlms3[[i]]]\[Conjugate].{{0,1/2},{1/2,0}}.ψs[nlms3[[j]]]r^2 Sin[θ],{r,0,∞},{θ,0,π},{ϕ,0,2π},Assumptions->a>0],{i,1,18},{j,1,18}];
(*Zeeman effect spin angular momentum (Subscript[S, x ])contribution to (Subscript[H, 3,αβ])*)
fs3s=Table[Sum[ClebschGordan[{nlms3[[i,2]],nlms3[[i,3]]},{1/2,nlms3[[i,4]]},{j1,nlms3[[i,3]]+nlms3[[i,4]]}]ClebschGordan[{nlms3[[j,2]],nlms3[[j,3]]},{1/2,nlms3[[j,4]]},{j2,nlms3[[j,3]]+nlms3[[j,4]]}]KroneckerDelta[j1,j2]KroneckerDelta[nlms3[[i,2]],nlms3[[j,2]]]KroneckerDelta[nlms3[[i,3]]+nlms3[[i,4]],nlms3[[j,3]]+nlms3[[j,4]]]En[3,nlms3[[i,2]],j1],{j1,Max[Abs[nlms3[[i,3]]+nlms3[[i,4]]],nlms3[[i,2]]-1/2],nlms3[[i,2]]+1/2},{j2,Max[Abs[nlms3[[j,3]]+nlms3[[j,4]]],nlms3[[j,2]]-1/2],nlms3[[j,2]]+1/2}],{i,1,18},{j,1,18}];
(*Fine structure contribution to (Subscript[H, 3,αβ])*)

Generate Table 2.1

MatrixForm[starkz2s+Sin[τ](gl zeemanLz2s+gs zeemanSz2s)+Cos[τ](gl zeemanLx2s+gs zeemanSx2s)+fs2s]

Evaluate dipole vector for |nlmlms> states (BraKet in far right of Eq. 2.13)

rxs=Table[Integrate[ψs[nlms2[[i]]]\[Conjugate].ψs[nlms3[[j]]]x r^2 Sin[θ],{r,0,∞},{θ,0,π},{ϕ,0,2π},Assumptions->a>0],{i,1,8},{j,1,18}]/a;
rys=Table[Integrate[ψs[nlms2[[i]]]\[Conjugate].ψs[nlms3[[j]]]y r^2 Sin[θ],{r,0,∞},{θ,0,π},{ϕ,0,2π},Assumptions->a>0],{i,1,8},{j,1,18}]/a;
rzs=Table[Integrate[ψs[nlms2[[i]]]\[Conjugate].ψs[nlms3[[j]]]z r^2 Sin[θ],{r,0,∞},{θ,0,π},{ϕ,0,2π},Assumptions->a>0],{i,1,8},{j,1,18}]/a;

Load numerical values for beam voltage, magnetic field strength, injection angle, fine structure energies and view direction to get a solution
I was unable to get Eigensystem to give an analytic result when including the fine-structure symbolically so substitute in values to get a numerical answer

isotope=2;(*1 for hydrogen, 2 for deuterium*)
volt=80000 e;(*Beam voltage in eV multiplied by e to give SI*)
Bf=2;(*Magnetic field strength (Tesla)*)
ρ=π/3;(*Angle between Magnetic field and beam velocity*)
τ=0;(*Possible microscopic electric field effects ignored for now*)

mi=1.6726*^-27;
Ef=Sqrt[2 volt/(isotope*mi)]Bf Sin[ρ];(*Motional field strength*)
e=1.6022*^-19;
a=5.2918*^-11;
μB=9.274*^-24;
gl=1;
gs=2.002319;
ϵ=3e a Ef;
γ=μB Bf;
(*deuterium fine structure energies taken from NIST database https://physics.nist.gov/PhysRefData/ASD/levels_form.html, specifically A.E.Kramida, At.Data Nucl.Data Tables 96, 586–644 (2010)*)
If[isotope==1,
E∞=13.59843449e;
E2s1=10.19881044 e;
E2p1=10.19880606e;
E2p3=10.19885143e;
E3s1=12.08749486e;
E3p1=12.08749356e;
E3p3=12.08750700e;
E3d3=12.08750698e;
E3d5=12.08751146 e;
,If[isotope==2,
E∞=13.60213452e;
E2s1=10.20158545 e;
E2p1=10.20158107e;
E2p3=10.20162644e;
E3s1=12.09078378e;
E3p1=12.09078247e;
E3p3=12.09079592e;
E3d3=12.09079589 e;
E3d5=12.09080037 e;
,Print["Only has hydrogen and deuterium energies built in"]]
]

Get eigenvalues and eigenvectors of matrix for n=2 then n=3 levels
Numerical eigenvectors are already normalised
Note there is no longer any degeneracy (need to drill down some more significant figures for |3,0,±2> like states) therefore eigenvectors are unique

{energy2s,vector2s}=Eigensystem[starkz2s+Sin[τ](gl zeemanLz2s+gs zeemanSz2s)+Cos[τ](gl zeemanLx2s+gs zeemanSx2s)+fs2s];(*Combine all contributions to (Subscript[H, 2,αβ]) before finding eigenvalues and eigenvectors*)
order2s={1,2,7,8,3,5,4,6};(*Give same ordering as in table 2.10*)
energy2s=energy2s[[order2s]];
1000(energy2s-E∞(1-1/4))/e(*n=2 splitting in meV*)
vector2s=vector2s[[order2s]];
{energy3s,vector3s}=Eigensystem[starkz3s+Sin[τ](gl zeemanLz3s+gs zeemanSz3s)+Cos[τ](gl zeemanLx3s+gs zeemanSx3s)+fs3s];(*Combine all contributions to (Subscript[H, 3,αβ]) before finding eigenvalues and eigenvectors*)
order3s={1,2,17,18,9,12,8,10,7,11,3,6,4,5,14,15,13,16};(*Give same ordering as in table 2.10*)
energy3s=energy3s[[order3s]];
1000(energy3s-E∞(1-1/9))/e(*n=3 splitting in meV*)
vector3s=vector3s[[order3s]];

Show these fine-structure states are like |n,k,m>L
To compare states with and without spin I first decompose fine structure states into nlm states with spin aligned + and anti-aligned - with the magnetic field. Pauli spin matrix in the x direction gives that this is simply a sum and difference of pairs with spin in z direction.
Take dot product between |n,k,m>L and fine-structure states to show how similar they are.

aligned2s=Table[(vector2s[[i,2j-1]]+vector2s[[i,2j]])/Sqrt[2],{i,{1,3,5,7}},{j,1,4}];
anti2s=Table[(vector2s[[i,2j-1]]-vector2s[[i,2j]])/Sqrt[2],{i,{2,4,6,8}},{j,1,4}];
Table[Abs[vector2[[i]].aligned2s[[i]]],{i,1,4}](*Circular states with spin alligned fine-structue states. See the 20±1 like states are different*)
Table[Abs[vector2[[i]].anti2s[[i]]],{i,1,4}](*Circular states with spin anti-alligned fine-structue states. See the 20±1 like states are different*)
Table[Abs[vector2L[[i]].aligned2s[[i]]],{i,1,4}](*Circular states with spin alligned fine-structue states. All states are similar*)
Table[Abs[vector2L[[i]].anti2s[[i]]],{i,1,4}](*Circular states with spin anti-alligned fine-structue states. All states are similar*)

aligned3s=Table[(vector3s[[i,2j-1]]+vector3s[[i,2j]])/Sqrt[2],{i,{1,3,5,7,9,11,13,15,17}},{j,1,9}];
anti3s=Table[(vector3s[[i,2j-1]]-vector3s[[i,2j]])/Sqrt[2],{i,{2,4,6,8,10,12,14,16,18}},{j,1,9}];
Table[Abs[vector3[[i]].aligned3s[[i]]],{i,1,9}](*Circular states with spin alligned fine-structue states. See the 30±2 and 3±1±1 like states are different*)
Table[Abs[vector3[[i]].anti3s[[i]]],{i,1,9}](*Circular states with spin anti-alligned fine-structue states. See the 30±2 and 3±1±1 like states are different*)
Table[Abs[vector3L[[i]].aligned3s[[i]]],{i,1,9}](*Linear states with spin alligned fine-structue states. All states are similar*)
Table[Abs[vector3L[[i]].anti3s[[i]]],{i,1,9}](*Linear states with spin anti-alligned fine-structue states. All states are similar*)

Evaluate dipole vector for |nkm(Subscript[>, C]) states (full Eq 2.13)

Dxs=vector2s\[Conjugate].rxs.Transpose[vector3s];
Dys=vector2s\[Conjugate].rys.Transpose[vector3s];
Dzs=vector2s\[Conjugate].rzs.Transpose[vector3s];

Characterise transition energy, intensity and type
This characterisation uses metrics that would work for pure Stark states

MatrixForm[Transpose[Table[Round[20(energy3s[[i]]-energy2s[[j]]+E∞(1/9-1/4))/ϵ]/10.,{j,1,8},{i,1,18}]]]
(*Splitting in units of ϵ/2*)

Dmag2=Dxs^2+Abs[Dys]^2+Dzs^2;
MatrixForm[Transpose[Round[5^14/(2^14*3^6)Dmag2]]](*Relative intenisty*)

MatrixForm[Transpose[Round[Abs[100Dzs^2/Dmag2]]]](*% of being like a traditional pi transition, bright transitions are either close to 100% or 0% if sigma like*)

MatrixForm[Transpose[Round[100(Dxs^2-Abs[Dys]^2)/Dmag2]]](*% of being like a (σ^B) transition (if positive) or like a (σ^v) transition (if negative), bright transitions are either close to ±100% or 0% if pi like*)

MatrixForm[Transpose[Round[200(Dxs*I Dys)/Dmag2]]](*% of being like a traditional (SuperPlus[σ]) transition (if positive) or a traditional (SuperMinus[σ]) transition (if negative)*)

Generate components of Table 2.10
Take the maximum of Dx^2, Dy^2 or Dz^2 as the characteristic intensity

MatrixForm[Transpose[Round[5^14/(2^14*3^6)Dxs^2]]](*sigmaB like if bright*)
MatrixForm[Transpose[Round[-5^14/(2^14*3^6)Dys^2]]](*sigmav like if bright*)
MatrixForm[Transpose[Round[5^14/(2^14*3^6)Dzs^2]]](*pi like if bright*)

Now consider microscopic electric fields

Repeat calculations with microscopic electric field such that τ=1degree
Notice increased splitting between almost degenerate states
Try τ=0.001 deg instead, notice the |3,0,±2> like states give circular sigma emission but the others don’t

τ=1π/180;(*set τ to be no zero*)
{energy2s,vector2s}=Eigensystem[starkz2s+Sin[τ](gl zeemanLz2s+gs zeemanSz2s)+Cos[τ](gl zeemanLx2s+gs zeemanSx2s)+fs2s];(*Combine all contributions to (Subscript[H, 3,αβ]) before finding eigenvalues and eigenvectors*)
order2s={1,2,7,8,3,5,4,6};(*Give same ordering as in table 2.10*)
energy2s=energy2s[[order2s]];
1000(energy2s-E∞(1-1/4))/e(*n=2 splitting in meV*)
vector2s=vector2s[[order2s]];
{energy3s,vector3s}=Eigensystem[starkz3s+Sin[τ](gl zeemanLz3s+gs zeemanSz3s)+Cos[τ](gl zeemanLx3s+gs zeemanSx3s)+fs3s];(*Combine all contributions to (Subscript[H, 3,αβ]) before finding eigenvalues and eigenvectors*)
If[τ>0,order3s={1,2,17,18,9,12,7,10,8,11,3,5,4,6,13,15,14,16},order3s={1,2,17,18,9,12,8,11,7,10,4,6,3,5,14,16,13,15}];(*Give same ordering as in table 2.10, dependent on the sign of τ*)
energy3s=energy3s[[order3s]];
1000(energy3s-E∞(1-1/9))/e(*n=3 splitting in meV*)
vector3s=vector3s[[order3s]];

Show with nonzero τ states are like |n,k,m>C
Same calculation as earlier with τ=0 but now the result changes

aligned2s=Table[(vector2s[[i,2j-1]]+vector2s[[i,2j]])/Sqrt[2],{i,{1,3,5,7}},{j,1,4}];
anti2s=Table[(vector2s[[i,2j-1]]-vector2s[[i,2j]])/Sqrt[2],{i,{2,4,6,8}},{j,1,4}];
Table[Abs[vector2[[i]].aligned2s[[i]]],{i,1,4}](*Circular states with spin alligned fine-structue states. All states are similar*)
Table[Abs[vector2[[i]].anti2s[[i]]],{i,1,4}](*Circular states with spin anti-alligned fine-structue states. All states are similar*)
Table[Abs[vector2L[[i]].aligned2s[[i]]],{i,1,4}](*Linear states with spin alligned fine-structue states. See the 20±1 like states are different*)
Table[Abs[vector2L[[i]].anti2s[[i]]],{i,1,4}](*Linear states with spin anti-alligned fine-structue states. See the 20±1 like states are different*)

aligned3s=Table[(vector3s[[i,2j-1]]+vector3s[[i,2j]])/Sqrt[2],{i,{1,3,5,7,9,11,13,15,17}},{j,1,9}];
anti3s=Table[(vector3s[[i,2j-1]]-vector3s[[i,2j]])/Sqrt[2],{i,{2,4,6,8,10,12,14,16,18}},{j,1,9}];
Table[Abs[vector3[[i]].aligned3s[[i]]],{i,1,9}](*Circular states with spin alligned fine-structue states. All states are similar*)
Table[Abs[vector3[[i]].anti3s[[i]]],{i,1,9}](*Circular states with spin anti-alligned fine-structue states. All states are similar*)
Table[Abs[vector3L[[i]].aligned3s[[i]]],{i,1,9}](*Linear states with spin alligned fine-structue states. See the 30±2 and 3±1±1 like states are different*)
Table[Abs[vector3L[[i]].anti3s[[i]]],{i,1,9}](*Linear states with spin anti-alligned fine-structue states. See the 30±2 and 3±1±1 like states are different*)

Evaluate dipole vector and characterise transition energy, intensity and type
Same calculations as earlier when τ=0 but now the σ transitions are circular not linear v,B style

Dxs=vector2s\[Conjugate].rxs.Transpose[vector3s];
Dys=vector2s\[Conjugate].rys.Transpose[vector3s];
Dzs=vector2s\[Conjugate].rzs.Transpose[vector3s];
MatrixForm[Transpose[Table[Round[20(energy3s[[i]]-energy2s[[j]]+E∞(1/9-1/4))/ϵ]/10.,{j,1,8},{i,1,18}]]]
(*Splitting in units of ϵ/2*)

Dmag2=Dxs^2+Abs[Dys]^2+Dzs^2;
MatrixForm[Transpose[Round[5^14/(2^14*3^6)Dmag2]]](*Relative intenisty*)

MatrixForm[Transpose[Round[Abs[100Dzs^2/Dmag2]]]](*% of being like a traditional pi transition, bright transitions are either close to 100% or 0% if sigma like*)

MatrixForm[Transpose[Round[100(Dxs^2-Abs[Dys]^2)/Dmag2]]](*% of being like a (σ^B) transition (if positive) or like a (σ^v) transition (if negative), bright transitions are all close to 0%*)

MatrixForm[Transpose[Round[200(Dxs*I Dys)/Dmag2]]](*% of being like a traditional (SuperPlus[σ]) transition (if positive) or a traditional (SuperMinus[σ]) transition (if negative), bright transitions are either close to ±100% or 0% if pi like*)