2.10 原子核の崩壊：問題 2.8 原子核の崩壊 a

前回までは単一核種の系を扱ってましたが、今回からは複数の核種の系、すなわちある核種が放射性崩壊してできた核種がさらに放射性崩壊するような系を見ていきます。

放射性崩壊する2核種の系

最初の種類の原子核の個数を ${ N_1(t) }$ 、崩壊によってできる放射性原子核の個数を ${ N_2(t) }$ とし、同様にそれぞれの崩壊定数を ${ \lambda_1,\, \lambda_2 }$ とします。　このとき、 ${ N_1(t),\,N_2(t) }$ の満たす微分方程式は以下のようになります：

　　 ${ \displaystyle \begin{align*} \frac{dN_1}{dt} &= -\lambda_1 N_1 \\[4mm] \frac{dN_2}{dt} &= \lambda_1 N_1 - \lambda_2 N_2 \end{align*} }$

この系は解析的に解けて、初期条件が ${ N_1(0) = N_{10},\,N_2(0) = N_{20} }$ のとき

　　 ${ \displaystyle \begin{align*} N_1(t) &= N_{10}\,e^{-\lambda_1 t} \\[2mm] N_2(t) &= \frac{\lambda_1}{\lambda_2-\lambda_1}N_{10}\,e^{-\lambda_1 t} + \left(N_{20} - \frac{\lambda_1}{\lambda_2-\lambda_1}N_{10}\right)\,e^{-\lambda_2 t} \end{align*} }$

特に、 ${ N_{20} = 0 }$ のとき、つまり最初2種目の原子核がないときは

　　 ${ \displaystyle \begin{align*} N_1(t) &= N_{10}\,e^{-\lambda_1 t} \\[2mm] N_2(t) &= \frac{\lambda_1}{\lambda_2-\lambda_1}N_{10}\left(e^{-\lambda_1 t} - e^{-\lambda_2 t}\right) \end{align*} }$

となります。　導出はこちらを参照。

問題

a.

⁷⁶Kr は電子を捕獲して半減期14.8時間で ⁷⁶Br に崩壊する

⁷⁶Br は電子を捕獲するか陽電子を放出して半減期16.1時間で ⁷⁶Se に崩壊する

試料には最初に純粋の ⁷⁶Kr が1g含まれていたとする

プログラムを修正して ⁷⁶Kr と ⁷⁶Se の量の1週間の時間変化を計算せよ。
b.

²⁸Mg はベータ線を出して半減期21時間で ²⁸Al に崩壊する*1

²⁸Al は陽電子を放出して半減期2.31分で ²⁸Si に崩壊する

時間の単位に分を使うと、ある程度の量の ²⁸Mg が崩壊するまでには多数回の反復計算が必要になる。　計算を速めるには、どのような簡単化の過程を置いたらよいか。
c. ²¹¹Rn は、下図に示すように、2つに枝分かれした過程で崩壊する。　必要ならばなんらかの近似をし、各同位体元素の量を時間の関数として計算せよ。　試料は初め 1μg の ²¹¹Rn からなっていたとする。

今回は問題 a をします。　問題 a に出てくる元素をもう少し詳しく見ると下表の通り：

元素記号	元素名	原子番号	周期・族	Wikipeida
Kr	クリプトン (krypton)	36	第4周期18族	wikipedia:クリプトン
Br	臭素 (bromine)	35	第4周期17族	wikipedia:臭素
Se	セレン (selenium)	34	第4周期16族	wikipedia:セレン

物理系

問題 a の物理系の簡単なドメイン分析をすると

状態変数
- 時間 ${ t }$ ： double
- Kr, Br, Se の核数 ${ N_\textrm{Kr},\,N_\textrm{Br},\,N_\textrm{Se} }$ ：それぞれ double
シミュレーション・パラメータ
- Kr, Br の崩壊定数 ${ \lambda_\textrm{Kr},\,\lambda_\textrm{Br} }$ ：それぞれ double

ってな感じですね。　今回は時間 ${ t }$ を BigDecimal ではなく double で定義してます。　これは後ほど。　で、この系が従う時間発展の微分方程式は

　　 ${ \displaystyle \begin{align*} \frac{dN_\textrm{Kr}}{dt} &= -\lambda_\textrm{Kr} N_\textrm{Kr} \\[4mm] \frac{dN_\textrm{Br}}{dt} &= \lambda_\textrm{Kr} N_\textrm{Kr} - \lambda_\textrm{Br} N_\textrm{Br} \\[4mm] \frac{dN_\textrm{Se}}{dt} &= \lambda_\textrm{Br} N_\textrm{Br} \end{align*} }$

