Model powłokowy

Rozwiązujemy numerycznie równanie Schrödingera

Wstęp

Model kroplowy i inne podobne modele makroskopowe, mimo niewątpliwych sukcesów, nie są w stanie opisać wszystkich własności jądra. W szczególności nie da się w nich ująć zjawisk, w których uczestniczą pojedyncze nukleony (efekty jednocząstkowe).

Co więcej, możemy oczekiwać, że wszystkie własności jądra powinny dac się zrozumieć i opisać mikroskopowo, jako wynik oddziaływania między nukleonami. Model powłokowy jest pierwszym ważnym krokiem w tym kierunku.

W przybliżeniu nierelatywistycznym i przyjmując tylko siły dwuciałowe między nukleonami, hamiltonian układu A nukleonów jest:

H=T+V=(∑_(i=1)^A)((Subscript[p, i]^2)/(2m))+(∑_(i>k=1)^A)(Subscript[V, ik])((Subscript[r, i])-(Subscript[r, k]))

Ponieważ nie znamy dobrze potencjału sił jądrowych, to możemy uprościć problem wprowadzając pole średnie (potencjał jednocząstkowy):

H=(∑_(i=1)^A)[((Subscript[p, i]^2)/(2m))+(Subscript[U, i])(r)]+((∑_(i>k=1)^A)(Subscript[V, ik])((Subscript[r, i])-(Subscript[r, k]))-(Subscript[U, i])(r))

Jeśli dobrze dobierzemy jednocząstkowe pole średnie U, to wyraz w prawym nawiasie będzie mały i będzie traktowany jako zaburzenie (oddziaływanie resztkowe). W pierwszym przybliżeniu pomijamy ten wyraz i budujemy jądro z nukleonów, które obsadzają stany jednocząstkowe w potencjale U.

Istotę i funkcjonowanie modelu powłokowego można poznać przyjmując za U potencjał fenomenologiczny, mniej lub bardziej realistyczny. Takie podejście stosuje się też w wielu prostych zastosowaniach modelu. Lepszym rozwiązaniem jest wyznaczenie pola średniego metodą Hartree-Fock’a, ale to jest temat na inna okazję.

Potencjały fenomenologiczne

W praktyce wykorzystuje się najczęściej dwa potencjały fenomenologiczne: oscylatora harmonicznego i Woodsa-Saxona.

Potencjał 3-wym. oscylatora harmonicznego (HO):

(Subscript[U, HO])(r)=-(Subscript[V, 1])+(1/2)m (ω^2)(r^2).

Ma rozwiązania analityczne, co sprawia, że odgrywa ważną rolę, jako źródło dogodnej ortonormalnej bazy stanów. Pozwala łatwo wprowadzić deformację (model Nilssona). Często dodaje się człon ∼(l^2), aby “wypłaszczyć dno” i uczynić go bardziej realistycznym.

Potencjał Woodsa-Saxona (WS):

(Subscript[U, WS])(r)=((-(Subscript[V, 0]))/(exp[((r-R)/a)]+1)) .

Najbardziej realistyczny z potencjałów, wymaga jednak obliczeń numerycznych.

Wymienione potencjały opisują jedynie główny wkład sił jądrowych. Do tego trzeba dodać odpychanie kulombowskie między protonami, oraz dodatkowy wyraz opisujący jądrowe oddziaływanie spin-orbita.

Pełny, realistyczny potencjał modelu powłokowego (sferycznie symetryczny):

U(r)=(Subscript[U, WS])(r)+(Subscript[V, C])(r)+(Subscript[V, LS])(r) L·S .

Równanie Schrödingera

Równanie na stany stacjonarne cząstki o masie m w zadanym potencjale U (równanie Schrödingera niezależne od czasu):

((-(ℏ^2))/(2m))(∇^2)(Subscript[ψ, n])(r) +U(r)(Subscript[ψ, n])(r)=(Subscript[E, n]) (Subscript[ψ, n])(r).

W przypadku sferycznie symetrycznym można dokonać separacji zmiennych kątowych (θ, φ) od r podstawiając:

(Subscript[ψ, n])(r) = (Subscript[R, n])(r)·f(θ,φ) .

Rozwiązaniem części kątowej są oczywiście harmoniki sferyczne: f(θ,φ)=((Subscript[Y, l])^m)(θ,φ) , czyli
funkcje własne orbitalnego momentu pędu l . Część radialna funkcji falowej spełnia wtedy równanie:

[((-(ℏ^2))/(2m)) ((Subscript[∇, r])^2)+(((ℏ^2)l(l+1))/(2((mr)^2)))+U(r)](Subscript[R, nl])(r)=(Subscript[E, nl]) (Subscript[R, nl])(r),

czyli

((-(ℏ^2))/(2m)) (1/r)((d^2)/((dr)^2))(Subscript[rR, nl])(r)=[(Subscript[E, nl])-(((ℏ^2)l(l+1))/(2((mr)^2)))-U(r)] (Subscript[R, nl])(r) .

Wprowadzamy nową funkcję (Subscript[u, nl])(r)=(Subscript[rR, nl])(r)  i porządkujemy:

((d^2)/((dr)^2))(Subscript[u, nl])(r)=-((2m)/(ℏ^2))[(Subscript[E, nl])-(((ℏ^2)l(l+1))/(2((mr)^2)))-U(r)] (Subscript[u, nl])(r) .

To równanie będziemy chcieli rozwiązać numerycznie, żeby uzyskać funkcje radialne i wartości energii stanów jednocząstkowych.

Rozgrzewka: studnia prostokątna

Na początek rozważymy podobne równanie jednowymiarowe, lecz prostsze, żeby zapoznać się z jedną z metod numerycznych.

Rozwiązujemy:

((d^2)/((dx)^2))ψ(x)=-((2m)/(ℏ^2))[E-V(x)] ψ(x)

z prostokątną jamą potencjału o skończonej głębokości:

