Der anharmonische Oszillator

Clear["Global`*"];

Der Anharmonische Oszillator hat den Hamiltonian
H=((p^2)/(2m))+((m (ω^2)(x^2))/2)+(gx^4)
Im Folgenden wollen wir dessen Grundzustandsenergie unter der Nutzung unterschiedlicher Methoden berechnen. Achtung: Für die folgenden Rechnungen setzen wir ℏ=1.

Störungstheorie

Der spontan naheliegendste Ansatz ist natürlich Störungstheorie. Hierfür lösen wir zunächst den harmonischen Oszillator mit Hamiltonian
(Subscript[H, 0]) =((p^2)/(2m))+ ((((mω)^2)(x^2))/2)=-((ℏ^2)/(2m))Δ+((((mω)^2)(x^2))/2).
Die zugehörigen Eigenfunktionen sind dann:

ϕ[n_,m_,ω_,x_]:=(((((m ω)/(π ))))^(1/4))(1/(Sqrt[(2^n)n!]))HermiteH[n,(Sqrt[m ω])x](E^(-(1/2)m ω (x^2)));
Plot[Evaluate@Table[ϕ[n,1,1,x],{n,0,3}],{x,-5,5}]

Und als Eigenenergie ergeben sich:

E0[n_,ω_]:=ω(n+(1/2));

Die Energiekorrektur erster Ordnung Störungstheorie ist dann der Erwartungswert:
((ΔE)^((1)))=(∫_(-∞)^∞)(Subscript[ϕ, n]^2)(x)(gx^4)dx
Dieser lässt sich analytisch berechnen und ergibt sich zu:

E1[n_,m_,ω_,g_]:=(g/(4(m^2)(ω^2)))(6((n^2)+n)+3);

Und damit die Grundzustandsenergie des anharmonischen Oszillators:

EPT[n_,m_,ω_,g_]:=Refine[E0[n,ω]+E1[n,m,ω,g],ω∈Reals&&m∈Reals&&ω>0&&m>0];
Table[EPT[n,1,1,g],{n,0,5}]

Im Modell unabhängiger Teilchen

Configuration Interaction

Zunächst müssen wir das in Störungstheorie berechnete Matrixelement für zwei unterschiedliche Zustände (Subscript[ϕ, n])(x) und (Subscript[ϕ, k])(x) berechnen, also
g<n|(x^4)|k>.
Auch das lässt sich wieder analytisch lösen:

Enk[n_,k_,m_,ω_,g_]:=(g/(4(m^2)(ω^2)))((6((n^2)+n)+3)KroneckerDelta[n,k]+2(2n+3)(Sqrt[(n+2)(n+1)])KroneckerDelta[n,k-2]+2(2n-1)(Sqrt[n(n-1)])KroneckerDelta[n,k+2]+(Sqrt[(n+4)(n+3)(n+2)(n+1)])KroneckerDelta[n,k-4]+(Sqrt[n(n-1)(n-2)(n-3)])KroneckerDelta[n,k+4]);

Dann ist die Hamilton’sche Matrix:

H0[ω_,confs_]:=DiagonalMatrix[Table[E0[n,ω],{n,confs}]];
H1[m_,ω_,g_,confs_]:=Table[Enk[n,k,m,ω,g],{n,confs},{k,confs}];
Hmat[m_,ω_,g_,confs_]:=H0[ω,confs]+H1[m,ω,g,confs];

Diese müssen wir nur noch diagonalisieren und bekommen die Energieeigenwerte. Da wir aus der analytischen Herleitung der CI-Matrix wissen, dass immer Konfigurationen mit Inkrementen von zwei interagieren, erhöhen wir die Größe der Basis direkt um zwei Zustände.

ϵ[m_,ω_,g_,num_]:=Table[Eigenvalues[Hmat[m,ω,g,Table[l,{l,0,i}]],-1][[1]],{i,0,num,2}];
ϵrel[m_,ω_,g_,num_]:=(1/2)Divide[Differences[ϵ[m,ω,g,num]],Plus[ϵ[m,ω,g,num],RotateLeft[ϵ[m,ω,g,num]]][[;;-2]]];

Nun wollen wir einmal für feste Werte (Subscript[m, 0]), (Subscript[ω, 0])und (Subscript[g, 0]) die Grundzustandsenergie berechnen und die Konvergenz bezüglich der Basisgröße für unterschiedliche Teilchenmassen.

m0=1;
ω0=1;
g0=0.1;
num0=20;
Print[MatrixForm[ϵ[m0,ω0,g0,num0]]];
Print[MatrixForm[ϵrel[m0,ω0,g0,num0]]];
ListLinePlot[{ϵrel[m0,ω0,g0,num0],ϵrel[2m0,ω0,g0,num0],ϵrel[3m0,ω0,g0,num0]}]

Numerische Lösung

Um die Mächtigkeit unserer Methode zu beurteilen, wollen wir zum Schluss das Problem numerisch lösen.

m0=1;
ω0=1;
g0=0.1;
NDEigenvalues[-(1/(2m0))D[ψ[x],{x,2}]+(((m0 ((ω0)^2)(x^2))/2)+g0 (x^4))ψ[x],ψ[x],{x,-50,50},1,Method->{"SpatialDiscretization"->{"FiniteElement",{"MeshOptions"->{"MaxCellMeasure"->0.001}}},"Eigensystem"->{"Arnoldi",MaxIterations->40000}}][[1]]

Fazit: Hey, CI ist gar nicht mal so schlecht, schon mit etwa 15 Basiselementen.