|
中性子非弾性散乱 –– 結晶場励起に対する中性子散乱関数S(kappa,omega)の計算 ––
-
固有関数の確認
Print["Eigenenergies (K) and eigenvectors:"]
Do[Print[i, Tab, vals[[i]], vecs[[i]]], {i, 1, 2*J + 1}]
-
行列要素の計算
jxelem=Table[Chop[N[ Conjugate[vecs[[i]]].(eljx.vecs[[j]])], 0.00001]
,{i,1,2*J+1},{j,1,2*J+1}];
jyelem=Table[Chop[N[ Conjugate[vecs[[i]]].(eljy.vecs[[j]])], 0.00001]
,{i,1,2*J+1},{j,1,2*J+1}];
jzelem=Table[Chop[N[ Conjugate[vecs[[i]]].(eljz.vecs[[j]])], 0.00001]
,{i,1,2*J+1},{j,1,2*J+1}];
jelem={jxelem,jyelem,jzelem}; MatrixForm[jxelem]
MatrixForm[jyelem]
MatrixForm[jzelem]
-
温度，磁気形状因子，Debye-Waller因子
(* Temperature, magnetic formfactor, Debye-Waller factor *)
temp = 20; formfac = 1; dwfac = 0;
z = Sum[Exp[-vals[[i]]/temp], {i,1,2*J+1}];
-
静帯磁率 chi_{alpha,beta}(omega=0)
– Note – N_{A}g^2mu_{B}^2/k_{B}/11.6056=g^2*0.0323278をかけると(emu/mol)単位になる．
固有値の差と温度の比(Ei-Ej)/Tがlevelepsより小さければ，iとjは縮退した状態とみなす．
(* Static susceptibility chi_{alpha,beta}(omega=0) in (1/meV) unit. *)
leveleps = 0.001;
chi0stat[alp_,bet_]:=(Sum[Exp[-vals[[i]]/temp]/z*(Sum[
If[Abs[vals[[i]]-vals[[j]]]/temp < leveleps,
(jelem[[alp,i,j]]*jelem[[bet,j,i]]+jelem[[bet,i,j]]*jelem[[alp,j,i]])/2/temp,
(jelem[[alp,i,j]]*jelem[[bet,j,i]]+jelem[[bet,i,j]]*jelem[[alp,j,i]])/(vals[[j]]-vals[[i]])]
,{j,1,2*J+1}]), {i,1,2*J+1}]
-Sum[Exp[-vals[[i]]/temp]/z*jelem[[alp,i,i]],{i,1,2*J+1}]*
Sum[Exp[-vals[[i]]/temp]/z*jelem[[bet,i,i]],{i,1,2*J+1}]/temp)*11.6056
Print["chi0stat=", MatrixForm[Table[Chop[chi0stat[i,j],0.000001],{i,1,3},{j,1,3}]]]
Print["chi_zz=",chi0stat[3,3]*g^2*0.0323278,"(emu/mol)"]
-
遷移強度
(* Table of transition strength, in 1/meV unit. *)
strength1=Chop[Table[ Table[ If[ Abs[vals[[i]]-vals[[j]]]/temp<leveleps,
Exp[-vals[[i]]/temp]/z*(jelem[[alp,i,j]]*jelem[[bet,j,i]]
+jelem[[bet,i,j]]*jelem[[alp,j,i]])/2/temp,
Exp[-vals[[i]]/temp]/z*(jelem[[alp,i,j]]*jelem[[bet,j,i]]
+jelem[[bet,i,j]]*jelem[[alp,j,i]])/(vals[[j]]-vals[[i]])]*11.6056
, {i,1,2*J+1}, {j,1,2*J+1}], {alp,1,3}, {bet,1,3}]];
MatrixForm[Chop[strength1[[2, 2]], 0.00001]] (* The second Curie term due to static component *)
strength2 = Chop[Table[
-Sum[Exp[-vals[[i]]/temp]/z*jelem[[alp, i, i]], {i, 1, 2*J + 1}]*
Sum[Exp[-vals[[i]]/temp]/z*jelem[[bet, i, i]], {i, 1, 2*J + 1}]/temp*11.6056
, {alp, 1, 3}, {bet, 1, 3}]];
MatrixForm[strength2] (* Sum of the transition strengths is equal to the static susceptibility. *)
MatrixForm[Table[Sum[
strength1[[alp,bet,i,j]], {i,1,2*J+1}, {j,1,2*J+1}], {alp,1,3}, {bet,1,3}]
+strength2]
-
動帯磁率を表すスペクトル関数
– Note – 実部と虚部は．互いにKramers-Kronigの関係で結ばれている．
(* Lorentzian spectral functions for dynamical susceptibility *)
specfcnIm[w_,w0_,gam_]:=w*(gam/(gam^2+(w+w0)^2)+gam/(gam^2+(w-w0)^2))/2;
specfcnRe[w_,w0_,gam_]:=(((w+w0)*w0+gam^2)/((w+w0)^2+gam^2)
+(-(w-w0)*w0+gam^2)/((w-w0)^2+gam^2))/2;
-
１イオン結晶場励起に対する動帯磁率 chi_{alpha,beta}(omega)
(* Imaginary and real parts of the dynamic susceptibility chi_{alpha,beta}(omega). *)
chi0Im[alp_,bet_,omega_,gamma_]:=Chop[(Sum[strength1[[alp,bet,i,j]]*specfcnIm[omega,(vals[[j]]-vals[[i]])/11.6056,gamma],{i,1,2*J+1},{j,1,2*J+1}]+strength2[[alp,bet]]*specfcnIm[omega,0,gamma])];
chi0Re[alp_,bet_,omega_,gamma_]:=Chop[Sum[strength1[[alp,bet,i,j]]*specfcnRe[omega,(vals[[j]]-vals[[i]])/11.6056,gamma],{i,1,2*J+1},{j,1,2*J+1}]+strength2[[alp,bet]]*specfcnRe[omega,0,gamma]];
-
動帯磁率 chi_{alpha,beta}(omega)の図示 (* Plot chi_{alpha,beta}(omega). *)
alpha=3; beta=3;
fig1=Plot[Evaluate[{chi0Re[alpha,beta,w,0.4], chi0Im[alpha,beta,w,0.4]}], {w,-10,10},
PlotRange->All,Frame->True, FrameStyle->Thickness[0.002], BaseStyle->Large,
FrameLabel->{"E (meV)","chi(E) (1/meV)"}, PlotStyle->{{Thick,Dashed},{Thick}},
AspectRatio->0.8, ImageSize->400, PlotStyle->{RGBColor[1,0,0], RGBColor[0,0,1]}]
Clear[alpha,beta];
 (Fig. 1)
-
１イオン結晶場励起に対する中性子散乱関数 S(kappa, omega) (* Neutron scattering function S(kappa,omega). *) kappa = {0, 0, 1}; kap = kappa/Sqrt[Sum[kappa[[i]]^2, {i, 1, 3}]]; sqomg0[omega_, gamma_] := ((g/2*formfac)^2*Exp[-2*dwfac]/(1-Exp[-omega*11.6056/temp]) *Sum[(KroneckerDelta[alp, bet]-kap[[alp]]*kap[[bet]])*chi0Im[alp, bet, omega, gamma] , {alp, 1, 3}, {bet, 1, 3}]) fig2 = Plot[Evaluate[sqomg0[w, 0.4]], {w, -5, 10}, PlotRange->All]
 (Fig. 2)
|
