| References
[1] T. Nagao and J. Igarashi, Phys. Rev. B 74, 104404 (2006); 82 024402 (2010); J. Phys. Soc. Jpn. 77 084710 (2008).
[2] S. W. Lovesey and S. P. Collins: ”X-Ray Scattering and Absorption by Magnetic Materials”, (Oxford).
共鳴X線（弾性）散乱の散乱強度（散乱断面積），偏光依存性，および散乱後の偏光状態は，σσ’, σπ’, πσ’, ππ’の4つの散乱過程についての構造因子（散乱振幅）が求まれば，それを2x2行列として表し，Stokesパラメータ(P1,P2,P3)を使った計算方法に当てはめることで，すべて計算することができる．ここに，Loveseyの教科書で与えられている式を用いる．また，構造因子の計算には，Nagao-Igarashiによって導かれた，多極子のrank別に記述された幾何学的構造因子を用いる．
Euler Rotationの定義 (* Euler rotation matrix *)
rotalp[alp_] := {{Cos[alp], -Sin[alp], 0}, {Sin[alp], Cos[alp], 0}, {0, 0, 1}}
rotbet[bet_] := {{Cos[bet], 0, Sin[bet]}, {0, 1, 0}, {-Sin[bet], 0, Cos[bet]}}
rotgam[gam_] := {{Cos[gam], -Sin[gam], 0}, {Sin[gam], Cos[gam], 0}, {0, 0, 1}}
rotEuler[alp_, bet_, gam_] := rotalp[alp].rotbet[bet].rotgam[gam];
X線の波数ベクトルと偏光ベクトルの定義 Bragg角thetaで入射(inc)および反射(out)するX線の波数ベクトルと偏光ベクトルを次のように定義する．実験室座標系での表記．偏光ベクトルはsigmaとpiを2要素として定義．いずれも大きさ1の単位ベクトル．ε_σ x ε_π = k の関係（ε_σからε_πへ右ねじを回したときにねじが進む向きが波数ベクトルk）になるようにする (ε_π x ε_σ = k にすると，F_σπとF_πσの符号が逆転し，円偏光を扱うときや，F_πσとF_ππの両方，あるいはF_σσとF_σπの両方が関与するときなどに結果が違ってくる）．
(* Wave vectors and polarization vectors *)
Clear[theta];
kinc0[theta_]: = {0, Cos[theta], -Sin[theta]};
kout0 [theta_]: = {0, Cos[theta], Sin[theta]};
epsinc0 [theta_]: = {{1, 0, 0}, {0, -Sin[theta], -Cos[theta]}};
epsout0 [theta_]: = {{1, 0, 0}, {0, Sin[theta], -Cos[theta]}};
計算モジュール xyz座標系立方テンソル表記でのX線テンソルを計算するためのqtE2[a_, b_]と計算モジュールcalcXrayVecsCubeを定義しておく．calcXrayVecsCubeを実行すると，入射X線の波数ベクトルkincと偏光ベクトルepsinc，散乱X線の波数ベクトルkoutと偏光ベクトルepsoutに対して，幾何学的散乱因子pfacZ0E1, pfacZ1E1, pfacZ2E1, pfacZ0E2, pfacZ1E2, pfacZ2E2, pfacZ3E2, pfacZ4E2が計算される．
(* N.I's X-ray E2 rank-2 tensor *)
qtE2[a_, b_] := {(3*a[[3]]*b[[3]] - a.b)/2,
(a[[1]]*b[[1]] - a[[2]]*b[[2]])*Sqrt[3]/2,
(a[[2]]*b[[3]] + a[[3]]*b[[2]])*Sqrt[3]/2,
(a[[3]]*b[[1]] + a[[1]]*b[[3]])*Sqrt[3]/2,
(a[[1]]*b[[2]] + a[[2]]*b[[1]])*Sqrt[3]/2}
(* Levi-Civita tensor *)
lceps = {{{0, 0, 0}, {0, 0, 1}, {0, -1, 0}}, {{0, 0, -1}, {0, 0, 0},
{1, 0, 0}}, {{0, 1, 0}, {-1, 0, 0}, {0, 0, 0}}};
(* Calculation Module , 2017.5.24 *)
calcXrayVecsCube := Module[{},
pfacZ0E1 = -(eout.einc)/Sqrt[3];
pfacZ1E1 = I*Cross[eout, einc]/Sqrt[2];
pfacZ2E1 = qtE2[eout, einc]*Sqrt[2/3];
pfacZ0E2 = ((kout.kinc)*(eout.einc) + (kout.einc)*(eout.kinc))/2/Sqrt[5];
pfacZ1E2 = ((eout.einc)*Cross[kout, kinc] + (kout.kinc)*Cross[eout, einc]
+ (kout.einc)*Cross[eout, kinc] + (eout.kinc)*Cross[kout, einc])*I/Sqrt[10]/2;
pfacZ2E2 = Chop[((eout.einc)*qtE2[kinc, kout] + (kout.kinc)*qtE2[einc, eout]
+ qtE2[Cross[kout, kinc], Cross[eout, einc]])*(-1/Sqrt[14])];
pfacZ3E2 = Flatten[Chop[{(Sum[Cross[kout, kinc][[i]]*qtE2[eout, einc][[i + 2]]
+ Cross[eout, einc][[i]]*qtE2[kout, kinc][[i + 2]] + Cross[kout, einc][[i]]*qtE2[eout, kinc][[i + 2]]
+ Cross[eout, kinc][[i]]*qtE2[kout, einc][[i + 2]], {i, 1, 3}])*(-I/Sqrt[2]/6),
Table[(Cross[kout, kinc][[i]]*eout[[i]]*einc[[i]] + Cross[eout, einc][[i]]*kout[[i]]*kinc[[i]]
+ Cross[kout, einc][[i]]*eout[[i]]*kinc[[i]] + Cross[eout, kinc][[i]]*kout[[i]]*einc[[i]]
)*(-I)/4*Sqrt[5/2] + 1/2*pfacZ1E2[[i]], {i, 1, 3}],
Table[(Cross[kout, kinc][[i]]*Sum[lceps[[i, j, k]]*(eout[[j]]*einc[[j]] - eout[[k]]*einc[[k]]), {j, 1, 3}, {k, 1, 3}]
+ Cross[eout, einc][[i]]*Sum[lceps[[i, j, k]]*(kout[[j]]*kinc[[j]] - kout[[k]]*kinc[[k]]), {j, 1, 3}, {k, 1, 3}]
+ Cross[kout, einc][[i]]*Sum[lceps[[i, j, k]]*(eout[[j]]*kinc[[j]] - eout[[k]]*kinc[[k]]), {j, 1, 3}, {k, 1, 3}]
+ Cross[eout, kinc][[i]]*Sum[lceps[[i, j, k]]*(kout[[j]]*einc[[j]] - kout[[k]]*einc[[k]]), {j, 1, 3}, {k, 1, 3}]
)*(-I)/8*Sqrt[3/2], {i, 1, 3}]}]];
pfacZ4E2 = Flatten[Chop[{
(5*(kout[[1]]*kinc[[1]]*eout[[1]]*einc[[1]] + kout[[2]]*kinc[[2]]*eout[[2]]*einc[[2]]
+ kout[[3]]*kinc[[3]]*eout[[3]]*einc[[3]])/2 - Sqrt[5]*pfacZ0E2)*Sqrt[2/15],
(-Sqrt[14]*(kout[[1]]*kinc[[1]]*eout[[1]]*einc[[1]]
+ kout[[2]]*kinc[[2]]*eout[[2]]*einc[[2]] - 2*kout[[3]]*kinc[[3]]*eout[[3]]*einc[[3]])
- 2*Sqrt[2/7]*((kout.kinc)*qtE2[eout, einc][[1]]+ (eout.einc)*qtE2[kout, kinc][[1]]
+ (kout.einc)*qtE2[eout, kinc][[1]]+ (eout.kinc)*qtE2[kout, einc][[1]]))/Sqrt[3]/2,
(Sqrt[14]*(kout[[1]]*kinc[[1]]*eout[[1]]*einc[[1]]- kout[[2]]*kinc[[2]]*eout[[2]]*einc[[2]])
- 2*Sqrt[2/21]*((kout.kinc)*qtE2[eout, einc][[2]] + (eout.einc)*qtE2[kout, kinc][[2]]
+ (kout.einc)*qtE2[eout, kinc][[2]] + (eout.kinc)*qtE2[kout, einc][[2]]))/2,
Table[(qtE2[kout, kinc][[i + 2]]*Sum[lceps[[i, j, k]]*(eout[[j]]*einc[[j]] - eout[[k]]*einc[[k]]), {j, 1, 3}, {k, 1, 3}]
+ qtE2[eout, einc][[i + 2]]*Sum[lceps[[i, j, k]]*(kout[[j]]*kinc[[j]] - kout[[k]]*kinc[[k]]), {j, 1, 3}, {k, 1, 3}]
+ qtE2[kout, einc][[i + 2]]*Sum[lceps[[i, j, k]]*(eout[[j]]*kinc[[j]] - eout[[k]]*kinc[[k]]), {j, 1, 3}, {k, 1, 3}]
+ qtE2[eout, kinc][[i + 2]]*Sum[lceps[[i, j, k]]*(kout[[j]]*einc[[j]] - kout[[k]]*einc[[k]]), {j, 1, 3}, {k, 1, 3}]
)/4/Sqrt[6], {i, 1, 3}],
Table[((7*kout[[i]]*kinc[[i]] - 3*(kout.kinc))*qtE2[eout, einc][[i + 2]]
+ (7*eout[[i]]*einc[[i]] - 3*(eout.einc))*qtE2[kout, kinc][[i + 2]]
+ (7*eout[[i]]*kinc[[i]] - 3*(eout.kinc))*qtE2[kout, einc][[i + 2]]
+ (7*kout[[i]]*einc[[i]] - 3*(kout.einc))*qtE2[eout, kinc][[i + 2]]
)/Sqrt[42]/2, {i, 1, 3}]}]];
]
散乱条件 X線のエネルギー，回折指数，two-theta，thetaの計算．これで，kinc0, kout0, epsinc0,epsout0が決まる．
ただし，結晶の逆格子基本ベクトルb1, b2, b3は別に定義しておく．結晶の構造因子も計算しておくとよい．
(* X-ray energy, wavenumber, incident k-vector *)
engy = 6.712;
lambda = 12.398/engy;
kxray = 2*Pi/lambda;
(* reflection index *)
{h, k, l} = {0, 0.75, 7};
kappa = b1*h + b2*k + b3*l;
tau = Sqrt[Sum[kappa[[i]]^2, {i, 1, 3}]];
tth = ArcSin[tau/2/kxray]*180/Pi*2; theta = ArcSin[tau/2/kxray];
Print["kappa=", kappa]
Print["tau=", tau, Tab, "tau/4Pi=", tau/4/Pi]
Print["two-theta=", tth]
着目原子に対するEuler Rotationの設定 （例）立方晶CeB6の(3/2, 3/2, 1/2)反射で[-1 1 0]//実験室系X軸，[-3 -3 -1]//実験室系Z軸，結晶軸abcをxyz軸にとるとき，次のようにする．うまく設定できていれば，testの結果がすべてゼロになるはず．
(* XYZ-vectors e1, e2, e3, expressed by the ionic basis xyz. Q // [3, 3, 1] *)
e10={-1, 1, 0}/Sqrt[2]; e30={-3, -3, -1}/Sqrt[19]; e20=Cross[e30, e10];
alpha=-Pi/2; beta=ArcCos[-1/Sqrt[19]]; gamma=-Pi/4;
Print["e10=", e10, ", e20=", e20, ", e30=", e30]
(* Euler angle test *)
Do[euangs = {alpha + psi*Pi/180, beta, gamma};
chi = psi*Pi/180;
{e1,e2,e3}={e10*Cos[chi]-e20*Sin[chi], e10*Sin[chi]+e20*Cos[chi], e30};
Print[psi, Tab,
N[{1, 0, 0}.rotEuler[euangs[[1]], euangs[[2]], euangs[[3]]] - e1],
N[{0, 1, 0}.rotEuler[euangs[[1]], euangs[[2]], euangs[[3]]] - e2], N[{0, 0, 1}.rotEuler[euangs[[1]], euangs[[2]], euangs[[3]]] - e3]]
,{psi, 0, 180, 30}];
秩序変数の構造因子 目的の反射に対する各ランク秩序変数の構造因子stfacdip, stfacquad, stfacoct, stfachexaを入力．
磁気双極子は{Jx,Jy,Jz}の3成分，電気四極子は{O20,O22,Oyz,Ozx,Oxy}の5成分，磁気八極子は{Txyz,Tax,Tay,Taz,Tbx,Tby,Tbz}の7成分，電気十六極子は{Ha,Heg1,Heg2,H1gx,H1gy,H1gz,H2gx,H2gy,H2gz}の9成分．
例：
stfacdip={1,0,0}
stfacquad={1,0,0,0,1}
stfacoct={1,0,0,0,1,0,0}
stfachexa={0,1,0,0,0,1,0,0,1}
σσ’, σπ’, πσ’, ππ’の4つの散乱過程についての構造因子の計算 E1共鳴のrank-1,2，E2共鳴のrank-1,2,3,4について，σσ’, σπ’, πσ’, ππ’の4つの散乱過程についての構造因子（散乱振幅）を計算する．計算は原子のxyz座標系について行う．
(* azimuthal angle *)
psi = 0;
(* Multipole scattering amplitudes for s-s', s-p', p-s', and p-p' *)
polif = {{0, 0}, {0, 90}, {90, 0}, {90, 90}};
polifcalc4 = Table[thpolinc = N[polif[[ii, 1]]*Pi/180];
thpolout = N[polif[[ii, 2]]*Pi/180];
epsinc = N[epsinc0[theta][[1]]*Cos[thpolinc] + epsinc0[theta][[2]]*Sin[thpolinc]];
epsout = N[epsout0[theta][[1]]*Cos[thpolout] + epsout0[theta][[2]]*Sin[thpolout]];
euangs = N[{alpha + psi*Pi/180, beta, gamma}];
kinc = kinc0[theta].rotEuler[euangs[[1]], euangs[[2]], euangs[[3]]];
kout = kout0[theta].rotEuler[euangs[[1]], euangs[[2]], euangs[[3]]];
einc = Chop[epsinc.rotEuler[euangs[[1]], euangs[[2]], euangs[[3]]]];
eout = Chop[epsout.rotEuler[euangs[[1]], euangs[[2]], euangs[[3]]]];
calcXrayVecsCube;
z1e1 = Chop[pfacZ1E1.stfacdip];
z2e1 = Chop[pfacZ2E1.stfacquad];
z1e2 = Chop[pfacZ1E2.stfacdip];
z2e2 = Chop[pfacZ2E2.stfacquad];
z3e2 = Chop[pfacZ3E2.stfacoct];
z4e2 = Chop[pfacZ4E2.stfachexa];
{z1e1, z2e1, z1e2, z2e2, z3e2, z4e2}, {ii, 1, 4}];
{z1e1, z2e1, z1e2, z2e2, z3e2, z4e2}の各要素にかかる係数 6種類の構造因子{z1e1, z2e1, z1e2, z2e2, z3e2, z4e2}それぞれにかかる係数を入力する．
エネルギー依存性f(E)=f'(E)+if''(E)を表す部分と考えてよい．全体の散乱振幅Gは
G=f11(E)*z1e1+f21(E)*z2e1+f12(E)*z1e2+f22(E)*z2e2+f32(E)*z3e2+f42(E)*z4e2
と書ける．エネルギー依存性はそれぞれの項で異なっていると考える．
(* factor for each amplitudes *)
facamp = {1, 1, 0, 1, 1, 0};
非共鳴Thomson散乱 非共鳴Thomson散乱があるときは，これを足す．
(*Non-resonant Thomson*)
fc=0;
gNRT = fc*{{1, 0}, {0, Cos[tth*Pi/180]}};
散乱振幅行列の表示 σσ’, σπ’, πσ’, ππ’の4つの散乱過程の散乱振幅を行列要素とする，2x2の散乱振幅行列Gを計算する．これが，散乱に関するすべての情報を含んでいる．
(*Scattering amplitude matrix*)
gmat = gNRT + {{
Sum[facamp[[i]]*polifcalc4[[1, i]], {i, 1, 6}],
Sum[facamp[[i]]*polifcalc4[[3, i]], {i, 1, 6}]},
{Sum[facamp[[i]]*polifcalc4[[2, i]], {i, 1, 6}],
Sum[facamp[[i]]*polifcalc4[[4, i]], {i, 1, 6}]}};
gbet = ComplexExpand[(gmat[[1, 1]] + gmat[[2, 2]])/2];
galp3 = ComplexExpand[(gmat[[1, 1]] - gmat[[2, 2]])/2];
galp1 = ComplexExpand[(gmat[[1, 2]] + gmat[[2, 1]])/2];
galp2 = ComplexExpand[I*(gmat[[1, 2]] - gmat[[2, 1]])/2];
galp = {galp1, galp2, galp3};
Print["gbet=",Chop[gbet],", galp=",Chop[galp]]
Print["Gmat=",
MatrixForm[
Chop[{{gbet + galp3, galp1 - I*galp2}, {galp1 + I*galp2, gbet - galp3}}]]]
アナライザー情報 (* Analyzer d-spacing *)
dtha=1.43167; (*Al_220: 1.43167, Mo_200: 1.5734, Cu_110: 1.27798*)
tha=ArcSin[2*Pi/dtha/2/kxray]; tthA=2*tha;
Print["thA=", tha*180/Pi, ", tthA=", tthA*180/Pi]
Print["cos(ttha)^2=", Cos[tthA]^2]
散乱強度，散乱後の偏光状態，アナライザー後の検出器での観測強度 入射X線の偏光状態をStokesパラメータpoli=(P1,P2,P3)で定義．以下の例は直線偏光（偏光角eta）を仮定．
phiAnaを入れると，アナライザーの後の検出器での強度intを計算．このルーチンが基本．
(*eta: incident linear polarization angle*)
p1 = Sin[2*eta*Pi/180]; p2 = 0; p3 = Cos[2*eta*Pi/180];
poli = {p1, p2, p3};
(*scattering cross-section*)
cross = Chop[ Conjugate[gbet]*gbet+Conjugate[galp].galp
+Conjugate[gbet]*(poli.galp)+(poli.Conjugate[galp])*gbet
+I*poli.Cross[Conjugate[galp],galp]];
(*Polarization of the scattered x-ray*)
crosspolf = Chop[Conjugate[gbet]*galp+Conjugate[galp]*gbet
-I*Cross[Conjugate[galp],galp]+(Conjugate[gbet]*gbet)*poli
-I*Conjugate[gbet]*Cross[poli,galp]+I*Cross[poli,Conjugate[galp]]*gbet
+Conjugate[galp]*(poli.galp)-Cross[Conjugate[galp],Cross[poli, galp]]];
polf = crosspolf/cross;
(*Intensity at the detector after analyzer*)
p3A = polf[[1]]*Sin[2*phiAna*Pi/180]+polf[[3]]*Cos[2*phiAna*Pi/180];
int = cross*(1-(1-p3A)/2*Sin[tthA]^2)
例1: phiA依存性 入射X線が直線偏光のとき，偏光角eta=0,30,60,90に対して，強度のphiA (POL) 依存性を計算する．
(*incident polarization angles, eta *)
etalist = {0, 30, 60, 90};
(*phiA dependence of the intensity*)
Do[eta = N[etalist[[ii]]];
p1 = Sin[2*eta*Pi/180]; p2 = 0; p3 = Cos[2*eta*Pi/180];
poli = {p1, p2, p3};
cross = Chop[ Conjugate[gbet]*gbet+Conjugate[galp].galp
+Conjugate[gbet]*(poli.galp)+(poli.Conjugate[galp])*gbet
+I*poli.Cross[Conjugate[galp],galp]];
crosspolf = Chop[Conjugate[gbet]*galp+Conjugate[galp]*gbet
-I*Cross[Conjugate[galp],galp]+(Conjugate[gbet]*gbet)*poli
-I*Conjugate[gbet]*Cross[poli,galp]+I*Cross[poli,Conjugate[galp]]*gbet
+Conjugate[galp]*(poli.galp)-Cross[Conjugate[galp],Cross[poli, galp]]];
polf = crosspolf/cross;
phiAnadep[etalist[[ii]]] = Table[{phiAna,
p3A = polf[[1]]*Sin[2*phiAna*Pi/180]+polf[[3]]*Cos[2*phiAna*Pi/180];
int = cross*(1-(1-p3A)/2*Sin[tthA]^2)}, {phiAna, -120, 120, 5}]
, {ii, 1, Length[etalist]}]
ListPlot[{phiAnadep[0], phiAnadep[30], phiAnadep[60], phiAnadep[90]},
PlotStyle -> PointSize[Medium], BaseStyle -> {Medium},
PlotRange -> {All, {0, All}}, AxesLabel -> {"phiA", ""}]
例2: eta依存性 入射X線が直線偏光のとき，散乱強度，散乱後の偏光状態の入射X線偏光角eta依存性を計算する．さらに，その結果を用いて，phiA (POL) =0,30,60,90のときの強度のeta依存性を計算する．
(*eta dependence of cross-section and polf*)
polidep = Table[
p1=Sin[2*eta*Pi/180]; p2=0; p3=Cos[2*eta*Pi/180];
poli = {p1, p2, p3};
cross = Chop[ Conjugate[gbet]*gbet+Conjugate[galp].galp +Conjugate[gbet]*(poli.galp)+(poli.Conjugate[galp])*gbet
+I*poli.Cross[Conjugate[galp],galp]];
crosspolf = Chop[Conjugate[gbet]*galp+Conjugate[galp]*gbet
-I*Cross[Conjugate[galp],galp]+(Conjugate[gbet]*gbet)*poli
-I*Conjugate[gbet]*Cross[poli, galp]+I*Cross[poli,Conjugate[galp]]*gbet +Conjugate[galp]*(poli.galp)-Cross[Conjugate[galp],Cross[poli, galp]]];
polf = crosspolf/cross;
{eta, cross, polf, Sum[polf[[i]]^2, {i, 1, 3}]}
, {eta, -100, 120, 5}];
(*Extract data and draw graph*)
policross=Table[{polidep[[i,1]],polidep[[i,2]]},{i,1,Length[polidep]}];
polidepp1=Table[{polidep[[i,1]],polidep[[i,3,1]]},{i,1,Length[polidep]}];
polidepp2=Table[{polidep[[i,1]],polidep[[i,3,2]]},{i,1,Length[polidep]}];
polidepp3=Table[{polidep[[i,1]],polidep[[i,3,3]]},{i,1,Length[polidep]}];
GraphicsArray[{ListPlot[policross,AxesLabel->{"eta",""},
PlotStyle->PointSize[Medium],BaseStyle->{Medium}],
ListPlot[{polidepp1,polidepp2,polidepp3},AxesLabel->{"eta",""},
PlotStyle->{PointSize[Medium],Directive[Green,PointSize[Medium]],
Directive[Red,PointSize[Medium]]},BaseStyle->{Medium}]}]
(*analyzer pol angles, phiA*)
phiAnalist={0,30,60,90};
(*eta dependence of intensity after POLana*)
Do[
phiAna=phiAnalist[[ii]]*Pi/180;
polidepphiA[phiAnalist[[ii]]]=Table[{polidep[[i,1]],
cross=polidep[[i,2]]; polf=polidep[[i,3]];
p3A=polf[[1]]*Sin[2*phiAna]+polf[[3]]*Cos[2*phiAna];
int=cross*(1-(1-p3A)/2*Sin[tthA]^2)},{i,1,Length[polidep]}];
,{ii,1,Length[phiAnalist]}]
(*Extract data and draw graph*)
ListPlot[{polidepphiA[0],polidepphiA[30],polidepphiA[60],polidepphiA[90]},
PlotStyle->{PointSize[Medium],Directive[Green,PointSize[Medium]],
Directive[Red,PointSize[Medium]],Directive[Gray,PointSize[Medium]]},
BaseStyle->{Medium},PlotRange->{All,{0,All}}]
|
