放射性崩壊の離散モデル

おことわり

この記事は、思いついたのでやってみたがうまくいかなかったというものです。なにか情報やちょっとしたアイデアなどありましたらご教示いただけると幸いです。

はじめはよくある微分方程式モデル

放射性崩壊はふつう微分方程式で記述される。放射性同位体の原子核数$N$の時間変化は、
\begin{aligned} \frac{dN}{dt}=-\lambda N \qquad\therefore\, N=N_0\exp(-\lambda t) \end{aligned}
と表される。ここで$N_0$は初期状態における原子核数、$\lambda$は崩壊定数で、単位時間で崩壊する放射性原子核の割合を意味する。よく用いられる量として、半減期$t_{1/2}$と平均寿命$\tau$がある。前者は原子核数が$\dfrac{N_0}{2}$となる時間、後者は$\dfrac{N_0}{e}$となる時間。
\begin{aligned} N_0\exp(-\lambda t_{1/2})=\frac{N_0}{2} \qquad&\therefore\, t_{1/2}=\frac{\log 2}{\lambda} \\ N_0\exp(-\lambda\tau)=\frac{N_0}{e} \qquad&\therefore\, \tau=\frac{1}{\lambda} \end{aligned}
このモデルは計算も扱いも簡単で使い勝手がよい一方で、原子核の数って離散的なのでは？という気もする。そこで離散モデルを組んでみた。のはいいものの、計算がわたしの手には負えなくなりました、という話です。

離散モデルの構築

初期時刻$t=0$に$N_0$個の放射性原子核があるとする。これらは時刻$t$から$t+1$の間に確率$\lambda$で崩壊する ($0<\lambda<1$)。時刻$t$における未崩壊の放射性原子核の数を$N(t)$とする。
ある放射性原子核が時刻$\tau$で崩壊したとき、その原子核の寿命は$\tau$であるとする。
はじめて放射性原子核の数が$\dfrac{N_0}{2}$を下回った時刻を半減期とし、$t_{1/2}$とする。なお、$N_0$は偶数とする。
全ての放射性原子核が崩壊し終えた時刻を崩壊の長さとよび$L$とする。

放射性原子核の個数を数えるようにモデル化した。時刻も離散的に扱う。またこれによって、現実的とは言えないものの『すべての原子核が崩壊する』という状況が考えられる。このモデルで求めたいものは、

• 寿命$\tau$の期待値
• 半減期$t_{1/2}$の期待値と最頻値
• 崩壊の長さ$L$の期待値と最頻値

平均寿命と半減期については微分方程式モデルと比較出来たらなおよい。

寿命$\tau$とその期待値

ある放射性原子核の寿命が$\tau$である確率を$p_\tau(\tau)$と書く。
\begin{aligned} p_\tau(\tau)=\lambda(1-\lambda)^{\tau-1} \end{aligned}
だから、$\tau$の期待値$E_\tau$は$\disp E_\tau=\lim_{m\to\infty}\sum_{\tau=1}^m \tau\cdot p_\tau(\tau)=\lim_{m\to\infty}\sum_{\tau=1}^m\tau\lambda(1-\lambda)^{\tau-1}$であって、この極限の中の総和$S$について
\begin{aligned} S&=\sum_{\tau=1}^m\tau\lambda(1-\lambda)^{\tau-1} \\ (1-\lambda)S&=\sum_{\tau=1}^m\tau\lambda(1-\lambda)^\tau=\sum_{\tau=2}^m(\tau-1)\lambda(1-\lambda)^{\tau-1} \\ \therefore\,\,\lambda S&=\lambda+\sum_{\tau=2}^m\lambda(1-\lambda)^{\tau-1}=\sum_{\tau=1}^m\lambda(1-\lambda)^{\tau-1} \\ \therefore\,\, S&=\sum_{\tau=1}^m (1-\lambda)^{\tau-1}=\frac{1-(1-\lambda)^m}{1-(1-\lambda)}\to \frac{1}{\lambda} \quad(m\to\infty) \end{aligned}
より、$E_\tau=\dfrac{1}{\lambda}$と求まった。これは微分方程式モデルでの議論と整合するんですね。面白い発見でした。

崩壊の長さ$L$について

