原子物理学: 遷移強度比計算機作った
pythonにqutipっていう便利なライブラリがあって、それでクレプシュ・ゴルダン係数が計算できるから3jとか6jとか難しいものを使わなくても単純な基底の変換と愚直な積分で光学遷移強度比とかを計算できそうなので作ってみた。
(超)微細構造までクレプシュ・ゴルダン分解して積分の形にするのは教科書とかに書いてあるようにやる。簡単に下に要約する。
$${\bm{r}}$$を位置ベクトル、$${\bm{p}}$$を運動量ベクトル、$${\bm{L}=\bm{r} \times \bm{p}}$$を軌道角運動量ベクトル、$${\bm{S}}$$をスピン、$${\bm{J}=\bm{L}+\bm{S}}$$を全軌道角運動量、$${M_J}$$を磁気副準位とする。原子の固有ケットは$${|\alpha ,J,M_J>}$$と表すことができる。双極子遷移行列の要素は、
$$
d_{eg}=< \alpha ',J'M'_J | \hat{\epsilon} \cdot e \bm{r}
| \alpha ,J,M_J >
$$
とあらわされる。ここで、$${e}$$は電荷素量、$${\hat{\epsilon}}$$は光の偏光ベクトル、$${\alpha ',J',M'J}$$は下準位の固有ベクトルの固有値を表す。固有ケット$${|\alpha ,J,M_J>}$$は$${L}$$と$${S}$$の波動関数とClebsh-Gordan係数$${C(L,M_L,S,M_S,J,M_J)}$$を用いて表すことができる。
$$
|\alpha ,J,M_J>=\sum_{M_L,M_S}C(L,M_L,S,M_S,J,M_J) |\alpha ,L,M_L> |S,M_S>
$$
$${|\alpha ,L,M_L>|S,M_S>\equiv |L,M_L:S,M_S>|\alpha>}$$と定義する。
これを上式に代入して、
$$
d_{eg}\propto\sum_{M_L,',M_S',M_L,M_S}C(L',M_L',S',M_S',J',M_J')C(L,M_L,S,M_S,J,M_J) \\
\times < L', L_M':S', M_S' | \epsilon _q | L, M_L:S, M_S > <\alpha '|r |\alpha>
$$
ここで、$${q=0,\pm1}$$であり、それぞれ$${\pi}$$、$${\sigma_{\pm}}$$偏光に対応する。
この前半部分のブラケットは球座標上で以下のように表現され、積分を実行できる。
$$
< L', L_M':S', M_S' | \epsilon _q | L, M_L:S, M_S > \\
= \delta (S', S) \delta (M_S', M_S) \int_{0}^{\pi} d \theta \int_{0}^{2 \pi} d \phi \\
\times Y^*_{L', M_L'} (\theta, \phi) Y _{L, M_L} (\theta, \phi) \sin(\theta) \cos(\theta-q\pi/2) \exp(qi\phi) \cos(-q\pi/4)
$$
ここで$${Y_{L,M_L}(\theta,\phi)}$$は球面調和関数、$${\delta}$$はクロネッカーのデルタである。
この計算を実行することで同じ$${S}$$、$${L}$$を持つ微細準位間の相対的な遷移強度を計算することができる。
同様の計算を超微細構造に対して行う(角運動量のクレプシュ・ゴルダン分解をさらにもう一回する)ことで同じ$${S}$$、$${L}$$、$${J}$$を持つ超微細準位間の相対的な遷移強度を計算することもできる。
ということで、この計算をpythonでコーディングする。
import numpy as np
from math import sqrt
from qutip import *
import numpy.linalg as LA
#import matplotlib.pyplot as plt
import scipy
import numpy
from scipy.integrate import quad,dblquad
from sympy.physics.quantum.cg import CG
import time
import tkinter
#right clic#
def make_menu(w):
global the_menu
the_menu = tkinter.Menu(w, tearoff=0)
the_menu.add_command(label="cut")
the_menu.add_command(label="copy")
the_menu.add_command(label="paste")
def show_menu(e):
w = e.widget
the_menu.entryconfigure("cut", command=lambda: w.event_generate("<<Cut>>"))
the_menu.entryconfigure("copy", command=lambda: w.event_generate("<<Copy>>"))
the_menu.entryconfigure("paste", command=lambda: w.event_generate("<<Paste>>"))
the_menu.tk.call("tk_popup", the_menu, e.x_root, e.y_root)
def complex_quadrature(func, a, b, c, d, **kwargs):
def real_func(x,y):
return numpy.real(func(x,y))
def imag_func(x,y):
return numpy.imag(func(x,y))
real_integral = dblquad(real_func, a, b, c, d, **kwargs)
imag_integral = dblquad(imag_func, a, b, c, d, **kwargs)
return (real_integral[0] + 1j*imag_integral[0], real_integral[1:], imag_integral[1:])
def drange(begin, end, step):
n = begin
while n+step <= end:
yield n
n += step
def krond(a,b):
if a==b:
return 1
else:
return 0
def d1(l0,ml0,l1,ml1,q):
return complex_quadrature(lambda theta,phi:np.conjugate(scipy.special.sph_harm(ml0,l0,theta,phi))*scipy.special.sph_harm(ml1,l1,theta,phi)*np.sin(phi)*np.cos(phi+q*np.pi/2)*np.exp(-q*1j*theta)*np.cos(q*np.pi/4),0,np.pi,0,2*np.pi)[0]
def d2(j0,mj0,l0,s0,j1,mj1,l1,s1):
d=0
for ms0 in drange(-s0,s0+1,1):
for ms1 in drange(-s1,s1+1,1):
for ml0 in drange(max(-l0,mj0-ms0),min(l0,mj0-ms0)+1,1):
for ml1 in drange(max(-l1,mj1-ms1),min(l1,mj1-ms1)+1,1):
for q in range(-1,1+1,1):
if ms0==ms1 and s0==s1:
d+=clebsch(s0,l0,j0,ms0,ml0,mj0)*clebsch(s1,l1,j1,ms1,ml1,mj1)*d1(l0,ml0,l1,ml1,q)
return d
def d3(f0,mf0,j0,i0,l0,s0,f1,mf1,j1,i1,l1,s1):
d=0
for mi0 in drange(-i0,i0+1,1):
for mi1 in drange(-i1,i1+1,1):
for mj0 in drange(max(-j0,mf0-mi0),min(j0,mf0-mi0)+1,1):
for mj1 in drange(max(-j1,mf1-mi1),min(j1,mf1-mi1)+1,1):
if mi0==mi1 and i0==i1:
d+=clebsch(i0,j0,f0,mi0,mj0,mf0)*clebsch(i1,j1,f1,mi1,mj1,mf1)*d2(j0,mj0,l0,s0,j1,mj1,l1,s1)
return d
def ButtonEvent(event):
ic = EditBox_nuclearspin.get()
#ic=input("nuclear spin I:")
i=float(ic)
sc = EditBox_totalspin.get()
#sc=input("total spin S:")
s=float(sc)
j0c = EditBox_totalJoflowerlevel.get()
#j0c=input("total J of lower level:")
j0=float(j0c)
l0c = EditBox_totalloflowerlevel.get()
#l0c=input("total L of lower level:")
l0=float(l0c)
j1c = EditBox_totalJofupperlevel.get()
#j1c=input("total J of upper level:")
j1=float(j1c)
l1c = EditBox_totallofupperlevel.get()
#l1c=input("total L of upper level:")
l1=float(l1c)
for f0 in drange(np.abs(j0-i),j0+i+1,1):
for f1 in drange(np.abs(j1-i),j1+i+1,1):
dd=0
end=0
start=time.time()
for mf0 in drange(-f0,f0+1,1) :
for mf1 in drange(-f1,f1+1,1) :
if np.conjugate(d3(f0,mf0,j0,i,l0,s,f1,mf1,j1,i,l1,s))*d3(f0,mf0,j0,i,l0,s,f1,mf1,j1,i,l1,s)>0.0001:
#print(mf0,"->",mf1,np.conjugate(d3(f0,mf0,j0,i,l0,s,f1,mf1,j1,i,l1,s))*d3(f0,mf0,j0,i,l0,s,f1,mf1,j1,i,l1,s))
dd+=np.conjugate(d3(f0,mf0,j0,i,l0,s,f1,mf1,j1,i,l1,s))*d3(f0,mf0,j0,i,l0,s,f1,mf1,j1,i,l1,s)
end=time.time()-start
if dd!=0:
print("total",":",f0,"->",f1,":",2.5*np.real(dd),"T=",end)
root = tkinter.Tk()
root.title("line strength")
root.geometry("400x220")
make_menu(root)
root.bind_class("Entry", "<Button-3><ButtonRelease-3>", show_menu)
Static_nuclearspin = tkinter.Label(text='nuclear spin I')
Static_nuclearspin.pack()
Static_nuclearspin.place(x=20, y=50)
Static_totalspin = tkinter.Label(text='total spin S')
Static_totalspin.pack()
Static_totalspin.place(x=20, y=70)
Static_totalJofupperlevel = tkinter.Label(text='total J of upper level')
Static_totalJofupperlevel.pack()
Static_totalJofupperlevel.place(x=20, y=90)
Static_totallofupperlevel = tkinter.Label(text='total L of upper level')
Static_totallofupperlevel.pack()
Static_totallofupperlevel.place(x=20, y=110)
Static_totalJoflowerlevel = tkinter.Label(text='total J of lower level')
Static_totalJoflowerlevel.pack()
Static_totalJoflowerlevel.place(x=20, y=130)
Static_totalloflowerlevel = tkinter.Label(text='total L of lower level')
Static_totalloflowerlevel.pack()
Static_totalloflowerlevel.place(x=20, y=150)
EditBox_nuclearspin = tkinter.Entry(width=15)
EditBox_nuclearspin.insert(tkinter.END,"4.5")
EditBox_nuclearspin.place(x=150, y=50)
EditBox_totalspin = tkinter.Entry(width=15)
EditBox_totalspin.insert(tkinter.END,"1")
EditBox_totalspin.place(x=150, y=70)
EditBox_totalJofupperlevel = tkinter.Entry(width=15)
EditBox_totalJofupperlevel.insert(tkinter.END,"3")
EditBox_totalJofupperlevel.place(x=150, y=90)
EditBox_totallofupperlevel = tkinter.Entry(width=15)
EditBox_totallofupperlevel.insert(tkinter.END,"2")
EditBox_totallofupperlevel.place(x=150, y=110)
EditBox_totalJoflowerlevel = tkinter.Entry(width=15)
EditBox_totalJoflowerlevel.insert(tkinter.END,"2")
EditBox_totalJoflowerlevel.place(x=150, y=130)
EditBox_totalloflowerlevel = tkinter.Entry(width=15)
EditBox_totalloflowerlevel.insert(tkinter.END,"1")
EditBox_totalloflowerlevel.place(x=150, y=150)
Button1 = tkinter.Button(text='excute', width=15)
Button1.bind("<Button-1>",ButtonEvent)
Button1.place(x=270, y=95)
Button2 = tkinter.Button(text='exit',command=root.quit, width=15)
Button2.place(x=270, y=140)
root.mainloop()出力結果
例として核スピン9/2の原子の、$${{}^3P_{2} \rightarrow {}^3D_{3}}$$軌道間光学遷移を計算してみた。
各超微細遷移の相対的な強度が分かる。
絶対値には意味はない。
$${\sigma _{\pm}}$$、$${\pi}$$偏光が均等に含まれる場合を計算した。
Tはその遷移の計算にかかった時間(秒)で、並列計算はしてない。
cpuはAMD Ryzen 9 5900X。
total : 2.5 -> 1.5 : 4.000000000000002 T= 53.95044946670532
total : 2.5 -> 2.5 : 3.142857142857144 T= 73.30421352386475
total : 2.5 -> 3.5 : 1.2571428571428576 T= 84.02154922485352
total : 3.5 -> 2.5 : 2.8571428571428594 T= 83.29183435440063
total : 3.5 -> 3.5 : 4.876190476190479 T= 103.7797338962555
total : 3.5 -> 4.5 : 3.4666666666666672 T= 115.88220715522766
total : 4.5 -> 3.5 : 1.8666666666666674 T= 112.02357363700867
total : 4.5 -> 4.5 : 5.515151515151516 T= 128.35754299163818
total : 4.5 -> 5.5 : 6.618181818181819 T= 135.29685020446777
total : 5.5 -> 4.5 : 1.0181818181818183 T= 131.53960728645325
total : 5.5 -> 5.5 : 5.012587412587415 T= 144.0786702632904
total : 5.5 -> 6.5 : 10.769230769230774 T= 147.45775246620178
total : 6.5 -> 5.5 : 0.3692307692307694 T= 143.09361577033997
total : 6.5 -> 6.5 : 3.2307692307692317 T= 149.975834608078
total : 6.5 -> 7.5 : 16.00000000000001 T= 153.03075170516968