放射性崩壊と乱数の不思議な関係 ~計算物理学I (朝倉書店)を参考にpythonを使って~

放射性原子核の崩壊のシミュレーション

　前回と前々回の記事に引き続き、今回の記事でも乱数を使用したプログラム、放射性崩壊のシミュレーションを作成します。ランダムな過程なのになぜか規則性が見えてくるのが面白いです。

放射性崩壊とは

　さて、放射性崩壊とはどんなものでしょうか。参考にしている計算物理学I(朝倉書店)のp.79には以下のように記載されています。

　原子核が外部からの影響無しに崩壊して他の種類の原子核に変わる物理的な過程が放射性崩壊である. 1個の原子核がある時間内に崩壊する確率は一定だが, いつ崩壊するかというのはランダムな事象である. 自発的に起きる崩壊(自発崩壊ともいう)では, どの原子核についてもそれが崩壊する時刻は常にランダムであるし他の原子核から影響を受けないから, その原子核がどれほど長く崩壊せずにいようと, 他の原子核が崩壊していようと, 他の原子核が崩壊していようと, 崩壊する確率は同じである.

　この文章のポイントは、それぞれの原子核が崩壊するかは外部の影響を受けないということです。各原子核間の相互作用などの外的な影響を全く考慮しなくてもよく、単位時間あたりに各原子核が崩壊したかどうかだけを数えていけばシミュレーションができあがるというわけです。また、放射性崩壊は確率的な事象なので、それぞれの原子核が単位時間の中で崩壊するかどうかの判定は乱数を用いて実施します。

放射性崩壊のシミュレーション

　さて、今回は数学的なことはさておいて、早速シミュレーションのコードから見ていきましょう。自作のコードは以下です。閾値$λ$を設定してそれよりも乱数の結果が高いか低いかで放射性崩壊が起こったかどうかを判定しています。

radioactive_decay

import random
import math
import numpy as np

#閾値lambdaN、最大経過時間time_max、原子核の最大個数のリストlist_num_max
#乱数の種seedを渡すことで、放射性崩壊のシミュレーション結果を返す。
#原子核の個数はあたりの対数表示
def radioactive_decay(lambdaN,time_max,list_num_max,seed):
    #乱数の種の指定
    random.seed(seed)

    #結果を格納するためのリストを作成。行列形式を使用していることに注意。
    num_max_no = list_num_max.shape[0]
    list_number = np.zeros((num_max_no,time_max+1))
    list_deltan = np.zeros((num_max_no,time_max+1)) 
    list_time = np.arange(0,time_max+1,1)

    #原子核の各最大個数に対してのシミュレーションの実施
    for i in range(num_max_no):
        number = list_num_max[i]
        list_number[i,0] = math.log(number)
        #単位時間あたりの崩壊個数の計算のループ
        for time in range(1, time_max+1):
            deltan = 0
            #単位時間での各原子一つ一つの崩壊の有無を調べるループ
            for atom in range(1,number+1):
                #0~1の間の乱数を発生し、閾値lambdaNを下回ったら
                #その原子は崩壊したものとして数える
                decay = random.random()
                if (decay < lambdaN):
                    deltan += 1
                else :
                    deltan += 0
            #単位時間後のまだ崩壊していない原子数を計算
            number = number - deltan
            #リストに単位時間後のまだ崩壊していない原子核数の対数を保存
            if number == 0:
                list_number[i,time] = 0
            else:
                list_number[i,time] = math.log(number)
            #単位時間で崩壊した原子核数の対数を保存
            if deltan == 0:
                list_deltan[i,time] = 0
            else:
                list_deltan[i,time] = math.log(deltan)
                
    return list_number, list_deltan, list_time

放射性崩壊のシミュレーション結果

　$λ = 0.005$とし、原子核の個数$N$を$10,100,1000,10000,100000$とした時のシミュレーション結果が以下です。あまり時間経過していない時(崩壊原子核が少ない時)は全て直線的なふるまいをしています。傾きも同じようです。

　実際のところ、各直線(?)の傾きはどんな値になるのでしょうか。$N = 10000, 100000$のときのグラフで、平均二乗誤差を最小にするような直線モデルでフィッティングしてみました。以下に示します。

$\lambda = 0.005, N = 10000$の時、
フィッティング関数$y = -0.005 x + 9.242$で、平均二乗誤差の平方根 $0.034$

$\lambda = 0.005, N = 100000$の時、
フィッティング関数$y = -0.005 x + 11.525$で、平均二乗誤差の平方根 $0.021$

　結果を見ると、閾値$\lambda$と直線モデルの傾きの絶対値は同じ値になるようです！本当でしょうか。

　確認のために、$\lambda = 0.01$として、上記と同様のシミュレーション、フィッティングを実施してみます。

$\lambda = 0.01, N = 10000$の時、
フィッティング関数$y = -0.010 x + 9.148$で、平均二乗誤差の平方根 $0.041$

$\lambda = 0.01, N = 100000$の時、
フィッティング関数$y = -0.010 x + 11.489$で、平均二乗誤差の平方根 $0.028$

したがって、閾値$\lambda = 0.01$でも、閾値とフィッティング直線の傾きの絶対値と一致しているようです！ただ放射性崩壊に閾値を設けただけなのに、崩壊の仕方に何かしらの法則性が生じたようです！驚き！

放射性崩壊のシミュレーションのまとめ

　さて、シミュレーション結果を見たところでまとめたいと思います。シミュレーションの結果、崩壊していない原子核数が十分に多い時、対数$\log(N(t))$は$t$の一次関数となることが分かりました。$N(0) = N$とすれば
　　　　　　　　$N(t) = Ne^{-\lambda t}$
と書けます。これは、有名な放射性崩壊の減衰を表す式と同じですね！閾値$\lambda$は実は崩壊定数と呼ばれるもので寿命$\tau$の逆数です。非常に簡単なモデルで放射性崩壊をなかなか良い感じに再現することができたようです。面白いですね～。
　また、数学的には、原子核数$N$が十分に大きければ放射能の差分方程式が微分方程式になるので指数関数的な減衰となるそうです。実際シミュレーション結果と理論が合致するのも面白いですね～。不思議でおもしれぇ～。また、何かやってみます。

参考書籍
R.H.Landau・他 (2018)『実践Pythonライブラリー 計算物理学I -数値計算の基礎/HPC/フーリエウェーブレット解析』 (小柳義夫・他訳) 朝倉書店