ラザフォード散乱における粒子の軌道【SageMath】【Python】

小出昭一郎『量子力学（Ⅰ）』の 5.2 を参考にラザフォード散乱における粒子の軌道を可視化します。用いる変数を次のように定義します。

$r : アルファ粒子と原子核の距離\\ \varphi : アルファ粒子と原子核を結ぶ線とx軸が成す角\\ b : 衝突パラメータ（x軸からの距離）$

$r,\varphi$ は以下の関係を満たします。

$\frac{1}{r}=\Bigl(\frac{\cos (\varphi +\delta)}{\sin (\pi +\delta)}-\frac{1}{\tan\delta}\Bigr) / b$

ここで $\tan\delta$ は

$\tan\delta=\frac{4\pi\epsilon_0mbv_0^2}{qQ}$

で定義されます。ほかの記号の意味については、下記のとおりです。

$\epsilon_0 : 真空の誘電率\\ m : アルファ粒子の質量\\ v_0 : 遠方での速さ\\ q : アルファ粒子の電荷\\ Q : 原子核（質量無限大）の電荷$

真空の誘電率はえらい小さい値なので単位をメートルからピコメートルに変換して

$\tan\delta=\frac{4\pi mbv_0^2}{qQ}\times 8.858$

を代わりに使います。

原点に原子核をおいて、原子核へ向かう粒子の初期座標 $(x, b)$ を指定します。すると $\varphi$ の初期値が

$\varphi_0=\arctan (b/x)$

で決まります。また、散乱角は次で定義されます。

$\theta = 2\arctan \Bigl(\frac{1}{\tan\delta}\Bigr)$

よって、粒子の軌道は $\varphi$ を $\varphi_0$ から $\theta$ まで変化させた場合の極座標 $(r,\varphi)$ をプロットすればよいことになります。実際には以下の変換を行ってからプロットします。

$x=r\cos\varphi\\ y=r\sin\varphi$

プログラム

以下のプログラムを orbits.sage として保存します。

m=1 # mass of alpha particle
q=1 # charge of alpha particle
Q=5 # charge of atomic neuleus
v0=1 # velocity of alpha particle

def orbits(x,b,N):
    l=[]
    K = 8.858*4*pi*m*b*v0^2/q/Q
    delta = arctan(K)
    theta = 2*arctan(1/K) # scattering angle
    phi0 = pi - arctan(abs(b/x))
    
    for phi in sxrange(phi0,theta,-0.01):
        r = b / (cos(phi+delta)/sin(pi+delta) - 1/K)
        x = r*cos(phi)
        y = r*sin(phi)
        l.append((x,y))

    return l

SageShell 起動して、プログラムを読み込みます。

sage: load('orbits.sage')

衝突パラメータ b を 0.01 から 0.5 まで 0.02 刻みで変化させ、軌道を line で可視化します。

sage: l = [orbits(-1.5,b,10) for b in sxrange(0.01,0.5,0.02)]
sage: sum([line(l[i],xmin=-1.5,xmax=1.5,ymin=0,ymax=1) for i in range(len(l))])
Launched png viewer for Graphics object consisting of 25 graphics primitives

結果は以下のとおりです。

f:id:riverta1992:20210110170028p:plain

次に b が負になる場合も考えます。座標 (x,y) と (x,-y) の軌道は同じになるので、orbits.sage の修正は b に関する部分だけで済みます。

m=1 # mass of alpha particle
q=1 # charge of alpha particle
Q=5 # charge of atomic neuleus
v0=1 # velocity of alpha particle

def orbits(x,b,N):
    l=[]
    K = 8.858*4*pi*m*abs(b)*v0^2/q/Q
    delta = arctan(K)
    theta = 2*arctan(1/K) # scattering angle
    phi0 = pi - arctan(abs(b/x))
    
    for phi in sxrange(phi0,theta,-0.01):
        r = abs(b) / (cos(phi+delta)/sin(pi+delta) - 1/K)
        x = r*cos(phi)
        y = r*sin(phi)*b/abs(b)
        l.append((x,y))

    return l

修正したプログラムをロードして、b が -0.8 から 0.8 までの場合の軌道をプロットします。

sage: load('orbits.sage')
sage: l = [orbits(-1.5,b,10) for b in [-0.8,-0.5,-0.3,-0.1,-0.07,-0.05,-0.03,-0.01,0.01,0.03,0.05,0.07,0.1,0.3,0.5,0.8]]
sage: sum([line(l[i],xmin=-1.5,xmax=1.5,ymin=-1,ymax=1,linestyle="--") for i in range(len(l))])+point((0,0),size=30,color="black")+text("Atomic neucleus",(0.35,0.07),color="black")
Launched png viewer for Graphics object consisting of 18 graphics primitives