V(x)=({{{, ({{0, dla -a/2<x<a/2}, {(Subscript[V, 0]), wszędzie poza tym}})}})

Szukamy stanów związanych w tej jamie potencjału, czyli takich że E<(Subscript[V, 0 ]).
Zadanie to można rozwiązać analitycznie, ale dla nas pełni rolę treningową. 
Bardzo ładna ilustracja tego przykładu, ze zbioru Wolfram Demonstration Project, jest tutaj.

Metoda “strzelania” (shooting)

Podejście to jest zainspirowane wykładem Daniela V. Schroedera (link punkt 7.), na który zwrócił moją uwagę Krzysztof Miernik.

Warunki brzegowe

Wiemy, że funkcja falowa musi znikać poza studnią. Założymy więc, że jej wartość jest równa zero w pewnym punkcie na lewo od studni, a jej pochodna w tym punkcie jest prawie równa zero. Będziemy szukać (przez zgadywanie, stąd nazwa “shooting”) takiej wartości energii, żeby funkcja znikała także po prawej stronie studni.

Naturalne jednostki

W tego typu problemach należy wybrać rozsądne jednostki, aby algorytm numeryczny miał do czynienia z małymi liczbami. Nie tylko przyspiesza to rachunki, ale czasem je wręcz umożliwia. 
Tu naturalną jednostką długości jest a - szerokość studni. Kładziemy też m/(ℏ^2)=1. Jednostką energii będzie wtedy (ℏ^2)/(m (a^2)). W tych jednostkach równanie przybiera postać:

((d^2)/((dx)^2))ψ(x)=-2 [E-V(x)] ψ(x) .

Głębokość studni (Subscript[V, 0]) jest wolnym parametrem, od którego będą zależeć rozwiązania.

Wybór (Subscript[V, 0])

Aby wartości liczbowe miały, mimo wszystko, jakieś odniesienie do jąder atomowych, głębokość studni określamy nastepująco. Bezwzględną głębokość studni ustalamy na 50 MeV. Jako m bierzemy masę nukleonu, 940 MeV/(c^2). Wybieramy teraz “promień” studni, czyli a/2, w femtometrach. Ustala to jednostkę energii: eunit = 197^2/(940·(a^2)) MeV.  Wówczas (Subscript[V, 0])=50/eunit.

Manipulacja 1: prostokątna studnia potencjału (1D)

Uwaga: wartość energii oraz jej skok kontrolują przyciski (Button). Związane z nimi zmienne są globalne. Trzeba pamiętać, żeby w tej samej sesji unikać powtarzania ich nazw.

Manipulate[
Module[
{eunit,v0,v,sol,psi,norm,xb=2.0,
studnia=
Graphics[{Thick, Darker[Gray],
Line[{{-1.8,2},{-0.5,2},{-0.5,-2},{0.5,-2},{0.5,2},{1.8,2}}]}]
},
eunit=197.^2/(940.*(2*radius)^2);(* jednostka energii w MeV*)
v0=50./eunit;
(* definicja potencjału *)
v[x_]:= If[Abs[x]<0.5,0,v0];
(* centralne polecenie: numeryczne rozwiązanie równania S. *)
sol=NDSolve[{psi''[x]==-2*(energy1-v[x])psi[x],psi[-xb]==0,psi'[-xb]==0.0001},psi,{x,-xb,xb}];
(* normalizacja rozwiązania *)
norm=NIntegrate[(psi[x])^2/.sol,{x,-xb*0.6,xb*0.6}];
(* wykres otrzymanego rozwiązania *)
graph=Plot[psi[x]/Sqrt[norm]/.sol,{x,-xb,xb},
PlotRange->{{-2,2},{-2.2,2.2}},
Frame->True,
FrameLabel->{"x [2R]","ψ(x)"},
Axes->{True,False},
BaseStyle->{Black,FontFamily->"Verdana",FontSize->18},PlotStyle->{Red,Thickness[0.005]},
ImageSize->500];
Grid[{
{Row[{Style[" Jednostka energii = ",Black,FontSize->18],
       Style[NumberForm[eunit,3],Blue,FontFamily->"Verdana", FontSize->18],
Style[" MeV ",Black,FontSize->18]}],
Show[graph,studnia]
}
,
{Row[{Style["Głębokość studni = ",Black,FontSize->18],
       Style[NumberForm[v0,3],Blue,FontFamily->"Verdana", FontSize->18],
Style[" j.e. = 50 MeV",Black,FontSize->18]}],
SpanFromAbove
},
{Graphics[Text[Style[StringForm["    "]]],AspectRatio->0.03,ImageSize->300],
SpanFromAbove (* dummy, to make a space *)
}
,
{Button[" Energy Up ",
energy1=energy1+step1,
BaseStyle->18,
ImageMargins->5,
ImageSize->All]
Button[" Step Up ",
step1=If[step1==10.,10,step1*10.],
BaseStyle->18,
ImageMargins->5,
ImageSize->All],
SpanFromAbove
}
,
{Button["Energy Down",
energy1=energy1-step1,
BaseStyle->18,
ImageMargins->5,
ImageSize->All]
Button["Step Down",
step1=If[step1==0.00001,0.00001,step1/10.],
BaseStyle->18,
ImageMargins->5,
ImageSize->All],
SpanFromAbove
}
,
{Row[{Style[" Krok energii = ",Black,FontSize->18],
       Style[step1,Blue,FontFamily->"Verdana", FontSize->18],
Style[" j.e. ",Black,FontSize->18]}],
SpanFromAbove
},
{Graphics[Text[Style[StringForm["    "]]],AspectRatio->0.03,ImageSize->200],
SpanFromAbove (* dummy, to make a space *)
}
,
{Row[{Style["Energia stanu = ",Black,FontSize->18],
       Style[NumberForm[energy1,6],Blue,FontFamily->"Verdana", FontSize->18],
Style[" j.e. ",Black,FontSize->18]}],
SpanFromAbove
},
{Row[{Style["Energia stanu = ",Black,FontSize->18],
       Style[NumberForm[energy1*eunit,6],Red,FontFamily->"Verdana", FontSize->18],
Style[" MeV ",Black,FontSize->18]}],
SpanFromAbove
},
{Button["Reset",
energy1=2.;
step1=1.,
BaseStyle->18,
ImageMargins->5,
ImageSize->All],
SpanFromAbove
}
}
,Alignment->Top,ItemSize->Full]
],
(* Kontrolki *)
{{radius,2.0,Style["Promień jamy R [fm]",20]},1,5,0.1, Appearance->"Labeled", LabelStyle->Directive[Blue,20]},
FrameLabel->{None,None, Style["Równanie Schrödingera dla skończonej, prostokątnej jamy potencjału (1 wym.)",Blue,24], None},
Initialization:>{
step1=1.,
energy1=2.
}
]

Potencjał Woodsa - Saxona (WS)

Wracamy do “prawdziwych” jąder atomowych i szukamy funkcji radialnych i energii jednocząstkowych dla realistycznego potencjału Woodsa-Saxona. W tym celu musimy sprecyzować potencjał: U(r)=(Subscript[U, WS])(r)+(Subscript[V, C])(r)+(Subscript[V, LS])(r) L·S . Wybieramy parametryzację “klasyczną” pochodzącą z Książki (Bohr, Mottelson) opisaną też w podręczniku J. Suhonena “From Nucleons to Nuclei” (Springer 2007).

Wyraz główny, (Subscript[U, WS])(r)

(Subscript[U, WS])(r)=((-(Subscript[V, 0]))/(exp[((r-R)/a)]+1)) ,

(Subscript[V, 0])=(51±33 ((N-Z)/A))MeV,   R=(Subscript[r, 0])(A^(1/3)),   (Subscript[r, 0])=1.27 fm,   a=0.67 fm.

w powyżym znak “+“ odnosi się do protonów, a “-“ do neutronów.

Człon kulombowski, (Subscript[V, C])(r)

(Subscript[V, C])(r)=(((Ze)^2)/(4(Subscript[πε, 0])))({{{, ({{((3-(((r/R))^2))/(2R)), r≤R}, {1/r, r>R}})}}) ,

Człon ten występuje tylko dla protonów. W obliczeniach korzystamy z tego, że
 (e^2)/(4(Subscript[πε, 0])) = 1.44 MeV·fm .

Spin-orbita, (Subscript[V, LS])(r) L·S

(Subscript[V, LS])(r)=(v^((0)))(((((Subscript[r, 0])/ℏ)))^2)(1/r)[(d/dr)(1/(exp[((r-R)/a)]+1)) ].

Szukane funkcje falowe są funkcjami o dobrze określonym całkowitym i orbitalnym momencie pędu, j i l, przy czym j = l+1/2 lub j = l-1/2 w zależności od ustawienia spinu nukleonu względem orbitalnego momentu pędu. Operator L·S działając na taką funkcję daje:

L·S (Subscript[ψ, nlj])=(ℏ^2)(1/2)[j(j+1)-l(l+1)-(3/4)] (Subscript[ψ, nlj]) .

Zbierając, do równania musimy wstawić:

(Subscript[V, LS])(r)L·S=(v^((0)))(((Subscript[r, 0])^2)/(2a r))((exp[((r-R)/a)])/(((exp[((r-R)/a)]+1))^2)) [j(j+1)-l(l+1)-(3/4)],

(v^((0)))=0.44 ·(Subscript[V, 0]) .

Człon centryfugalny, (Subscript[V, L]):

(Subscript[V, L])(r)=(((ℏ^2)l(l+1))/(2m (r^2)))

formalnie, nie należy do potencjału, ale zmienia postać równania na funkcje radialne, więc ma znaczenie i warto zobaczyć jego wpływ na potencjał “efektywny”.

W literaturze występuje kilka innych parametryzacji potencjału Woodsa-Saxona. Jedna z nich, dość nowoczesna i chyba bardzo dobra, opisana jest w pracy: 
N. Schwierz, I. Wiedenhöver, and A. Volya, “Parameterization of the Woods-Saxon Potential for Shell-Model Calculations”, arxiv, 2007.

Manipulacja 2 : kształt potencjału WS

Manipulate[
Module[
{
(* Wartości parametrów klasycznych *)
V0=51.0 (* MeV *),
R0 = 1.27 (* fm *),
a = 0.67 (* fm *),
kappa = 33.0/51.,
lambda = 0.44,
(* zmienne lokalne *)
hc=197.3,
mred,R,RSO,n,j,V1,VN,VC,VL,VLS,Vtot
},
(* definicja potencjału Woodsa-Saxona według SWV *)
mred=940.0; (* (A-1)/A; *)
R=R0*A^(1/3);
RSO = R;
n=A-Z;
j=l+sign*(1/2);
V1=V0*(1+iso*kappa*(n-Z)/A);
VN[r_]:=-V1/(1+Exp[(r-R)/a]);
VC[r_]:=If[iso==1,1,0]*If[r<R,1.44*Z*(3-(r/R)^2)/(2*R),1.44*(Z-1)/r];
VL[r_]:= (hc^2/(2*mred))*l*(l+1)/r^2;
VLS[r_]:=If[LsOn,-lambda*V1*R0^2/(2*a*r)*  Exp[(r-RSO)/a]/((Exp[(r-RSO)/a]+1)^2)*(j*(j+1)-l*(l+1)-3/4),0];
Vtot[r_]:=VN[r]+VC[r]+VL[r]+VLS[r];
(* Plot command *)
gtot=Plot[Vtot[r],{r,0.2,21},PlotRange->{{0,21},{-60,50}},PlotStyle->{Black,Thickness[0.008]}];
gN=Plot[VN[r],{r,0.2,21}, PlotRange->{{0,21},{-60,50}},PlotStyle->{Blue,Thickness[0.005]}];
gC=Plot[VC[r],{r,0.2,21},PlotRange->{{0,21},{-60,50}},PlotStyle->{Red,Thickness[0.005]}];
gL=Plot[VL[r],{r,0.2,21},PlotRange->{{0,21},{-60,50}},PlotStyle->{Magenta,Thickness[0.005]}];
gLS=Plot[VLS[r],{r,0.2,21},PlotRange->{{0,21},{-60,50}},PlotStyle->{Orange,Thickness[0.005]}];
Show[gN,gC,gL,gLS,gtot,
AxesOrigin->{0.0,0.0},
BaseStyle->Directive[FontFamily->"Helvetica",FontSize->16,AxesStyle->Thick],
AxesLabel->{"r [fm]","Energia [MeV]"},AspectRatio->0.9,ImageSize->500]
],
(* Kontrolki *)
{{Z,71,Style["Z",18]},2,100,1,Appearance->"Labeled",LabelStyle->Directive[Red,18]},
(**)
{{A,151,Style["A",18]},Z,250,1,Appearance->"Labeled",LabelStyle->Directive[Black,18]},
(**)
Delimiter,
(**)
{{l,0,Style["L",18]},{0,1,2,3,4,5,6},ControlType->RadioButtonBar,Appearance->"Labeled",LabelStyle->Directive[Magenta,18]},
(**)
Delimiter,
(**)
{{iso,1,Style["Cząstki",18]},{1->Style[" protony   ",Red,16],-1->Style[" neutrony ",Blue,16]}},
Delimiter,
Row[{Control[{
{LsOn,False,Style["LS ",18]},{True,False}
}],
Control[{
{sign,1 ,""},
{1->Style[" j=l+1/2 ",17],
-1->Style[" j = l-1/2 ",17]},
Enabled->LsOn&&l>0
}]}],
(**)
ControlPlacement->Left,
ContinuousAction->False,
FrameLabel->{None,None, Style["Potencjał Woodsa-Saxona, część radialna, parametryzacja 'klasyczna'",Blue,24], None}
]

Potencjał WS, metoda “strzelania”

Przystępujemy do rozwiązywania radialnej części równania Schrödingera z potencjałem WS. Stosujemy metodę “strzelania”. Przypominamy pełne równanie:

((d^2)/((dr)^2))(Subscript[u, nlj])(r)=-((2m)/(ℏ^2))[(Subscript[E, nlj])-(((ℏ^2)l(l+1))/(2((mr)^2)))-U(r)] (Subscript[u, nlj])(r) .

Podobnie jak poprzednio obieramy “naturalne” jednostki. Narzucającą się jednostką długości jest promień jądra R.  Przyjmujemy (ℏ^2)/m=1, wtedy jednostką energii jest (((ℏc))^2)/((mc^2) (R^2)) . Za (mc^2) przyjmujemy 940 MeV (choć powinna być masa zredukowana!).

Manipulacja 3: metoda “strzelania”

Aby ułatwić zapisywanie wyników, tworzymy ich listę. Przycisk Zapisz dopisze do niej wartości aktualnych parametrów. Uwaga: Jeśli lista już zawiera jakieś dane, to wykonanie poniższej komendy skasuje wszystko co w niej było!

listaWynikow={}

Manipulate[
Module[{
(* Wartości parametrów "klasycznych"*)
V0=51.0 (* MeV *),
R0 = 1.27 (* fm *),
a0 = 0.67 (* fm *),
kappa = 33.0/51.,
lambda = 0.44,
hc=197.3,
(* pozostałe zmienne lokalne *)
A,mred,a,Vtot,norm,psi,g, norm2
},
(* definicja potencjału Woodsa-Saxona według SWV *)
(* przechodzimy do jednostek "naturalnych", R=1 *)
A=n+Z;
mred=940.0;(* (A-1)/A;*)
R=R0*A^(1/3);(*jednostka długości w fm*)
eunit=(hc)^2/(mred*R^2);(* jednostka energii w MeV *)
a=a0/R;
Vtot[r_]:=Module[{V1,VN,VC,VL,VLS},
V1=(V0/eunit)*(1+iso*kappa*(n-Z)/A);
VN=-V1/(1+Exp[(r-1)/a]);
VC=If[iso==1,1,0]*(1.44*Z/(eunit*R))*If[r<1,(3-r^2)/2,1.0/r];
VL=(1/2)*l*(l+1)/r^2;
VLS=If[LsOn,
If[l>0,
-(lambda/2)*(V1/eunit)*(R0/R)^2/(r*a)*  Exp[(r-1)/a]/((Exp[(r-1)/a]+1)^2)*
If[sign>0,l,-(l+1)],
0],
0
];
VN+VC+VL+VLS
];
vmin=First[FindMinimum[{Vtot[x],0.001≤x≤1},{x,0.5}]];
solu=NDSolve[{psi''[x]==-2*(energy2-Vtot[x])*psi[x],psi[0.0001]==0, psi'[0.0001]==1.0},psi,{x,0.0001,5}];
norm=NIntegrate[(psi[x])^2/.solu,{x,0.001,3}];
g=Plot[psi[x]/Sqrt[norm]/.solu,{x,0.001,4},
PlotRange->{{0,4},{-2,2}},
Frame->True,
Axes->{True,False},
FrameLabel->{"r/R","r·ψ(r)"},
BaseStyle->{Black,FontFamily->"Verdana",FontSize->18},PlotStyle->{Red,Thickness[0.005]},
ImageSize->600,
Epilog->{{Dashed,Line[{{0,0},{4,0}}]}}];
Row[{
Column[{
Row[{Style[" Jednostka długości, R = ",Black,FontSize->17],
       Style[NumberForm[R,3],Blue,FontFamily->"Verdana", FontSize->17],
Style[" fm",Black,FontSize->17]}],
Row[{Style[" Jednostka energii = ",Black,FontSize->17],
       Style[NumberForm[eunit,3],Blue,FontFamily->"Verdana", FontSize->17],
Style[" MeV ",Black,FontSize->17]}],
Row[{Style["Głębokość studni = ",Black,FontSize->17],
       Style[NumberForm[vmin,5],Blue,FontFamily->"Verdana", FontSize->17]}],
Row[{Button[" Energy Up ",
energy2=energy2+step2,
BaseStyle->17,
ImageMargins->5,
ImageSize->All]
Button[" Step Up ",
step2=If[step2==10.,10,step2*10.],
BaseStyle->17,
ImageMargins->5,
ImageSize->All]}],
Row[{Button["Energy Down",
energy2=energy2-step2,
BaseStyle->17,
ImageMargins->5,
ImageSize->All]
Button["Step Down",
step2=If[step2==0.0000001,0.0000001,step2/10.],
BaseStyle->17,
ImageMargins->5,
ImageSize->All]}],
Row[{Style[" Krok energii = ",Black,FontSize->17],
       Style[step2,Blue,FontFamily->"Verdana", FontSize->17],
Style[" j.e. ",Black,FontSize->17]}],
Row[{Style["Energia stanu = ",Black,FontSize->17],
       Style[NumberForm[energy2,6],Blue,FontFamily->"Verdana", FontSize->17],
Style[" j.e. ",Black,FontSize->18]}],
Row[{Style["Energia stanu = ",Black,FontSize->17],
       Style[NumberForm[energy2*eunit,6],Red,FontFamily->"Verdana", FontSize->17],
Style[" MeV ",Black,FontSize->17]}],
Row[{Button["Reset",
energy2=vmin;step2=1.,
BaseStyle->17,
ImageMargins->5,
ImageSize->All],
Spacer[20],
Button["Zapisz",
wavefun=Table[First[Flatten[psi[x]/.solu]],{x,0.01,3.01,0.02}];
norm2=R/50.*Sum[wavefun[[i]]^2,{i,1,Length[wavefun]}];
wavefun=Table[{(0.01+i-1)*R/50,wavefun[[i]]/Sqrt[norm2]},{i,1,Length[wavefun]}];
AppendTo[listaWynikow,{Z,n,n+Z,iso,l,LsOn,sign,eunit,energy2,eunit*energy2,wavefun}],
BaseStyle->17,
ImageMargins->5,
ImageSize->All],
Style[Length[listaWynikow],Black,FontFamily->"Verdana", FontSize->17]
}]
},Alignment->Center,Spacings->1],(* secondary controls *)
Show[g,ImageSize->{500,350}] (* plot panel *)
}]
], (* End of Module *)
(*Control panel*)
Panel[
Column[{
Row[{
Control[{
{Z,54,Style["Z",18]},2,100,1,Appearance->"Labeled",LabelStyle->Directive[Red,18]}],
Spacer[20],
Control[{
{n,82,Style["N",18]},2,150,1,Appearance->"Labeled",LabelStyle->Directive[Blue,18]}],
Spacer[20],
Style["A = ",18],Style[Dynamic[n+Z],18]
}],(* 1st row of controls end *)
Row[{
Control[{
{iso,-1,Style["Cząstki",18]},{1->Style[" protony   ",Red,16],-1->Style[" neutrony ",Blue,16]}
}],
Spacer[40],
Control[{
{l,0,Style["L",18]},{0,1,2,3,4,5,6},ControlType->RadioButtonBar,Appearance->"Labeled",LabelStyle->Directive[Magenta,18]
}],
Spacer[30],
Control[{
{LsOn,False,Style["LS ",18]},{True,False}
}],
Control[{
{sign,1 ,""},
Which[
l==0,{1->Style["  (Subscript[s, 1/2])   ",17],-1->Style["  (Subscript[s, 1/2])   ",17]},
l==1,{1->Style[" (Subscript[p, 3/2])    ",17],-1->Style[" \!\(\*SubscriptBox[\(p\), \(\(1/2\)\(\\ \\ \)\)]\)   ",17]},
l==2,{1->Style[" (Subscript[d, 5/2])    ",17],-1->Style[" (Subscript[d, 3/2])     ",17]},
l==3,{1->Style["  (Subscript[f, 7/2])    ",17],-1->Style["  (Subscript[f, 5/2])     ",17]},
l==4,{1->Style[" (Subscript[g, 9/2])    ",17],-1->Style[" (Subscript[g, 7/2])     ",17]},
l==5,{1->Style[" (Subscript[h, 11/2])  ",17],-1->Style[" (Subscript[h, 9/2])     ",17]},
l==6,{1->Style["  (Subscript[i, 13/2])  ",17],-1->Style["  (Subscript[i, 11/2])    ",17]}
],
ControlType->Setter,
Enabled->LsOn
}]
}](* 2nd row of controls end *)
},Alignment->Center,Spacings->1]
],(* Control panel end*)
FrameLabel->{None,None, Style["Równanie Schrödingera dla potencjału Woodsa-Saxona ",Blue,24], None},
Initialization:>{
step2=1.,
energy2=-20.
}
]

listaWynikow[[1]]

Zapisujemy listę wyników na dysk. Najpierw przechodzimy do właściwego folderu

SetDirectory["d:/Pfutzner/Teach/Mathematica4Nuclear/05_ModelPowlokowy"]

Directory[]

i robimy Export. Lepiej zachować tę instrukcję wykomentowaną, aby zabezpieczyć się przed niekontrolowanym wykonaniem jej, co nadpisze istniejący plik np. plikiem pustym.

Export["Wave_136Xe_p0f72.dat",listaWynikow[[3,11]],"Table"]

Wyniki przykładowe: 136Xe

Klika chwil klikania-strzelania pozwoliło uzyskać energie stanów jednocząstkowych neutronów i protonów dla 136Xe (Z=54, N=82). Aby ujawnić wpływ oddziaływania spin-orbita rachunki wykonano dwa razy: bez tego oddziaływania i z nim.

Wyniki bez członu spin-orbita:

neutronsWS={-38.81, -33.69, -27.62, -25.62, -20.76, -17.96, -13.25, -10.10, -9.045};

protonsWS={-35.82, -30.97,-24.96,-22.33,-18.03,-14.44,-10.31,-6.19,-4.71};

Wyniki z członem spin orbita :

neutronsWSLS={-38.81,-33.91,-33.27,-28.26,-26.72,-25.65,-22.00,-19.20,-18.34,-17.21,-15.25,-10.94,-10.89,-9.088,-8.83};

protonsWSLS={-35.82,-31.27,-30.42,-25.80,-23.83,-22.36,-19.61,-16.09,-14.91,-13.51,-12.84,-7.38,-7.26,-4.76,-4.58};

Dla pełnego porównania dodamy jeszcze energie z potencjału oscylatora harmonicznego (bez członu kulombowskiego i bez oddziaływania spin-orbita). Dane są one wyrażeniem:

((Subscript[E, nl])^HO)=-(Subscript[V, 1]) +(N+(3/2))ℏω =-(Subscript[V, 1]) +(2n+l+(3/2))ℏω,  n=0,1,2,3 ,...

Częstość oscylatora, ℏω, można oszacować (Blomqvist-Molinari):

ℏω=(45 (A^(-1/3))+25(A^(-2/3))) MeV.

Jako głębokość potencjału bierzemy (Subscript[V, 1])=50 MeV.

A=136;
homega=45.*A^(-1/3)-25.*A^(-2/3)

pureHO=Table[-50+(n+3/2)*7.805,{n,0,4}]

Rysujemy poziomy:

gHO=Graphics[Table[{Black,Thick,Line[{{20,#},{27,#}}]}]&/@pureHO];
gNWS=Graphics[Table[{Blue,Thick,Line[{{30,#},{37,#}}]}]&/@neutronsWS];
gNWSLS=Graphics[Table[{Red,Thick,Line[{{40,#},{47,#}}]}]&/@neutronsWSLS];
gPWS=Graphics[Table[{Blue,Thick,Line[{{10,#},{17,#}}]}]&/@protonsWS];
gPWSLS=Graphics[Table[{Red,Thick,Line[{{0,#},{7,#}}]}]&/@protonsWSLS];

I przedstawiamy na współnym rysunku, pracowicie opisując wszystko ręcznie :-(

Show[gHO,gNWS,gNWSLS,gPWS,gPWSLS,
Frame->{{True,True},{False,False}},
FrameLabel->{None,"Energia [MeV]"},FrameTicks->{{All,All},{None,None}},
PlotRange->{{-5,52},{-40.5,-1}},
LabelStyle->Directive[FontFamily->"Helvetica",FontSize->20],AspectRatio->0.6,PlotLabel->Style["Stany jednocząstkowe dla 136Xe",Darker[Blue],24],
ImageSize->1000,
Epilog->{
{Thin,Black, Dashed, Line[{{37,neutronsWS[[2]]},{40,neutronsWSLS[[2]]}}]},
{Thin,Black, Dashed, Line[{{37,neutronsWS[[2]]},{40,neutronsWSLS[[3]]}}]},
{Thin,Black, Dashed, Line[{{37,neutronsWS[[3]]},{40,neutronsWSLS[[4]]}}]},
{Thin,Black, Dashed, Line[{{37,neutronsWS[[3]]},{40,neutronsWSLS[[5]]}}]},
{Thin,Black, Dashed, Line[{{37,neutronsWS[[5]]},{40,neutronsWSLS[[7]]}}]},
{Thin,Black, Dashed, Line[{{37,neutronsWS[[5]]},{40,neutronsWSLS[[8]]}}]},
{Thin,Black, Dashed, Line[{{37,neutronsWS[[6]]},{40,neutronsWSLS[[9]]}}]},
{Thin,Black, Dashed, Line[{{37,neutronsWS[[6]]},{40,neutronsWSLS[[10]]}}]},
{Thin,Black, Dashed, Line[{{37,neutronsWS[[7]]},{40,neutronsWSLS[[11]]}}]},
{Thin,Black, Dashed, Line[{{37,neutronsWS[[7]]},{40,neutronsWSLS[[13]]}}]},
{Thin,Black, Dashed, Line[{{37,neutronsWS[[8]]},{40,neutronsWSLS[[12]]}}]},
{Thin,Black, Dashed, Line[{{37,neutronsWS[[8]]},{40,neutronsWSLS[[15]]}}]},
{Thin,Black, Dashed, Line[{{10,protonsWS[[2]]},{7,protonsWSLS[[2]]}}]},
{Thin,Black, Dashed, Line[{{10,protonsWS[[2]]},{7,protonsWSLS[[3]]}}]},
{Thin,Black, Dashed, Line[{{10,protonsWS[[3]]},{7,protonsWSLS[[4]]}}]},
{Thin,Black, Dashed, Line[{{10,protonsWS[[3]]},{7,protonsWSLS[[5]]}}]},
{Thin,Black, Dashed, Line[{{10,protonsWS[[5]]},{7,protonsWSLS[[7]]}}]},
{Thin,Black, Dashed, Line[{{10,protonsWS[[5]]},{7,protonsWSLS[[8]]}}]},
{Thin,Black, Dashed, Line[{{10,protonsWS[[6]]},{7,protonsWSLS[[9]]}}]},
{Thin,Black, Dashed, Line[{{10,protonsWS[[6]]},{7,protonsWSLS[[10]]}}]},
{Thin,Black, Dashed, Line[{{10,protonsWS[[7]]},{7,protonsWSLS[[11]]}}]},
{Thin,Black, Dashed, Line[{{10,protonsWS[[7]]},{7,protonsWSLS[[12]]}}]},
{Thin,Black, Dashed, Line[{{10,protonsWS[[8]]},{7,protonsWSLS[[13]]}}]},
{Thin,Black, Dashed, Line[{{10,protonsWS[[8]]},{7,protonsWSLS[[14]]}}]},
Text[Style[" protony ",Black,22],{8,-3}],
Text[Style[" neutrony ",Black,22],{38,-3}],
Text[Style[" WS+SO ",Darker[Red],20],{3.5,-40.3}],
Text[Style[" WS ",Darker[Blue],20],{13.5,-40.3}],
Text[Style[" HO ",Black,20],{23.5,-40.3}],
Text[Style[" WS ",Darker[Blue],20],{33.5,-40.3}],
Text[Style[" WS+SO ",Darker[Red],20],{43.5,-40.3}],
Text[Style[" N ",Black,17],{19.2,-40}],
Text[Style["0",Black,17],{19.2,-38.5}],
Text[Style["1",Black,17],{19.2,-30.5}],
Text[Style["2",Black,17],{19.2,-22.7}],
Text[Style["3",Black,17],{19.2,-14.9}],Text[Style["4",Black,17],{19.2,-7.01}],
Text[Style["0s",Black,17],{27,-38.5+0.8}],
Text[Style["0p",Black,17],{27,-30.5+0.8}],
Text[Style["1s 0d",Black,17],{27,-22.7+0.8}],
Text[Style["1p 0f",Black,17],{27,-14.9+0.8}],Text[Style["2s 1d 0g",Black,17],{27,-7.01+0.8}],
Text[Style["0s",Blue,17],{38,neutronsWS[[1]]}],Text[Style["0p",Blue,17],{38,neutronsWS[[2]]},Background->White],
Text[Style["0d",Blue,17],{38,neutronsWS[[3]]},Background->White],
Text[Style["1s",Blue,17],{38,neutronsWS[[4]]}],
Text[Style["0f",Blue,17],{38,neutronsWS[[5]]},Background->White],
Text[Style["1p",Blue,17],{38,neutronsWS[[6]]},Background->White],
Text[Style["0g",Blue,17],{38,neutronsWS[[7]]},Background->White],
Text[Style["1d",Blue,17],{38,neutronsWS[[8]]},Background->White],
Text[Style["2s",Blue,17],{38,neutronsWS[[9]]}],
Text[Style["0s",Blue,17],{9.2,protonsWS[[1]]}],Text[Style["0p",Blue,17],{9.2,protonsWS[[2]]},Background->White],
Text[Style["0d",Blue,17],{9.2,protonsWS[[3]]},Background->White],
Text[Style["1s",Blue,17],{9.2,protonsWS[[4]]}],
Text[Style["0f",Blue,17],{9.2,protonsWS[[5]]},Background->White],
Text[Style["1p",Blue,17],{9.2,protonsWS[[6]]},Background->White],
Text[Style["0g",Blue,17],{9.2,protonsWS[[7]]},Background->White],
Text[Style["1d",Blue,17],{9.2,protonsWS[[8]]},Background->White],
Text[Style["2s",Blue,17],{9.2,protonsWS[[9]]}],
Text[Style["0(Subscript[s, 1/2])",Red,17],{48.6,neutronsWSLS[[1]]}],
Text[Style["0(Subscript[p, 3/2])",Red,17],{48.6,neutronsWSLS[[2]]-0.2}],
Text[Style["0(Subscript[p, 1/2])",Red,17],{48.6,neutronsWSLS[[3]]+0.2}],
Text[Style["0(Subscript[d, 5/2])",Red,17],{48.6,neutronsWSLS[[4]]}],
Text[Style["0(Subscript[d, 3/2])",Red,17],{48.6,neutronsWSLS[[5]]}],
Text[Style["1(Subscript[s, 1/2])",Red,17],{48.6,neutronsWSLS[[6]]+0.1}],
Text[Style["0(Subscript[f, 7/2])",Red,17],{48.6,neutronsWSLS[[7]]}],
Text[Style["0(Subscript[f, 5/2])",Red,17],{48.6,neutronsWSLS[[8]]-0.3}],
Text[Style["1(Subscript[p, 3/2])",Red,17],{48.6,neutronsWSLS[[9]]+0.2}],
Text[Style["1(Subscript[p, 1/2])",Red,17],{48.6,neutronsWSLS[[10]]}],
Text[Style["0(Subscript[g, 9/2])",Red,17],{48.6,neutronsWSLS[[11]]}],
Text[Style["0(Subscript[g, 7/2])",Red,17],{48.6,neutronsWSLS[[13]]+0.2}],
Text[Style["1(Subscript[d, 5/2])",Red,17],{48.6,neutronsWSLS[[12]]-0.5}],
Text[Style["2(Subscript[s, 1/2])",Red,17],{48.6,neutronsWSLS[[14]]-0.2}],
Text[Style["1(Subscript[d, 3/2])",Red,17],{48.6,neutronsWSLS[[15]]+0.5}],
Text[Style["0(Subscript[s, 1/2])",Red,17],{-1.4,protonsWSLS[[1]]}],
Text[Style["0(Subscript[p, 3/2])",Red,17],{-1.4,protonsWSLS[[2]]-0.2}],
Text[Style["0(Subscript[p, 1/2])",Red,17],{-1.4,protonsWSLS[[3]]+0.2}],
Text[Style["0(Subscript[d, 5/2])",Red,17],{-1.4,protonsWSLS[[4]]}],
Text[Style["0(Subscript[d, 3/2])",Red,17],{-1.4,protonsWSLS[[5]]}],
Text[Style["1(Subscript[s, 1/2])",Red,17],{-1.4,protonsWSLS[[6]]+0.1}],
Text[Style["0(Subscript[f, 7/2])",Red,17],{-1.4,protonsWSLS[[7]]}],
Text[Style["0(Subscript[f, 5/2])",Red,17],{-1.4,protonsWSLS[[8]]-0.3}],
Text[Style["1(Subscript[p, 3/2])",Red,17],{-1.4,protonsWSLS[[9]]+0.2}],
Text[Style["1(Subscript[p, 1/2])",Red,17],{-1.4,protonsWSLS[[10]]}],
Text[Style["0(Subscript[g, 9/2])",Red,17],{-1.4,protonsWSLS[[11]]}],
Text[Style["0(Subscript[g, 7/2])",Red,17],{-1.4,protonsWSLS[[12]]-0.5}],
Text[Style["1(Subscript[d, 5/2])",Red,17],{-1.4,protonsWSLS[[13]]+0.2}],
Text[Style["2(Subscript[s, 1/2])",Red,17],{-1.4,protonsWSLS[[15]]+0.5}],
Text[Style["1(Subscript[d, 3/2])",Red,17],{-1.4,protonsWSLS[[14]]-0.4}],
Text[Style["2",Black,20],{43.5,-36.5}],
Text[Style["8",Black,20],{43.5,-31}],
Text[Style["20",Black,20],{43.5,-24}],
Text[Style["28",Black,20],{43.5,-20.8}],
Text[Style["50",Black,20],{43.5,-13.5}],
Text[Style["2",Black,20],{3.5,-34}],
Text[Style["8",Black,20],{3.5,-29}],
Text[Style["20",Black,20],{3.5,-21.5}],
Text[Style["28",Black,20],{3.5,-18.2}],
Text[Style["50",Black,20],{3.5,-10.5}]
}
]

Rozwiązanie w bazie oscylatora harmonicznego

Stany własne i energie jednocząstowe równania S. z potencjałem WS można znacznie szybciej i wygodniej, choć mniej dokładnie, wyznaczyć inną, bardzo ogólną i ważną metodą.

Przedstawiamy szukane rozwiązania jako kombinację liniową znanych funkcji, które tworzą układ zupełny, ortonormalny. Jako te bazę wybieramy rozwiązania oscylatora harmonicznego:

(Subscript[R, nlj])(r) =(∑_ν)((Subscript[A, ν])^((nlj)))(Subscript[g, νl])(r),  (∑_ν)(([((Subscript[A, ν])^((nlj)))])^2)=1

Tworzymy macierz, której elementy są zadane przez:

(Subscript[M, ν'ν])=⟨ν' |(Subscript[h, lj])|ν⟩ ≡(∫_0^∞)(Subscript[g, ν'l])(r)(Subscript[h, lj ])(Subscript[g, νl])(r) (r^2)dr ,

gdzie (Subscript[h, lj]) jest hamiltonianem, dla którego szukamy rozwiązań. Zgodnie z wcześniejszym, jego postać to:

(Subscript[h, lj])=((-(ℏ^2))/(2m))[ ((Subscript[∇, r])^2)-((l(l+1))/(r^2))]+(Subscript[U, WS])(r)+(Subscript[V, C])(r)+(Subscript[V, LS])(r) ·((ℏ^2)/2)[j(j+1)-l(l+1)-(3/4)] .

Aby obliczyć elementy macierzowe musimy “zadziałać” h na funkcje (Subscript[g, nl]). Wypiszmy najpierw r-nie oscylatora harmonicznego, które spełniają te funkcje:

[((-(ℏ^2))/(2m))[ ((Subscript[∇, r])^2)-((l(l+1))/(r^2))]-(Subscript[V, 1])+(1/2)m (ω^2)(r^2)](Subscript[g, nl])(r)=((Subscript[E, nl])^HO)(Subscript[g, nl])(r)
=[-(Subscript[V, 1]) +(2n+l+(3/2))ℏω](Subscript[g, nl])(r)

Wynik działania pierwszego wyrazu (Subscript[h, lj]) na funkcje (Subscript[g, nl]) możemy “wziąć” z równania oscylatora! W ten sposób unikamy dwukrotnego różniczkowania przy obliczaniu elementów macierzowych. Stosując ten trick dostajemy nastepujący wzór na elementy macierzy M:

⟨ν' |(Subscript[h, lj])|ν⟩ =(∫_0^∞)(Subscript[g, ν'l])(r)(Subscript[g, νl])(r)
[((ℏω)/2)(4ν+2l+3-((r^2)/(b^2)))+(Subscript[U, WS])(r)+(Subscript[V, C])(r)+(Subscript[V, LS])(r) ·((ℏ^2)/2)[j(j+1)-l(l+1)-(3/4)]] (r^2)dr ,

gdzie

(b^2)=((((ℏc))^2)/(((mc)^2)·ℏω)) ,     ℏω=(45 (A^(-1/3))-25(A^(-2/3))) MeV.

Całki te łatwo obliczymy numerycznie. Gdy już utworzymy macierz M, to rozwiązujemy jej problem własny (diagonalizujemy). Wartości własne są szukanymi energiami, a wektory własne to współczynniki rozkładu (Subscript[A, ν]), które zadają postać funkcji falowej. Przybliżenie tej metody polega na tym, że musimy zadać skończony wymiar tej macierzy (przy nieskończonym wynik byłby dokładny!), a to jest ograniczone mocą obliczeniową.

Funkcje falowe oscylatora harmonicznego

Musimy jeszcze jawnie wypisać postać funkcji (Subscript[g, nl])(r). Oto ona:

(Subscript[g, nl])(r) =(Sqrt[((2n!)/((b^3)Γ(n+l+(3/2))))])((((r/b)))^l)(e^(-(r^2)/2(b^2)))(Subscript[L, n]^((l+(1/2))))((r^2)/(b^2)) ,

gdzie (Subscript[L, n]^k) to stowarzyszone wielomiany Laguerre’a. Mathematica je zna:

LaguerreL[1,l+1/2,x]

Zobaczmy jak te funkcje (Subscript[g, nl])(r) wyglądają.

gnl[A_,n_,l_,r_]:=Module[{b,homega},
homega=45.*A^(-1/3)-25.*A^(-2/3);
b=197.3/Sqrt[940.0*homega];
Sqrt[(2*n!/(b^3*Gamma[n+l+3/2]))]*(r/b)^l*Exp[-r^2/(2*b^2)]*LaguerreL[n,l+1/2,r^2/b^2]
]

Rysujemy zestaw dla n=0 i kilku wartości l:

Plot[{
gnl[16,0,0,r]*r,
gnl[16,0,1,r]*r,
gnl[16,0,2,r]*r,
gnl[16,0,3,r]*r,
gnl[16,0,4,r]*r
},{r,0,10},
PlotRange->{{0,10},{0,0.8}},
Frame->True,
FrameLabel->{"r [fm]","r ψ(r)"},
BaseStyle->{Black,FontFamily->"Verdana",FontSize->18},
PlotStyle->{
Directive[Black, Thick],
Directive[Blue, Thick],
Directive[Magenta, Thick],
Directive[Orange, Thick],
Directive[Red, Thick]
},
PlotLegends->{"s","p","d","f","g"},
ImageSize->600,
Epilog->{Text[Style["\!\(\*SuperscriptBox[\(\), \(16\)]\)O  n = 0",Black,24],{7,0.6}]}
]

Oraz zestaw dla l=0 i kilku wartości n:

Plot[{
gnl[16,0,0,r]*r,
gnl[16,1,0,r]*r,
gnl[16,2,0,r]*r,
gnl[16,3,0,r]*r
},{r,0,10},
PlotRange->{{0,10},{-0.8,0.8}},
Frame->True,
FrameLabel->{"r [fm]","r ψ(r)"},
BaseStyle->{Black,FontFamily->"Verdana",FontSize->18},
PlotStyle->{
Directive[Black, Thick],
Directive[Blue, Thick],
Directive[Magenta, Thick],
Directive[Red, Thick]
},
PlotLegends->{"0","1","2","3"},
ImageSize->600,
Epilog->{Text[Style["(^16)O l = 0",Black,24],{7,0.6}]}
]

Powinniśmy też sprawdzić regułę ortonormalności względem n:

(∫_0^∞)(Subscript[g, n'l])(r)(Subscript[g, nl])(r) (r^2)dr =(Subscript[δ, nn']).

NIntegrate[gnl[16,1,2,x]*gnl[16,1,2,x]*x^2,{x,0,Infinity}]

Manipulacja 4: radialne funkcje HO

Dobrze też zobaczyć przebieg funkcji (Subscript[g, nl])(r) dla dowolnych parametrów, w tym A. Warto też przypomnieć, że metodą “strzelania” wyznaczaliśmy nie ostateczną funkcję radialną, tylko (Subscript[u, nlj])(r)=r (Subscript[R, nlj])(r) . Dla lepszego porównania i uwidocznienia różnicy rysujemy oba warianty: (Subscript[g, nl]) i r (Subscript[g, nl]) .

Manipulate[
Plot[gnl[A,n,l,r]*If[fselect,r,1],{r,0,20},
PlotRange->{{0,20},All},
Frame->True,
Axes->{True,False},
GridLines->Automatic,
GridLinesStyle->Dashed,
Filling->Axis,
FrameLabel->{"r [fm]","f(r)"},
BaseStyle->{Black,FontFamily->"Verdana",FontSize->20},PlotStyle->{If[fselect,Red,Blue],Thickness[0.005]},
ImageSize->600],
{{fselect,True ,Style["f(r) =",20]},{True->Style[" r·(Subscript[g, nl])(r) ",Red,20],False->Style[" (Subscript[g, nl])(r) ",Blue,20]},ControlType->SetterBar,Background->Gray,LabelStyle->Directive[Magenta,18]},
Spacer[20],
{{A,40,Style["A",20]},10,240,1, Appearance->"Labeled", LabelStyle->Directive[Blue,20]},
Spacer[30],
{{n,0,Style["n",20]},{0->Style[" 0 ",20],1->Style[" 1 ",20],2->Style[" 2 ",20],3->Style[" 3 ",20],4->Style[" 4 ",20]},ControlType->SetterBar,Background->Orange,LabelStyle->Directive[Magenta,18]},
Spacer[40],
{{l,0,Style["l",20]},{0,1,2,3,4,5,6},ControlType->RadioButtonBar,Appearance->"Labeled",LabelStyle->Directive[Magenta,20]},
SaveDefinitions->True,
FrameLabel->{None,None, Style["Radialna część funkcji falowej oscylatora harmonicznego",Blue,24], None}
]

Rozwinięcie w bazie HO

Budujemy procedurę, która dla zadanych wartości N, Z, j i l oraz rodzaju cząstek (neutrony/protony) oblicza elementy macierzy M (mtx) , a następnie rozwiązuje problem własny.

matrixHO[n_,Z_,iso_,j_,l_]:=Module[
(* liczaba neutronów, protonów, iso: +1=protony, -1=neutrony *)
{
(* Wartości parametrów potencjału WS (Suhonen) *)
V0=51.0 (* MeV *),
R0 = 1.27 (* fm *),
(* R0SO = 1.16 *) (* fm *),
a = 0.67 (* fm *),
kappa = 33.0/51.0,
lambda = 0.44,
hc=197.3,
A,R,RSO, mred,isocor,homega,b,V1,VWS,hWS,VC,VLS,gnl,mtx
},
A=n+Z;
mred=940.0 (* (A-1)/A*931.49 *);
R=R0*A^(1/3);
RSO =R (* R0SO*A^(1/3) *);
(* isocor=If[n==Z,3,If[n>Z,iso*(n-Z+1)+2,iso*(n-Z-1)+2]];*)
homega=45.*A^(-1/3)-25.*A^(-2/3);
b=hc/Sqrt[mred*homega];
(* potencjał WS bez centryfugi *)
hWS[r_]:=Module[{},
V1=V0*(1+iso*kappa*(n-Z)/A);
VWS=-V1/(1+Exp[(r-R)/a]);
VC=If[iso==1,1,0]*1.44*Z*If[r<R,(3-(r/R)^2)/(2*R),1.0/r];
VLS=-lambda*V1*R0^2/(2*a*r)*  Exp[(r-RSO)/a]/((Exp[(r-RSO)/a]+1)^2)*(j*(j+1)-l*(l+1)-3/4);
VWS+VC+VLS
];
(* rozwiązanie dla HO *)
gnl[nu_,r_]:=Module[{},
Sqrt[(2*nu!/(b^3*Gamma[nu+l+3/2]))]*(r/b)^l*Exp[-r^2/(2*b^2)]*LaguerreL[nu,l+1/2,r^2/b^2]
];
(* obliczenie elementu macierzowego *)
mtx=Table[NIntegrate[gnl[f,r]*gnl[i,r]*(homega/2 * (4*i+2*l+3-(r/b)^2)+hWS[r])*r^2,{r,0,Infinity}],{i,0,5},{f,0,5}];
Eigensystem[mtx]
]

sol=matrixHO[82,54,1,7/2,4]
(* Select[Sort[sol[[1]]],#<2.0&]*)

sol[[1]]

sol[[2,5]]

Porównanie metody “strzelania” z diagonalizacją “HO” dla 136Xe: energie stanów

neutronsWSLS={-38.81,-33.91,-33.27,-28.26,-26.72,-25.65,-22.00,-19.20,-18.34,-17.21,-15.25,-10.94,-10.89,-9.088,-8.83};

neutronsHO={-38.81,-33.91,-33.28,-28.24,-26.74,-25.65,-21.96,-19.24,-18.32,-17.22,-15.19,-10.95,-10.91,-9.08,-8.86};

protonsWSLS={-35.82,-31.27,-30.42,-25.80,-23.83,-22.36,-19.61,-16.09,-14.91,-13.51,-12.84,-7.38,-7.26,-4.76,-4.58};

protonsHO={-35.82,-31.26,-30.43,-25.77,-23.86,-22.35,-19.56,-16.14,-14.87,-13.52,-12.76,-7.45,-7.21,-4.74,-4.61};

g1PHO=Graphics[Table[{Blue,Thick,Line[{{0,#},{7,#}}]}]&/@protonsHO];
g1PWS=Graphics[Table[{Red,Thick,Line[{{10,#},{17,#}}]}]&/@protonsWSLS];
g1NWS=Graphics[Table[{Red,Thick,Line[{{20,#},{27,#}}]}]&/@neutronsWSLS];
g1NHO=Graphics[Table[{Blue,Thick,Line[{{30,#},{37,#}}]}]&/@neutronsHO];

Show[g1PHO,g1PWS,g1NWS,g1NHO,
Frame->{{True,True},{False,False}},
FrameLabel->{None,"Energia [MeV]"},FrameTicks->{{All,All},{None,None}},
PlotRange->{{-5,42},{-40.9,-1}},
LabelStyle->Directive[FontFamily->"Helvetica",FontSize->20],AspectRatio->0.6,PlotLabel->Style["Stany jednocząstkowe dla 136Xe",Darker[Blue],24],
ImageSize->1000,
Epilog->{
Text[Style[" protony ",Black,22],{8,-3}],
Text[Style[" neutrony ",Black,22],{28,-3}],
Text[Style[" HO ",Darker[Blue],20],{3.5,-40.3}],
Text[Style[" Shooting ",Darker[Red],20],{13.5,-40.3}],
Text[Style[" Shooting ",Darker[Red],20],{23.5,-40.3}],
Text[Style[" HO ",Darker[Blue],20],{33.5,-40.3}]
}
]

Porównanie metody “strzelania” z diagonalizacją “HO” dla 136Xe: funkcje falowe

Przycisk "Zapisz" w Manipulacji 3 (slajd 11) rejestrował m.in. postać funkcji falowej w zależności od r w fm i znormalizowanej. Wyznaczając stany dla 136Xe, zapisałem na dysku 3 funkcje falowe. Porównany je teraz z wynikiem drugiej metody.

Wczytujemy dane z dysku, przykładowe funkcje z metody “strzelania”:

SetDirectory["d:/Pfutzner/Teach/Mathematica4Nuclear/05_ModelPowlokowy"]

FileNames["Wave*.dat"]

n1d32=Import["Wave_136Xe_n1d32.dat","Table"](* (Subscript[ν1d, 3/2]) *);
p0f72=Import["Wave_136Xe_p0f72.dat","Table"](* (Subscript[π0f, 7/2]) *);
p2s12=Import["Wave_136Xe_p2s12.dat","Table"](* (Subscript[π2s, 1/2]) *);

Wykresy dla stanu ν 1(Subscript[d, 3/2])

gSHOOTn1d32=ListPlot[n1d32,
PlotStyle->{Red,PointSize[0.008]},
ImageSize->500]

Jeszcze raz stosujemy metodę HO, tym razem zwracając uwagę na współczynniki rozkładu.

sol=matrixHO[82,54,-1,3/2,2]

Stan 1(Subscript[d, 3/2]) jest na 5-tej pozycji. Współczynniki rozkładu dla niego:

coefn1d32=sol[[2,5]]

Obliczamy funkcje składowe (aby je narysować) oraz ich sumę, która reprezentuje funkcję falową. Pamiętamy, że w metodzie “strzelania” wyznaczalismy (Subscript[u, nlj])(r)=r (Subscript[R, nlj])(r), podczas gdy metoda HO daje  (Subscript[R, nlj])(r). Mnożymy więc rozwiązania metody HO przez r. Okazuje się, że trzeba jeszcze dodać czynnik -1 (konwencja fazy).

tabHOn1d32=Table[{-gnl[136,n,2,r]*r*coefn1d32[[n+1]]},{n,0,5}];
solHOn1d32=Sum[-gnl[136,n,2,r]*coefn1d32[[n+1]]*r,{n,0,5}];

gHOn1d32Sum=Plot[solHOn1d32,{r,0,20},
Frame->True,
PlotRange->{{0,14},{-0.55,.55}},
FrameLabel->{"r [fm]","r·(Subscript[R, njl])(r)"},
PlotStyle->Directive[Black,Thickness[0.004]],
BaseStyle->{Black,FontFamily->"Verdana",FontSize->20},
GridLines->Automatic,
GridLinesStyle->Dashed,
ImageSize->700,
Epilog->{Text[Style["(^136)Xe",Black,30],{9,0.3}]}
]

gHOn1d32=Plot[tabHOn1d32,{r,0,20},
Frame->True,
PlotRange->{{0,14},{-0.55,.55}},
FrameLabel->{"r [fm]","r·(Subscript[R, njl])(r)"},
PlotStyle->{Directive[Magenta,Thickness[0.003]],
Directive[Blue, Thickness[0.003]],
Directive[Orange, Thickness[0.003]],
Directive[Green,Thickness[0.003]],
Directive[Brown,Thickness[0.003]],
Directive[Pink,Thickness[0.003]]},
BaseStyle->{Black,FontFamily->"Verdana",FontSize->20},
GridLines->Automatic,
GridLinesStyle->Dashed,
ImageSize->700,
Epilog->{Text[Style["(^136)Xe",Black,30],{9,0.3}]}
]

Show[gHOn1d32Sum,gHOn1d32,gSHOOTn1d32]

To samo dla π 0(Subscript[f, 7/2])

gSHOOTp0f72=ListPlot[p0f72,
PlotStyle->{Red,PointSize[0.008]},
ImageSize->500]

sol=matrixHO[82,54,1,7/2,3]

Stan 0(Subscript[f, 7/2]) jest na 4-tej pozycji. Współczynniki rozkładu dla niego:

coefp0f72=sol[[2,4]]

tabHOp0f72=Table[{gnl[136,n,3,r]*r*coefp0f72[[n+1]]},{n,0,5}];
solHOp0f72=Sum[gnl[136,n,3,r]*coefp0f72[[n+1]]*r,{n,0,5}];

gHOp0f72Sum=Plot[solHOp0f72,{r,0,20},
Frame->True,
PlotRange->{{0,14},{-0.1,.7}},
FrameLabel->{"r [fm]","r·(Subscript[R, njl])(r)"},
PlotStyle->Directive[Black,Thickness[0.004]],
BaseStyle->{Black,FontFamily->"Verdana",FontSize->20},
GridLines->Automatic,
GridLinesStyle->Dashed,
ImageSize->700,
Epilog->{Text[Style["(^136)Xe",Black,30],{9,0.54}]}
]

gHOp0f72=Plot[tabHOp0f72,{r,0,20},
Frame->True,
PlotRange->{{0,14},{-0.1,.7}},
FrameLabel->{"r [fm]","r·(Subscript[R, njl])(r)"},
PlotStyle->{Directive[Magenta,Thickness[0.003]],
Directive[Blue, Thickness[0.003]],
Directive[Orange, Thickness[0.003]],
Directive[Green,Thickness[0.003]],
Directive[Brown,Thickness[0.003]],
Directive[Pink,Thickness[0.003]]},
BaseStyle->{Black,FontFamily->"Verdana",FontSize->20},
GridLines->Automatic,
GridLinesStyle->Dashed,
ImageSize->700,
Epilog->{Text[Style["(^136)Xe",Black,30],{9,0.54}]}
]

Show[gHOp0f72Sum,gHOp0f72,gSHOOTp0f72]

... i dla π 2(Subscript[s, 1/2])

gSHOOTp2s12=ListPlot[p2s12,
PlotStyle->{Red,PointSize[0.008]},
ImageSize->500]

sol=matrixHO[82,54,1,1/2,0]

Stan 2(Subscript[s, 1/2]) jest na 6-tej pozycji. Współczynniki rozkładu dla niego:

coefp2s12=sol[[2,6]]

tabHOp2s12=Table[{gnl[136,n,0,r]*r*coefp2s12[[n+1]]},{n,0,5}];
solHOp2s12=Sum[gnl[136,n,0,r]*coefp2s12[[n+1]]*r,{n,0,5}];

gHOp2s12Sum=Plot[solHOp2s12,{r,0,20},
Frame->True,
PlotRange->{{0,14},{-0.5,.6}},
FrameLabel->{"r [fm]","r·(Subscript[R, njl])(r)"},
PlotStyle->Directive[Black,Thickness[0.004]],
BaseStyle->{Black,FontFamily->"Verdana",FontSize->20},
GridLines->Automatic,
GridLinesStyle->Dashed,
ImageSize->700,
Epilog->{Text[Style["(^136)Xe",Black,30],{11,0.3}]}
]

gHOp2s12=Plot[tabHOp2s12,{r,0,20},
Frame->True,
PlotRange->{{0,14},{-0.5,.6}},
FrameLabel->{"r [fm]","r·(Subscript[R, njl])(r)"},
PlotStyle->{Directive[Magenta,Thickness[0.003]],
Directive[Blue, Thickness[0.003]],
Directive[Orange, Thickness[0.003]],
Directive[Green,Thickness[0.003]],
Directive[Brown,Thickness[0.003]],
Directive[Pink,Thickness[0.003]]},
BaseStyle->{Black,FontFamily->"Verdana",FontSize->20},
GridLines->Automatic,
GridLinesStyle->Dashed,
ImageSize->700,
Epilog->{Text[Style["(^136)Xe",Black,30],{11,0.3}]}
]

Show[gHOp2s12Sum,gHOp2s12,gSHOOTp2s12]

Zestawienie funkcji falowych

Manipulate[
TabView[{
{1,Style[" ν 1(Subscript[d, 3/2]) ",Black,22]->If[gselect,Show[gHOn1d32Sum,gSHOOTn1d32],Show[gHOn1d32Sum,gHOn1d32]]
},
{2,Style[" π 0(Subscript[f, 7/2]) ",Black,22]->If[gselect,Show[gHOp0f72Sum,gSHOOTp0f72],Show[gHOp0f72Sum,gHOp0f72]]
},
{3,Style[" π 2(Subscript[s, 1/2]) ",Black,22]->If[gselect,Show[gHOp2s12Sum,gSHOOTp2s12],Show[gHOp2s12Sum,gHOp2s12]]
}
},Dynamic[t]],
Delimiter,
{{gselect,True ,Style[" Wyświetl też ",20]},{True->Style[" wyniki Shoting ",Red,20],False->Style[" składowe HO ",Blue,20]},ControlType->SetterBar,Background->Gray,LabelStyle->Directive[Magenta,18]},
{{t,1},None},
SaveDefinitions->True,
FrameLabel->{None,None, Style["Radialna postać funkcji falowej z metody oscylatora harmonicznego",Blue,25], None}
]