となります。　Se の核数の微分方程式も載せましたが大丈夫でしょうかね？

時間の単位

2つの核種の崩壊定数が異なるので、半減期も違ってきます。　半減期が異なる場合、短い方の半減期と同程度の大きさの単位で時間を測ると都合がよいとのことなので、シミュレーションを書く前に「2.10 原子核の崩壊：問題 2.7 単一核種の崩壊 c」でやったように時間の単位を変更してやりましょう*2。　問題より、各核種の半減期は

　　 ${ \displaystyle \begin{align*} T_{1/2\, \textrm{Kr}} &= 14.8 \,\textrm{h}, & T_{1/2\, \textrm{Br}} = 16.1 \,\textrm{h} \end{align*} }$

なので、短い方の Kr の半減期を単位にして時間を測りましょう。　その時間を ${ \tau }$ （タウ）とすると

　　 ${ \displaystyle \begin{align*} \tau &= \frac{t}{T_{1/2\,\textrm{Kr}}} = \frac{\lambda_\textrm{Kr}}{\log 2}t & \left(\because T_{1/2 \, \textrm{Kr}} = \frac{\log 2}{\lambda_\textrm{Kr}}\right) \end{align*} }$

と定義されます。　 ${ t }$ から ${ \tau }$ への書き換えは、以前行ったのと同じような式の書き換えなので途中は省略して、結果だけを書くと

　　 ${ \displaystyle \begin{align*} \frac{dN_\textrm{Kr}}{d\tau} &= -N_\textrm{Kr} \log 2 \\[4mm] \frac{dN_\textrm{Br}}{d\tau} &= \left(N_\textrm{Kr} - \lambda N_\textrm{Br}\right)\log 2 \\[4mm] \frac{dN_\textrm{Se}}{dt} &= \lambda N_\textrm{Br}\log 2 \end{align*} }$

となります。　ただし

　　 ${ \displaystyle \begin{align*} \lambda = \frac{\lambda_\textrm{Br}}{\lambda_\textrm{Kr}} = \frac{T_{1/2\,\textrm{Kr}}}{T_{1/2\,\textrm{Br}}} \end{align*} }$

です。　ふむ、結構綺麗にまとまりましたね。

シミュレーション・スクリプト

以上を踏まえてシミュレーション・スクリプト。　シミュレーションの反復条件は1週間と指定されてますが、Kr の半減期で測った時間ではキレイな数値にならないので、今までのように時間を BigDecimal にする意味はないでしょう。　したがって、時間も double 型にします。　また、tau と書くのは面倒なので、この tau を改めて t としてスクリプトを書くことにします。

@GrabResolver('https://github.com/waman/groops-core/tree/master/repo')
@Grab('org.waman.groops:groops-core:0.2-beta')

import org.waman.groops.builder.SimulationBuilder
import org.waman.groops.simulation.system.*

double tKr = 14.8d, tBr = 16.1d    // 半減期
double aweek = 24 * 7 / tKr    // 一週間

new SimulationBuilder().simulation(type:'composite'){    // 複合シミュレーション
    simulation{
        system(new RadioactiveDecaySystem(), dt:0.01d, lambda:tKr/tBr)
        init(t:0d, nKr:1d, nBr:0d, nSe:0d)

        count(total:10)    // 10回の反復
    }

    // これらの反復条件、データ出力は10回の反復ごとに行われる（複合シミュレーション）
    continueWhile{ it.t <= aweek }
    console(outputTiming:FINAL)
    png(fileName:'problem2_8a.png', overwrite:true){
        chart(title:'Problem 2.8 a', domainLabel:"Time [$tKr hour]", rangeLabel:'Mass [g]'){
            point(label:'Kr', x:'t', y:'nKr')
            point(label:'Br', x:'t', y:'nBr')
            point(label:'Se', x:'t', y:'nSe')
        }
    }
}.simulate()

class RadioactiveDecaySystem extends AnnotationObservableSystem{

    @State double t, nKr, nBr, nSe
    @Parameter double dt, lambda

    static double LOG2 = Math.log(2d)

    // 今回はオイラー法で。
    @Override
    void evolve(Map<String, Object> params) {
        double dnKr =          nKr * LOG2 * dt
        double dnBr = lambda * nBr * LOG2 * dt

        nKr -= dnKr
        nBr -= dnBr - dnKr
        nSe += dnBr

        t += dt
    }
}