$N_0$個の放射性原子核について、その崩壊の長さが$L$となる確率を$p_{L:N_0}(L)$と書く。崩壊の長さが$L$となるのは、$N_0$個すべての放射性原子核の寿命が$L$以下である場合から、すべての寿命が$L-1$以下である場合を除いた場合となる。一方、ある放射性原子核の寿命が$\tau$以下となるのは、時刻$\tau$までに崩壊が起こらない場合の余事象だから、その確率は$1-(1-\lambda)^\tau$。したがって、
\begin{aligned} p_{L:N_0}(L)=(1-\mu^L)^{N_0}-(1-\mu^{L-1})^{N_0} \end{aligned}
である。ただし$\mu=1-\lambda$とした。これより期待値$E_L$は、$\disp E_L=\lim_{m\to\infty}\sum_{L=1}^m Lp_{L:N_0}(L)=\lim_{m\to\infty}\sum_{L=1}^m L\Bigl[(1-\mu^L)^{N_0}-(1-\mu^{L-1})^{N_0}\Bigr]$であって、この極限の中の総和について
\begin{aligned} &\sum_{L=1}^m L\Bigl[(1-\mu^L)^{N_0}-(1-\mu^{L-1})^{N_0}\Bigr] \\ =&\sum_{L=1}^m L(1-\mu^L)^{N_0}-\sum_{L=0}^{m-1} (L+1)(1-\mu^L)^{N_0} \\ =&\,m(1-\mu^m)^{N_0}-\sum_{L=1}^{m-1} (1-\mu^L)^{N_0} \\ =&\,m(1-\mu^m)^{N_0}-\sum_{L=1}^{m-1}\sum_{k=0}^{N_0}\combi{N_0}{k}1^{N_0-k}(-\mu^L)^k \\ =&\,m(1-\mu^m)^{N_0}-\sum_{L=1}^{m-1} \left(1+\sum_{k=1}^{N_0}\combi{N_0}{k}(-1)^k(\mu^k)^L\right)\\ =&\,m(1-\mu^m)^{N_0}-(m-1)-\sum_{k=1}^{N_0}\combi{N_0}{k}(-\mu)^k\left(\sum_{L=1}^{m-1} (\mu^k)^{L-1}\right) \\ =&\,m\left(1+\sum_{j=1}^{N_0}\combi{N_0}{j}1^{N_0-j}(-\mu^m)^j\right)-(m-1)-\sum_{k=1}^{N_0}\combi{N_0}{k}(-\mu)^k\frac{1-(\mu^k)^{m-1}}{1-\mu^k} \\ =&\, 1+\sum_{j=1}^{N_0}\combi{N_0}{j}(-\mu^m)^j-\sum_{k=1}^{N_0}\combi{N_0}{k}(-\mu)^k\frac{1-(\mu^k)^{m-1}}{1-\mu^k} \\ \to&\, 1-\sum_{k=1}^{N_0}\combi{N_0}{k}\frac{(-\mu)^k}{1-\mu^k}=E_L \quad (m\to\infty) \end{aligned}
と求まった。
わたしはこの総和が求められずにいます。$p_{L:N_0}$と$E_L$がもっともらしいことは、例えば Desmos で描画してみるとわかります。(ただし$N_0$をあまり大きくしすぎると$E_L$の計算で数値計算による誤差がひどくなります。)最頻値も興味のあるところですが、微分はできても離散的な$L$に対して$p_{L:N_0}$の最大値を求めることはうまくできませんでした。

半減期$t_{1/2}$について

時刻$t_{1/2}-1$までに崩壊した原子核の個数を$i$、時刻$t_{1/2}$に崩壊した原子核の個数を$n$とすると、このようなことが起こる確率は、
\begin{aligned} \combi{N_0}{i}\cdot p_{L:i}(t_{1/2}-1)\times\combi{N_0-i}{n}\lambda^n \end{aligned}
となる。ただし、便宜上$p_{L:i}(0)=0$とする。$N_0$個の放射性原子核について、その半減期が$t_{1/2}$となる確率を$p_{1/2:N_0}(t_{1/2})$と書く。これは、$t_{1/2}>1$の場合には
\begin{aligned} p_{1/2:N_0}(t_{1/2})&=\sum_{i,n}\combi{N_0}{i}\cdot p_{L:i}(t_{1/2}-1)\times\combi{N_0-i}{n}\lambda^n \\ &=\sum_{i=0}^{N_0/2-1}\sum_{n=N_0/2-i}^{N_0-i} \combi{N_0}{i}\cdot\Bigl[(1-\mu^{t_{1/2}-1})^i-(1-\mu^{t_{1/2}-2})^i\Bigr]\cdot\combi{N_0-i}{n}\lambda^n \end{aligned}
であり、$t_{1/2}=1$の場合には
\begin{aligned} p_{1/2:N_0}(1)=\sum_{n=N_0/2}^{N_0} \combi{N_0}{n}\lambda^n \end{aligned}
離散的に考えたことで、必ずしもちょうど半分崩壊するとは限らなくなり、それによって複雑極まりない場合分けと総和計算に直面しています。正直な話、先が見えず疲れてしまったため、ここで一度打ち止めます。

はじめにも書きましたが、この件についてなにか情報やアイデアがありましたらお知らせいただけると大変うれしいです。