原子核のアルファ崩壊系列

2015 年 2 月
福井大学 工学部 物理工学科
松岡 弘和

目次
第1章

序論

3

第2章
2.1

2.2
2.3

自然の階層構造素粒子 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1.1 フェルミ粒子 . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1.2 ボース粒子 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1.3 ハドロン . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
原子核 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
原子 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4
4
4
5
6
7
7

第3章
3.1
3.2

原子核
構成要素 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
原子核の質量公式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

8
8
9

第4章
4.1
4.2

原子核の放射性崩壊
放射性同位体の性質 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
α・β・γ 崩壊 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2.1 質量数を変える崩壊様式 . . . . . . . . . . . . . . . . . . . . . .
4.2.2 核種を変えるが質量数を変えない崩壊様式 . . . . . . . . . . . .
4.2.3 核種を変えない崩壊様式 . . . . . . . . . . . . . . . . . . . . . .
放射性崩壊系列 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3.1 α 崩壊系列 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3.2 β 崩壊系列 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

11
11
11
11
12
12
13
13
16

第5章
5.1
5.2
5.3
5.4
5.5

放射性崩壊の時間発展を記述する連立常微分方程式とその解1 種類の原子核の崩壊を表す微分方程式とその一般解 . . . . . . . . . .
1 種類の親核から 1 種類の娘核への崩壊 . . . . . . . . . . . . . . . . . .
親核の崩壊先が 2 つに分岐する場合 . . . . . . . . . . . . . . . . . . . .
親核が娘核, 娘核が孫核へと 2 段階にわたり崩壊する場合 . . . . . . . .
親核が複数の経路を経て同一の娘核に崩壊する場合 . . . . . . . . . . .

21
21
22
23
25
28

第6章
6.1
6.2
6.3

計算結果
数値計算プログラム作成の概略 . . . . . . . . . . . . . . . . . . . . . .
α 崩壊系列の簡単な計算例 . . . . . . . . . . . . . . . . . . . . . . . . .
超長寿命核種を起点とする放射平衡 . . . . . . . . . . . . . . . . . . . .

33
33
34
38

4.3

1

6.4
第7章

短寿命な超重核を起点とする α 崩壊系列 . . . . . . . . . . . . . . . . .

46

まとめ

51

参考文献

52

謝辞

53

付録

54

2

第1章

序論

物理工学科で学んだことの集大成としての卒業研究のテーマとして私は原子核の α（アルファ）崩壊系列の時間発展の計算という課題を選んだ。この課題は, 数学の観点からは, 時間をかけて学んできた微分積分や線形代数の知識を活用できる線形 1 階連立
常微分方程式の解法という問題とみることができるし, 情報技術という観点からは, 私が比較的得意分野だと感じている教科である C 言語でのプログラミングの応用課題となっている。工学的視点からは 2011 年の東日本大震災で発生した東京電力福島第一原
子力発電所の炉心溶融事故以降社会的関心の高まった放射性物質の話題である。そして, 物理学としては私が興味をもつ原子核や素粒子といった物質のミクロな構造に奥では関連している課題として興味の持てるものである。こういった理由からこの課題を
私の卒業研究の課題として選んだのである。
本研究は大沼・森近両氏による 2011 年度の卒業研究 [1] で開発された原子核の放射性
崩壊系列の時間発展を連立常微分方程式を求める C 言語プログラムを発展させる。具
体的には, 半減期や崩壊分岐比等の大量のデータをプログラム中に記述する方法を改良する。大沼・森近両氏は, このプログラムを質量数 90〜140 の核分裂生成物の β（ベータ）崩壊系列に対して使用したが, 本研究ではアクチナイド核の α 崩壊系列に対して使用する。β 崩壊系列については, 崩壊系列の下流に進むにつれて核種の半減期が長くなることが原因で, 系列の全放射線量はベクレル単位で時間に反比例することが大沼・森近両氏により見出された。本研究では, α 崩壊系列で, 系列の開始点として, 核図表上で周囲に配置される核より顕著に半減期の大きな核種 (238 U 等) を使うと, これは放射平衡と呼ばれる, β 崩壊系列とは異なった崩壊の様態を示すことがよく知られているが, まず, そのことを確かめる。次に, 崩壊の開始点を, もっと重く不安定な（半減期の短い）核にとると系列の核種分布や総放射線量はどのような挙動を示すかを調べる。例として日本の理化学研究所が発見した原子番号 113 番元素の放射性崩壊系列の経時変化等を示す。
本論文では, まず, 素粒子, 原子核, 原子核の放射性崩壊についての基礎知識を整理してまとめを示す。次に定数係数連立 1 階常微分方程式の解析的厳密解について詳しく説明し, 典型的な状況設定をいくつか行って, その解の挙動をグラフ化し, 性質について論じる。最後に, 研究の主目的である α 崩壊系列の計算結果を示し, その挙動について論じる。付録として本研究のために新規に開発した C 言語プログラムや, 改良した C 言語プログラムのプログラム・リストを採録する。

3

第2章
2.1

自然の階層構造

素粒子

素粒子 (elementary particle) とは, 現在発見されている中では最小の物質構造であり,基本粒子とされているものをさす言葉である。
18 世紀にドルトンにより初めて (哲学的でなく) 科学的な原子論が成立し, その後 19
世紀末までは基本的な粒子は原子と考えられていた。
しかし, 1897 年のトムソンによる電子の発見および 1911 年のラザフォードによる原子核の発見により, 原子が物質構造の最小単位ではないということが発見された。その後の陽子と中性子の発見, さらに湯川秀樹の中間子論等によって素粒子の存在が確定され今日の素粒子物理学に至っている。
さて, 素粒子は大きくフェルミ粒子 (Fermion), ボース粒子 (Boson) にわけられる。フェルミ粒子及びボース粒子とは, その多粒子状態の波動関数が粒子座標の交換に対して完
全反対称及び完全対象である粒子として定義される。その重要な帰結として, フェルミ粒子は一つの状態には一つの粒子しか入れない。これをパウリの禁止則とよぶ。一方,
ボース粒子は何個でも同一の状態に入ることができる。特に, 非常に多くの個数の粒子が一状態に入った場合をボース・アインシュタイン凝縮とよぶ。
実験的には, スピンの大きさが ℏ の半整数 ( 21 , 32 , 52 ,・・・) 倍であればフェルミ粒子, 整数(0,1,2,・・・) 倍であればボース粒子と判定できる。この関係は「スピンと統計の関係」と呼ばれ, 適当な仮定の下に相対論的場の理論に基づき証明することができる。
また, クォークの複合粒子であるハドロン (hadron) については別途 2.1.3 で述べる。

2.1.1

フェルミ粒子フェルミ粒子 (Fermion) は物質の主要な構成要素であり, クォーク (quark) とレプトン(lepton) に大別される。以下でその 2 つについて順次述べる。
クォークは電子の ± 13 または ± 23 倍の電荷を持つ粒子である。なお,1 個の電子の持つ電荷の −1 倍を素電荷と呼び, 記号 e で表す。なお, 同じ記号を用いる数学定数の自然対数の底 (= 2.718 · · · ) とは無関係である。その最新の測定値は e = 1.602176565(35) × 10−19
クーロンである (丸括弧内の「35」は末尾の桁の「65」の含む誤差の標準偏差を表す[2]）。この記号 e を用いると, クォークの電荷は ± 31 e または ± 23 e である。
クォークは up, down, strange, charm, bottom, top の 6 種類が存在することが実験的に確認されており, それぞれ, u, d, s, c, b, t の記号で表される。クォークは似た性質を持
つもの同士の 2 つずつの組み合わせに分けることができる。この組み合わせは「世代(generation)」とよばれ, 表 2.1 に示されている。今現在においてクォークを単独で取り4

表 2.1: フェルミオンの一覧表。なお, これらのほかに, それぞれの反粒子が存在する。
名称
記号
世代
電荷
質量

クォーク

レプトン

up

u

第 1 世代+2e3

1
2.3+0.7−0.5 MeV

down
charm
strangetopbottomdcstb

第 1 世代
第 2 世代
第 2 世代
第 3 世代
第 3 世代

−1
e
3
+2
e
3
−1
e
3
+2
e
3
−1
e
3

4.8+0.7−0.3 MeV
1.275 ± 0.025GeV
95 ± 5MeV

電子
電子ニュートリノe
νe

第 1 世代
第 1 世代-1e00.510998928±0.000000011 MeV

ミュー粒子
ミュー・ニュートリノ
タウ粒子
タウ・ニュートリノ

µ
νµ
τ
ντ

第 2 世代
第 2 世代
第 3 世代
第 3 世代

-1e
0
-1e
0

105.6583715±0.0000035 MeV

1776.82 ± 0.16MeV

出すことはできていないが, 実験においてはジェット (そろった方向に発生する多数のハドロンの束) によって存在を確認されている。
レプトンはクォークとは異なり単独で粒子として存在する。レプトンと呼ばれるものの中には電子 (electron), 電子ニュートリノ (electron neutrino), ミュー粒子 (muon), ミュー・
ニュートリノ (muon neutrino), タウ粒子 (tauon, tau), タウ・ニュートリノ (tau neutrino)がある。それぞれ e, νe , µ, νµ , τ, ντ という記号で表される。また, クォークと同じように性質によって 2 つずつの組み合わせに分けることができる。この組み合わせを世代と呼ぶ。
レプトンの間には弱い相互作用が働き, 原子炉等における中性子や不安定原子核の β崩壊に関与する。また, ミュー粒子は宇宙から降りそそぐ放射線, すなわち宇宙線の 2 次粒子として, ニュートリノは太陽の中心部での核反応や超新星爆発に伴って放出されるものとしても観測することができる。ニュートリノについては岐阜県の神岡鉱山をはじめとして世界数箇所に観測施設が存在している。これらクォークやレプトンは基本粒子と考えられている。なお, フェルミ粒子の諸性質のまとめを表 2.1 の下半分に示す。

2.1.2

ボース粒子

粒子間に力が働くのは, 力を媒介する粒子が存在するからだと考えられている。この粒子がボース粒子 (Boson) であり, 前節のフェルミ粒子と同様に, これらのボース粒子
も基本粒子であるとされている。ボース粒子には, 電磁気力を媒介する粒子である光子
1eV とは, 電子ボルトとも呼ばれるエネルギーの単位であり, 1V の電位差において, １電気素量の粒子を加速させた時に得られる運動エネルギーのことである。なお, 1eV=1.602176565(35)×10−19 J となる。

5

名称
陽子
中性子

英名
protonneutron

パイ中間子

pion

表 2.2: 主なハドロンの一覧表
記号 電荷
質量
p
+e 938.272046±0.000021MeV
n0939.565379±0.000021MeV+
π+e139.57018±0.000MeV0π
0134.9766±0.0006MeV−
π-e
139.57018±0.000MeV(photon), 弱い力を媒介する粒子であるウィーク・ボソン (weak boson), 強い力を媒介する粒子であるグルーオン (gluon) が知られている。更に, 重力を媒介する粒子である重力子 (graviton) も存在すると考えられているが, 現在のところはまだ実験的に観測されたことがない。
また, ヒッグス粒子 (Higgs particle, Higgs Boson) とよばれるボース粒子も存在が仮定されてきたが, 2012 年に欧州原子核研究機構 (CERN) によって発見された粒子がこのヒッグス粒子であると考えられている。Higgs 粒子は物質のもつ質量の起源であると考えられている。

2.1.3

ハドロン

ハドロン (hadron) とは, 前述した基本粒子のうちのクォークが強い相互作用によって結びつくことで構成された複合粒子をさすことばであり, バリオン (baryon) とメソン
（またはメゾン, meson) の 2 つに大別することができる。
バリオンは重粒子ともよばれる。バリオンは 3 つのクォークにより構成され, 代表的なものとして陽子と中性子とが挙げられる。陽子は u クォーク 2 個と d クォーク 1 個から構成され, 中性子は u クォーク１個と d クォーク 2 個から構成される複合粒子だと考えられている。１個のバリオンを構成するクォークの組み合わせは, そのバリオンの総電荷量が常に電子の電荷の整数倍になるものに制限されている。そのようなクォークの可能な組み合わせで, クォーク概念が登場する以前の時代に素粒子であるとされていた粒子をすべて説明することができる。このことを陽子を例にとって確認してみると,陽子は電 荷が + 32 e の u クォークが 2 個と − 13 e の d クォークが１個から成るのだから, その総電荷量を計算してみると,

+2
−1
e×2+
e=1·e
3
3

(2.1)

という式の通り, ちゃんと素電荷 e の整数倍 (この場合は 1 倍) になっている。
メソンは中間子ともよばれる。メソンは 2 個のクォークにより構成され, 代表的なものとしてはパイ中間子 (pion) が挙げられる。パイ中間子も基本粒子であるクォークによって構成されており, 電荷が +e, 0, −e で 3 種類が存在し, それぞれ π+ , π0 , π− という記号で表される。主なハドロンの諸性質を表 2.2 にまとめておく。

6

2.2

原子核

原子核 (atomic nucleus) は, 2.1 節において既に述べたように, ラザフォードによって発見された。原子核は陽子と中性子から構成される。原子核は原子の中心に存在し, 原子の質量のほとんどを占める。例えば水素原子の場合, 99.95 ％を占める。水素の原子
核は陽子 1 個で構成されており, 直径は 10−15 m 程である。これは次節で述べる原子と比較すると 4 桁も小さい。原子核については 3 章で更に詳しく述べる。

2.3

原子

原子 (atom) は古代ギリシャ時代の「これ以上分解することができない基本粒子」という意味から生まれた言葉だが, 今日では原子核と（素粒子である）電子とに更に分解できることがわかっている。原子核の発見前は, 原子の構造のモデルとして, 1903 年にトムソンによるプラム・プディングモデル, 同じ年に長岡半太郎による土星型モデルが考案されていた。1911 年のラザフォードによる原子核の発見後は, ラザフォード本人による惑星モデルも考案されたが, 1913 年にボーアによって考案された前期量子論的な原子モデルが最も画期的であり, 今日でも半古典的量子化法の成功例として教育的な価値のある描像として活用されている。この原子モデルでは, 原子は原子核の周りを電子が回っていることにより形作られている。原子の大きさは 10−10 ｍ程度とされている。
この長さは, 非 SI 単位のオングストロームを使うと, ちょうど 1Å である。

7

第3章
3.1

原子核

構成要素

原子核 (atomic nucleus) は前述したように陽子 (proton) と中性子 (neutron) から構成されている。この事実は, 1932 年にチャドウィックによる中性子の発見により判明した。
陽子の個数を原子番号 (atomic number), 陽子と中性子の個数の和を質量数 (mass number) と呼ぶ。一般に, 質量数を A, 中性子の個数（中性子数, neutron number）を N, 原子
番号を Z という記号を用いて表す。なお, 陽子と中性子を総称して核子 (nucleon) という。これらの 3 数の間の関係を数式で表すと,

A=Z+N

(3.1)

となる。したがって, 原子核の種類は A, Z, N のうちの 2 つの数だけで指定することができ, 通常は質量数 A と原子番号 Z で指定する。そのように指定された原子核は, Z 個の陽子と A − Z 個の中性子よりなり, Ze の正電荷を持つ。
質量数 A と原子番号 Z によって指定されるそれぞれの原子核のことを核種とよぶ。
質量数 A, 原子番号 Z, 中性子数 N, 元素記号 X の核種は ZA XN のように表記する (脚注 1参照。）。
具体例として, 原子炉の燃料として利用されている質量数が 235 のウラン原子核を例にとると, ウラン原子は原子番号 92 なので, 中性子数は 235 − 92 = 143 であり, 235
92 U143という記号で表記される。
原子の化学的性質は, 陽子の個数（即ち原子番号）だけでほぼ決まってしまう。したがって, 陽子の個数が, 原子番号と呼ばれ, 元素名 (element name), 元素記号 (element
symbol) と一対一に対応するのである。元素を化学的性質を考慮して表に配置したもの
が周期表 (periodic table) である。今日までに 115 個の元素の存在が実験的に確認されている。
同じ元素の中でも原子核内の中性子数が異なる原子を同位体（isotope）とよぶ。なお, 中性子数が同じものを同調核（isotone）, 質量数が同じものを同重核（isobar）と呼ぶ。同位体の例を挙げると, 水素には 11 H0 , 21 H1 , 31 H2 の 3 つの同位体がある。同位体の中には安定して存在する安定同位体と放射性物質である放射性同位体がある。放射性同位体は, 時間とともに放射性崩壊を起こし, 放射線を出すとともに別種の原子核に変わる。放射性同位体は, 広く自然界に存在する天然放射性同位体と加速器や原子炉によって生成される人工放射性同位体に大別することができる。
1

但し, 元素記号の右下の添字は化学式では原子の個数を表すので読む人が意味を混同しないように注意して用いなければならない。

8

3.2

原子核の質量公式

原子核の質量公式 (mass formula) について述べる前に, まず, 質量とエネルギーの間に成り立つ関係, 即ち, アインシュタインによって見出された質量とエネルギーの等価性を表す下記の式が成り立つことを述べておく必要がある。

E = mc2(3.2)

ただし, c は光速であり, c = 299792458m·s−1 である。[2](3.2) 式により, 質量 m の物質はエネルギー E をもつと考えることができる。この式
は, 単に, それまでは別種の量と考えていたエネルギーと質量とを, 換算率を c2 だと定義して, 今後は比較してよいことにしよう, と提唱しているだけに過ぎないように思えるかもしれないが, 実はそれは正しい理解ではない。物質が反応をし, 別の物質になることを考えると,∆E = ∆mc2(3.3)となる。(3.3) 式は, 物質の反応で吸収・放出されるエネルギー ∆E と質量の変化 ∆m を関係付ける有用な法則である。
(3.2) 式を用い, 原子核の結合エネルギー（binding energy）を質量から求めるための式を導こう。結合エネルギーとは, 原子核をばらばらな核子の集まりに分解するのに必要な最小エネルギーのことである。ばらばらに離れて静止している Z 個の陽子と N 個の中性子のエネルギーを E1 とすると, (3.2) 式を用いてE1 = ZM p c2 + N Mn c2(3.4)となる。これら A 個の核子が集まって原子核を作ったときのエネルギーを E2 とすると,核子の質量を MA として,E 2 = M A c2
(3.5)
となる。このとき, 原子核となることで余るエネルギーが結合エネルギーであり, 記号
B であらわすと,E1 = E2 + B(3.6)∴ B = E1 − E2
= (Z M p + N Mn − MA,Z )c2 = ∆Mc2

(3.7)

が成り立つ。(3.7) 式における B が核種の全結合エネルギーである。B は正の値を取る。
また, ∆M を質量欠損と呼ぶ。結合エネルギー B を質量数 A で割ったもの, 即ち AB を核子１個あたりの結合エネルギーという。
(3.7) 式を用いると質量が既知の核種の B が計算できる。また, 核反応で放出・吸収されるエネルギーからも B が計算できる。ヴァイツゼッカーとベーテはこうして得られた多数の（現在では三千個にもなる）核種の B を N, Z, A の簡単な数式で表すことに成功した。それが原子核の質量公式である。以下が, Z 個の陽子と N 個の中性子を持つ原子の質量を求める式である。
2

MA,Z c2 = N Mn c2 + ZM p c2 + Zme c2 − av A + a s A 3 + ac
9Z2
(N − Z)2
δ
+
a
+ 1
a
1
4A
A3
A2

(3.8)

(3.8) 式において me は電子の質量を表す記号であり, 既に表 2.1 に記した通り,
me = 0.510998928 ± 0.000000011MeV(3.9)

である。パラメータ a s ,av ,ac ,aa ,δ の値は, 多数の原子核の束縛エネルギーへの最小２乗法によるフィッティングの結果, 下記の値が得られている。









δ=







av = 15.67MeV(3.10)

a s = 17.23MeV(3.11)

ac = 0.714MeV(3.12)

aa = 93.15MeV(3.13)− 11.2 MeV(ZとNが偶数の場合 (偶偶核))
0 MeV(Aが奇数の場合 (奇偶核))(3.14)+ 11.2 MeV(ZとNが奇数の場合 (奇奇核))(3.8) 式の第 1 項は体積項, 第 2 項は表面項, 第 3 項はクーロン項, 第 4 項は非対称項, 第 5
項は対結合エネルギー項とそれぞれよばれる。それぞれの項についての詳しい説明は
文献 [4] の 18 ページから始まる 2.3 節を参照されたい。

10

第4章
4.1

原子核の放射性崩壊

放射性同位体の性質3.1 で触れた放射性同位体は不安定な状態にあり, 粒子線や電磁波を放出して安定な同位体になろうとする性質を持つ。こういった性質を持つ原子核を含んだ物質を東日
本大震災以来よく聞くようになった放射性物質 (radioative material) と呼ぶ。この現象
を原子核の放射性崩壊 (radioactive decay) という。
一般的傾向として, 原子核は重くなればなるほど陽子同士のクーロン斥力が強くなる
ため不安定となり, α 崩壊と呼ばれる放射性崩壊を起こしやすくなる。また陽子と中性子の個数の差が大きくなるほど不安定となり, β 崩壊と呼ばれる放射性崩壊を起こしやすくなる。
放射性物質の質量が半分になるまでにかかる時間を半減期 (halflife) という。物質によって半減期は長短さまざまである。あまりにも半減期の長いものは放射性崩壊が起こらないと見なしてよい。

4.2

α・β・γ 崩壊

原子核の放射性崩壊は大きく 3 つに分けることができ, α 崩壊・β 崩壊・γ（ガンマ）崩壊という。α 崩壊は質量数（核子の個数）を変える崩壊様式, β 崩壊は核種を変える
が質量数（核子の個数）を変えない崩壊様式, γ 崩壊は核種（陽子の個数と中性子の個数の組み合わせ）を変えない崩壊様式である。この 3 つについて以下で述べる。

4.2.1

質量数を変える崩壊様式

重い原子核は α 線を出して崩壊するものが多い。この崩壊を α 崩壊という。α 線の正体は 42 He2 である。原子核 X の α 崩壊を反応式で表すと (4.1) 式のようになる。
4A−4A
Z XN →Z−2 XN−2 +2 He2(4.1)

この崩壊により質量数（核子の個数）が 4 だけ減少する。
陽子 2 個, 中性子 2 個をばらばらに放出するのではない理由は結合エネルギーにある。
4 個の核子がばらばらな状態を基準とした場合の 42 He2 の結合エネルギーは約 28 MeV[3]である。このエネルギーの分だけ陽子 2 個, 中性子 2 個をばらばらに放出するよりも残りの原子核が安定になるためである。

11

また, 中性子が過多な場合には, 中性子を 1 個放出する原子核も存在する。例として,
7
2 He5 はこの現象を起こす。さらに数は少なくなるが, 陽子を 1 個放出するものも存在する。こちらは 43 Li1 が例としてあげられる。なお, この 2 例は半減期がきわめて短い。
28

この他に α 粒子よりも重い Mg などの核を放出する崩壊様式や, 同程度の重さの 2個の核にちぎれる自発的核分裂という現象も存在し, それらも「質量数を変える崩壊様式」に含まれる。

4.2.2

核種を変えるが質量数を変えない崩壊様式

安定な原子核と比較して, 陽子や中性子のどちらかが多い場合, その原子核は陽子が中性子に, あるいは中性子が陽子に変化するとともに,β 線を放出するなどしてより安定な原子核へと崩壊する。この崩壊を β 崩壊という。
β 線の正体は電子または陽電子である。電子を放出する場合を β− 崩壊, 陽電子を放出する場合を β+ 崩壊という。
また, 核が核外にある軌道電子を捕獲することがある。この現象を軌道電子捕獲 (electron capture, EC と略す) という。原子核 X の β− 崩壊,β+ 崩壊, EC を反応式で表現するとそれぞれ（4.2）式, （4.3）式, （4.4）式のようになる。
AA−
Z XN →Z+1 XN−1 + e + ν¯e

(4.2)

A
A
+
Z XN →Z−1 XN+1 + e + νe

(4.3)

A
−
A
Z XN + e →Z−1 XN+1 + νe

(4.4)

上記のとおり,β 崩壊では質量数に変化は起こらないが原子番号が増減し, 同時に電子
ニュートリノ νe または反電子ニュートリノ ν¯e を放出する。

4.2.3

核種を変えない崩壊様式

4.2.1 や 4.2.2 で生じた原子核は多くが励起状態とよばれるエネルギーが高い状態にある。このためより低い状態に移行するために光子を放出する。これを γ 崩壊と言い,放出される光子は γ 線とよばれている。光子を放出する理由は原子核は基底状態に近いほうがより安定するからである。この光子のエネルギーはエネルギー準位差に相当するエネルギーである。
γ 崩壊は,α 崩壊や β 崩壊と違い, 核種（陽子の個数と中性子の個数の組み合わせ）を変えない崩壊様式である。また, γ 崩壊は電磁相互作用に起因する崩壊様式であるが, α崩壊は強い相互作用, β 崩壊は弱い相互作用に起因する崩壊様式である。
γ 線のエネルギーが低いとき, そして特に, γ 線のもつ角運動量が大きいときは, 原子核の周囲に存在する電子にγ線が吸収されて, γ 線の代わりにその電子が原子から跳ね飛ばされる現象が起きやすい。この跳ね飛ばされた電子を内部転換 (internal conversion)
12または内部電子転換と呼ぶ。ただし, この場合の γ 線は仮想光子であるにもかかわらず,
本来の γ 線が再吸収されるという誤解をもたれ易いため, 核子が電子を直接跳ね飛ばすという描像で説明されることも多い。
α 崩壊, β 崩壊の終状態が娘核の基底状態でない場合は, ほとんどの場合, 即座に基底状態に至るまで γ 崩壊を繰り返す。これを即発 γ 線と呼ぶ。しかし, 励起状態の中には,特異的に長寿命のものがあり, 核異性体 (nuclear isomer) と呼ばれる。核異性体から基底状態に至るまでの γ 崩壊は, 寿命の分だけ遅れて放射される。本研究の主題である放射性崩壊系列の計算では, 核異性体は基底状態は別の核のように扱う必要がある。

4.3

放射性崩壊系列

地球上には, 非常に多くの核種が存在している。まず, 安定核 (stable nucleus) と不安
定核 (unstable nucleus) に分けることができる。安定核とは自発的に崩壊を起こさない核種であり, 約 300 種類存在している。
不安定核は崩壊を起こす核種であり, 放射性核種 (radionuclide, radioactive nuclide) と呼ばれている。核種の存在・不存在の定義には任意性があるが, 核構造理論の立場では
約 1 万種が存在する（強い相互作用による自己束縛系としての構造計算の対象になりえる）と考えることが多い。放射性核種は, さらに人工放射性核種 (artificial radionuclide)と天然放射性核種 (natural radionuclide) に分けられる。
人工放射性核種は原子力発電や核爆発実験などで生成される, もともとは存在していなかった放射性核種のことである。天然放射性核種はさらに, 原始放射性核種 (primordialradionuclide) と宇宙線生成放射性核種 (cosmogenic nuclide) とに分けられる。
原始放射性物質は地球が作られたときから存在している放射性物質である。これらの原始放射性物質は非常に長い半減期を持っており, 18 種類 [5] 存在している。宇宙線
生成放射性核種は宇宙線と大気低層の気体が核反応を起こし, 崩壊したものであり, 約
20 種類 [5] 存在している。
放射性崩壊系列 (radioactive decay sequence) とは, この中の原始放射性物質が崩壊していく過程のことである。原始放射性物質には放射性崩壊系列に属するものと, 属さないものが存在しており, 例として前者にはウランやトリウム, 後者にはカリウム 40 がある。まずは α 崩壊系列について述べる。

4.3.1

α 崩壊系列α 崩壊系列としては 4 系列が存在している。そのうちの 3 系列は天然に存在するが 1
系列は存在していない。
天然に存在する 3 系列は, ウラン 238 を先祖とするウラン系列, トリウム 232 を先祖とするトリウム系列, ウラン 235 を先祖とするアクチニウム系列であり, それぞれ, 質量
数が整数 n を用いて 4n + 2, 4n, 4n + 3 と表される核種からなる崩壊系列である。次々に崩壊し最終的にそれぞれ質量数が 206, 208, 207 の鉛の安定な同位体になる。

13

天然に存在しない 1 系列は, ネプツニウム 237 のネプツニウム系列である。これは質量数が 4n + 1 の核種からなる系列であるが, ネプツニウム 237 の半減期が短いので, 地球ができてから現在までの間に全て崩壊して消えてしまい, 天然には存在しない。
これら 4 つの崩壊系列を表 4.1〜表 4.4 に記す。

14

表 4.1: ウラン系列 (4n + 2 系列)。 崩壊様式の欄で, SF は spontaneus fission(自発的核分裂), IT は isometric transision(異性体遷移, 即ち寿命の長い励起状態から基底状態への γ 遷移) を表す。なお複数の崩壊チャンネルが競合して崩壊が分岐している場合は,娘核を矢印 → の右側に示した。表は次ページにも継続している。
元素 質量数
陽子数 中性子数
半減期
崩壊様式
9
U
238
92
146
4.468 × 10 yαSF, 5.5 × 10−5 %
Th
234
90
144
24.10d
β−
Pa
234
91
143
6.70h
β−
234m
91
143
1.159mβ− ,99.84%
IT, 0.16%
5
U
234
92
142
2.455 × 10 yαSF, 1.6 × 10−9 %Mg, 1 × 10−11 %Ne, 9 × 10−12 %
Th
230
90
140
7.54 × 104 yα
24Ne, 6 × 10−11 %SF, 4 × 10−12 %
Ra
226
88
138
1.6 × 103 yα
14C, 3.2 × 10−9 %
Rn
222
86
136
3.8235d
α
Po
218
84
134
3.098m
α,99.98% →Pb
β− ,0.02% →At
At
218
85
133
1.5s
α,99.90% →Bi
β− ,0.10% →Rn
Rn
218
86
132
35ms
α
Pb
214
82
132
26.8m
β−
Bi
214
83
131
19.9m
β− ,99.98% →Po
α,0.02% →Tl
Po
214
84
130
164.3µs
α
Pb
212
82
130
10.64h
β−
Bi
212
83
129
60.55mβ− ,64.06%
α,35.94%
Po
212
84
128
0.299µs
α
Hg
210
80
130
300ns
β−
Tl
210
81
129
1.30m
β− →210 Pb
β− n,3.2 × 10−9 % →209 Pb
Pb
210
82
128
22.20yβ− →Biα,7.0 × 10−3 %→Hg
Bi
210
83
127
5.012dβ− →Poα,1.3×10−4 %→Tl

15

元素
Po
Pb
Bi
Tl
Pb
Hg
Tl
Pb

質量数
210
209
209
208
208
206
206
206

陽子数
84
82
83
81
82
80
81
82

4.3.2

β 崩壊系列

中性子数
126
127
126
127
126
126
125
124

半減期
138.376d
3.253h

3.053m8.32m
4.202m

崩壊様式
α
β−
安定
β−
安定
β−
β−
安定

核分裂反応で生まれた生成物 (fission product) は同じ質量数の安定核と比較して中性子が過多であるため, β− 崩壊を起こして中性子を陽子に変換していき, 安定核の陽子：中性子比率に近づこうとする。これが β 崩壊系列であり, 質量数毎に系列が存在しているため, 核分裂生成物に関するものに限っても数十もの系列が存在する。例として, 質
量数が 137 の β 崩壊系列を表 4.5 に示す。
文献 [1] では, β 崩壊系列の顕著な特徴として, 総崩壊率が時間に反比例して減少することが挙げられおり, また, その原因が１回崩壊する毎に半減期が大幅に増大することにあることが明確に説明されている。

16

元素
Th

質量数
232

Ra
Ac
Th

228
228
228

Ra

224

Rn
Po
Pb
Bi

220
216
212
212

Po
Pb

212
210

Bi

210

Po
Tl
Pb
Hg
Tl
Pb

210
208
208
206
206
206

表 4.2: トリウム系列 (4n 系列)。
陽子数 中性子数
半減期
崩壊様式
10
90
142
1.40 × 10 yαSF,1.1 × 10−9 %
88
140
5.75y
β−
89
139
6.15h
β−
90
138
1.9116yα20O,1.0 × 10−11 %88
136
3.6319dα
14C,4.0 × 10−9 %
86
134
55.6s
α
84
132
0.145sα8213010.64hβ−83
12960.55mβ− ,64.06 ％ →Po
α,35.94 ％ →Tl
84
128
0.229µs
α−82
12822.20yβ →Biα,1.9 × 10−6 %→Hg
83
127
5.012dβ− →Poα,1.3×10−4 %→Tl84
126
138.376d
α81
1273.053m
β−
82
126
安定
80
126
8.32m
β−
81
125
4.202m
β−
82
124
安定

17

元素
U

質量数
235

Th

231

Pa

231

Ra
Ac

227
227

Th
Fr

227
223

Ra

223

At

219

Rn
Pb
Bi
Po

219
215
215
215

At
Pb
Bi

215
211
211

Po
Pb
Bi
Hg
Tl
Pb

211
209
209
207
207
207

表 4.3: アクチニウム系列 (4n + 3 系列)。
陽子数 中性子数
半減期
崩壊様式
8
92
143
7.038 × 10 yαSF,7.0 × 10−9 %28Mg,8.0 × 10−10 %Ne,8.0 × 10−10 %90
14125.52hβ−α,4.0 × 10−1191
140
3.276 × 104 yαSF,2.0 × 10−11 %
88
139
42.2m
β−
89
138
21.772yβ− , 98.62 ％→Th
α,1.38 ％ →Fr
90
137
18.68dα−8713622.00mβ →Raα,(6.0 × 10−3 ％)→At
88
135
11.43dα14C,8.9 × 10−8 %
85
134
56s
α,97.00 ％ →Bi
β− ％→Rn
86
133
3.96s
α
82
133
147s
β−
83
132
7.6m
β−
84
131
1.781ms
α→Pb−
β ,(2.3 × 10−4 ％)→At
85
130
0.10ms
α
82
129
36.1m
β−
84
127
2.14m
α,99.72%→Tl
β− ,0.28%→Po84
127
0.516sα82
1273.253h
β−
83
126
安定
80
127
2.9m
β−
81
126
4.77m
β−
82
125
安定

18

元素
Np

質量数
237

Pa
U

233
233

Th
Ra
Ac

229
225
225

Fr

221

Ra

221

At

217

Rn
Bi

217
213

Po
Bi

213
211

Po
Tl
Pb
Bi
Tl
Pb
Hg
Tl

211
209
209
209
207
207
205
205

表 4.4: ネプツニウム系列 (4n + 1 系列)。
陽子数 中性子数
半減期
崩壊様式
6
93
144
2.144 × 10 yαSF,2.0 × 10−10 %91
142
26.975dβ−92
1411.592 × 105 yα24Ne,9 × 10−10 %SF,6 × 10−11 %28Mg,1 × 10−13 %
90
139
7392y
α
88
137
14.9d
β−
89
136
10.0d
α
14
C,4 × 10−12 %
87
134
4.9m
α→At
−
β ,0.10%→Ra
88
133
28s
β−
14
C,1 × 10−12 %8513232.3ms
α,99.99%→Bi−
β ,7.0 × 10−3 %→Bi86
131
0.54msα−8313045.59mβ ,97.80%→Po
α,2.20%→Tl84
129
3.72µs
α
83
128
2.14m
α,99.72%→Tl
β− ,0.28%→Po84
127
0.516sα811282.161mβ−82
1273.253h
β−
83
126
安定
81
126
4.77m
β−
82
125
安定
80
125
5.14m
β−
81
124
安定

19

元素
Sn
Sb
Te
I
Xe
Cs
Ba

質量数
137
137
137
137
137
137
137

表 4.5: 質量数が 137 の β 崩壊系列。
陽子数 中性子数 半減期 崩壊様式
50
87
190ms
β−
51
86
492ms
β−
52
85
2.49s
β−
53
84
24.5s
β−
54
83
3.818mβ−55
82
30.08y
β−
56
81

20

第5章

放射性崩壊の時間発展を記述す
る連立常微分方程式とその解本章では原子核の崩壊を表す方程式の一般解を求める。最も単純な場合から複雑な場合まで, 5.1 から 5.5 の 5 通りのケースについて, それぞれ一節を当てて説明する。この 5 通りの場合についての考察を一般化することで, 次章で述べる数十の核種からなるような崩壊系列にも適用できる汎用解法が自然と導かれる。この汎用解法については,
大沼・森近の卒業論文 [1] の 18〜19 頁の 3.2.2 節の数式 (3.10)〜(3.16) により完全に記述されているので, 本論文では説明を繰り返さない。本章執筆の目的は, 汎用解法の典型的な具体例を論じることで, 一般論だけを述べた大沼・森近の卒業論文 [1] を補うよ
うな説明を与えることである。

5.1

1 種類の原子核の崩壊を表す微分方程式とその一般解

放射性崩壊による原子核の個数の時間的変化について考える。N 個の原子核のうち単位時間に λN 個が崩壊する場合を考えると, N は時間 t に関する 1 階の常微分方程式

dN
= −λN
dt

(5.1)

に従う。(5.1) 式における λ は崩壊定数とよばれ, 1 個の原子核が単位時間当たりに崩壊する確率を表している。初期条件は,N(t0 ) = N0

(5.2)

ととる。(5.1) 式の形は変数分離型とよばれ, 容易に解を求めることができる。まず, (5.1)
式の両辺を N で割ると下式を得る。

1 dN
= −λ
N dt
次に (5.3) 式を時間 t に関して t = t0 から t = t1 まで定積分を行う。
∫ t1
∫ t1
1 dN
= −λ
dt
t0 N dt
t0

(5.3)

(5.4)

(5.4) 式における左辺は置換積分法を用い, 積分定数を t から N に変形する。したがって,
積分区間は t0 ≦ t ≦ t1 から N0 ≧ N ≧ N(t1 ) に写像される。したがって, (5.4) 式は以下の21

ようになる。

∫ N(t1 )N0

1
1
dN = −λ[t]t=t
t=t0
N

1)
⇔ [log N]N=N(t= −λ(t1 − t0 )N=N0⇔ log N(t1 ) − log N0 = −λ(t1 − t0 )⇔ log N(t1 ) = log N0 − λ(t1 − t0 )(5.5)

となる。ここで, (5.5) 式の両辺の exp をとると,

elog N(t1 ) = elog N0 −λ(t1 −t0 )⇔ N(t1 ) = elog N0 e−λ(t1 −t0 )⇔ N(t1 ) = N0 e−λ(t1 −t0 )(5.6)となる。次に, (5.6) 式において, t1 = t と置き換え, さらに式を簡略化するため t0 = 0 とすると解として,N(t) = N0 e−λt(5.7)
が得られる。(5.7) が 1 つの原子核の崩壊の時間変化を表す数式である。

5.2

1 種類の親核から 1 種類の娘核への崩壊

前節の例は, 崩壊の親核しか考慮していないので崩壊系列とは言えない。本節では,崩壊系列の第一の例として, 最も単純な崩壊系列である 1 種類の親核から 1 種類の娘核137−への崩壊を取り上げる。具体的には, 親核 137
55 Cs82 が娘核 56 Ba83 へと半減期 30.1 年で β崩壊する場合を考える。
137
137
55 Cs82 の個数を N1 (t) 個, 56 Ba83 の個数を N2 (t) 個とすると, それらの時間発展は下記の 2 本の微分方程式で記述される。

dN1= −λN1
dt

(5.8)

dN2
= λN1
(5.9)
dt
初期条件として, t0 = 0, N1 (0) = N0 , N2 (0) = 0 とする。
(5.8) 式は 5.1 節と同じように変数分離型であるので, 同じように解くことができ, その解は,
N1 (t) = N0 e−λt(5.10)となる。
次に, dtd (N1 + N2 ) = 0 より N1 + N2 =(一定) となり, 親核と娘核の個数の和は時間変化
せず一定値を保つことがわかる。その一定値は t = 0 での値 N0 + 0 = N0 と等しいので,

N2 (t) + N1 (t) = N022

(5.11)

となり, (5.10) を代入すると,

N2 (t) = N0 (1 − eλt )(5.12)

を得る。半減期の T 12 ＝ 30.1 年から崩壊定数を計算すると,

λ=ln 2 0.693≈= 7.30 × 10−10 s−1
T 1/2 30.1y(5.13)

である。
以上の手順を踏んで求まった解 (5.10) 式および (5.12) 式の時間経過をグラフで表したものが図 5.1 である。
N1
N2

1

0.8

N

0.6

0.4

0.2

0
0

10

20

30

40
50
time[year]

60

70

80

90

図 5.1: 137
55 Cs82 の崩壊の時間経過。赤色の実線は親核の個数（初期値を 1 に規格化）, 緑
色の破線は娘核の個数を表す。

5.3

親核の崩壊先が 2 つに分岐する場合

第 2 の例として, 40
19 K21 の場合を考える。この崩壊は枝分かれを起こし, 分岐比が 10.72−％の確率で軌道電子捕獲を起こし 40
18 Ar22 に崩壊し安定状態になり, 89.28 ％の確率で β91崩壊し 4020 Ca20 に崩壊し安定状態になる。半減期は T 2 = 1.277 × 10 年である。
4040
40
19 K21 の個数を N1 (t) 個, 18 Ar22 の個数を N2 (t) 個, 20 Ca20 の個数を N3 (t) 個とする。
このときの微分方程式は, それぞれ N1 (t), N2 (t), N3 (t) の時間発展を記述する。

dN1= −λN1dt
23(5.14)

dN2= λ1→2 N1(5.15)dt
dN3= λ1→3 N1(5.16)
dtとなる。初期条件として, t0 = 0, N1 (0) = N0 , N2 (0) = N3 (0) = 0 とする。
(5.14) 式は (5.1) 式と同様に変数分離型であるため, 同じように解くことができ, その解は,
N1 (t) = N0 e−λt(5.17)となる。また, N2 (t) は,∫ t

dN2 (t)dt + N2 (0)dt
0∫ t
dN1 (t)= λ1→2dt
dt
0

N2 (t) =となる。おなじく, N3 (t) は,(5.18)∫ t

dN3 (t)dt + N3 (0)dt
0∫ t
dN1 (t)= λ1→3dt
dt
0N3 (t) =(5.19)となる。
ここで, λ1→2 , λ1→3 は, 部分崩壊定数といい,と定義される。

λ1→2 = λ1 × 0.1072(5.20)

λ1→3 = λ1 × 0.8928(5.21)(5.18) 式, (5.19) 式により, N2 (t) および N3 (t) の割合は,N2 (t) : N3 (t) = λ1→2 : λ1→3

(5.22)となる。
ここで, dtd (N1 + N2 + N3 ) = 0 より, N1 + N2 + N3 = 一定であるから, 3 核種の個数の和は時間変化せず一定値 N0 を保つ。よって, それぞれの核種の個数の時間変化を表す解は,

N1 (t) = N0 eλt
λ1→2{N0 − N1 (t)}λ1
λ1→2 λt= N0e
λ1

N2 (t) =λ1→3{N0 − N1 (t)}λ1
λ1→3 λt= N0e
λ1

N3 (t) =(5.23)(5.24)(5.25)(5.26)(5.27)となる。
以上のようにして求まった解の (5.23) 式, (5.25) 式, (5.27) 式の時間変化をグラフに表したものが図 5.2 である。

24

N1
N2
N3

1

0.8

N

0.6

0.4

0.2

0
0

5000

10000

15000

20000
time[year]

25000

30000

35000

40
図 5.2: 40
19 K21 の崩壊の時間変化。赤色の実線は 19 K21 の個数（時刻 0 で 1 に規格化), 緑
40
色の破線は 40
18 Ar22 の個数, 青色の点線は 20 Ca20 の個数を表す。

5.4

親核が娘核, 娘核が孫核へと 2 段階にわたり崩壊する場
合第 3 の例として, 90
38 Sr52 を出発点とする崩壊系列を取り上げる。この崩壊は 2 段階に
90
90
−
1
わたり崩壊を起こし, まず 90
38 Sr52 が半減期 T 2 = 28.79 年で 39 Y51 に β 崩壊し, 次に 39 Y51−
が半減期 T 12 = 64.10 時間で 90
40 Zr50 に β 崩壊し安定核となる。
90
90
時刻 t における 90
38 Sr52 の個数を N1 (t) 個, 39 Y51 の個数を N2 (t) 個, 40 Zr50 の個数を N3 (t)個とする。初期値をそれぞれ N1 (0) = N0 , N2 (0) = 0, N3 (0) = 0 とする。また, 崩壊定数は 1 段階目を λ1 , 2 段階目を λ2 とすると,

λ1 =

ln 2
0.693
≈
= 7.269 × 10−10 s−1
T 1/2 28.79y(5.28)

λ2 =

ln 2
0.693
≈
= 3.004 × 10−6 s−1
T 1/2 64.10h(5.29)およびとなる。
N1 , N2 , N3 の従う連立 1 階常微分方程式は, λ1 , λ2 を使って以下のように表される。

 


 N1 (t)   −λ1 0 0   N1 (t) 
 


d 
 N2 (t)  =  λ1 −λ2 0   N2 (t) 
(5.30)
 
 

dt 
N3 (t)0
λ2 0
N3 (t)

25





 N1 (t)  −λ1 0 0 




あるいは,  N2 (t)  = ⃗n(t),  λ1 −λ2 0  = Λ とおくと, (5.30) 式は




N3 (t)
0λ2 0d⃗n(t) = Λ⃗n(t)dt(5.31)⃗n(t) = eΛt⃗n(0)(5.32)とも表すことができる。
これの解はとかける。eΛt を求めるには, まず Λ の固有値を求める。行列 Λ の固有値を ϵ, 固有ベク
トルを ⃗v とすると, 固有値方程式 Λ⃗v = ϵ⃗v が成り立つ。
行列 Λ の特性方程式 (または固有方程式ともいう) は固有値を ϵ, I を単位行列とするとdet(ϵI − Λ) = 0(5.33)

 ϵ + λ1
0
0 


det  −λ1 ϵ + λ2 0  = (ϵ + λ1 )(ϵ + λ2 )ϵ = 0


0
−λ2 0

(5.34)である。したがって,となるので, ϵ = −λ1 , −λ2 , 0 が固有値である。ϵ1 = −λ1 , ϵ2 = −λ2 , ϵ3 = 0 と書くことにする。
 vi1 
固有値 ϵi (1 ≤ i ≤ 3) に対応する固有ベクトル ⃗vi =  vi2  は Λ⃗vi = ϵ i v⃗i を満たすゼロ


vi3
ベクトルではないベクトルである。それらを求めた結果は



 

 0  0  1 


 
 −λ1 ⃗v1 =  λ1 −λ2  , ⃗v2 =  1  , ⃗v2 =  0 (5.35)  λ2 −11
λ1 −λ2である。なお, これらにゼロでない任意の定数を乗じたものもやはり固有ベクトルである。
行列 P を P = (v⃗1 v⃗2 v⃗3 ), 即ち Pi j = v ji と定義する。detP = 1 , 0 なので, 逆行列 P−1が存在する。
P−1 P = I = (e⃗1 , e⃗2 , e⃗3 )(5.36)    1  0  0    であるから (ただし e⃗1 =  0 , e⃗2 =  1 , e⃗3 =  0  とした。)
 
 
 
0
0
1

P−1 P = (P−1 v⃗1 , P−1 v⃗2 , P−1 v⃗3 ) = (e⃗1 , e⃗2 , e⃗3 )26(5.37)

となる。したがって, P−1 v⃗i = e⃗i (1 ≤ i ≤ 3) が成立する。これを用い相似変換を行うと,P−1 ΛP = P−1 (Λv⃗1 , Λv⃗2 , Λv⃗3 )= P−1 (ϵ 1 v⃗1 , ϵ 2 v⃗2 , ϵ 3 v⃗3 )= (ϵ 1 P−1 v⃗1 , ϵ 2 P−1 v⃗2 , ϵ 3 P−1 v⃗3 )= (ϵ 1 e⃗1 , ϵ 2 e⃗2 , ϵ 3 e⃗3 ) ϵ 1 0 0 =  0 ϵ 2 0 
0 0 ϵ3(5.38)となり, 対角行列であることが示せる。D = P−1 ΛP と書けば両辺に左から P, 右から P−1
を乗じて Λ = PDP−1 であるから, n ≤ 0 に対し,

Λn = (PDP−1 )n= PDP−1 PDP−1 · · · PDP−1= PD(P−1 P)D(P−1 P) · · · (P−1 P)DP−1= PDIDI · · · IDP−1= PDn P−1

(5.39)(5.40)および n ϵ 1
0
0 


Dn =  0 ϵ 2 n 0 
00 ϵ 3n(5.41)

が成立するので,
Λt

e

∞
∑
1

Λn tn
n!
n=0
 ∑∞ 1 n=0 n! (ϵ 1 t)n

0
0

 −1∑∞ 1n
 P= P 
00
n=0 n! (ϵ 2 t)

∑∞ 1

n 
0
0
n=0 n! (ϵ 3 t) ϵt e 1 0
0 
t


=  0 eϵ 2 0 

t 

0
0 eϵ 3=(5.42)(5.43)となる。

27

したがって, 解は ⃗n(t) = eΛt⃗n(0) であるから,








 N1 (t)  1
0 0   e−λ1 t
0 0   N0 
0 0   1  
 1
 N2 (t)  =  λ−λ−λ1
  0 e−λ2 t 0   λ λ−λ  0 
1
0
1
0

 

 1λ2 2   1 2N3 (t)
−1
1
0
0
1
0
1
−1
1
λ1 −λ2
N0 e−λ1 t


λ1
−λ2 t
−λ1 t

N0 λ1 −λ2 (e
−e )
= 
(5.44)
(
)


λ2
λ1
−λ1 t
−λ2 t
N0 λ1 −λ2 e− λ1 −λ2 e+1(5.45)
を得る。
以上の計算を経て得られた解の (5.44) 式の表す時間変化をグラフで表したものが図5.3 である。しかし, 図 5.3 では N2 (t) の値が小さすぎて読み取れない。その原因は娘核の半減期が, 親核の半減期や孫核の半減期（無限大）と比較して桁違いに短いからである。そこで, 両軸を対数目盛で表したグラフも図 5.4 として示しておく。
N1
N2
N3

1

0.8

N

0.6

0.4

0.2

0
0

10

20

30

40
50
time[year]

60

70

80

90

90
図 5.3: 90
38 Sr52 を起点とする崩壊系列の時間変化。赤色の実線は 38 Sr52 の個数（時刻 0 で90
1 に規格化), 緑色の破線は 90
39 Y51 の個数, 青色の点線は 40 Zr50 の個数を表す。

5.5

親核が複数の経路を経て同一の娘核に崩壊する場合

28

1

0.1

N1
N2
N3

0.01

N

0.001

0.0001

1e-0051e-0061e-0071e-008
0.001

0.01

0.1

1

10

time[year]

図 5.4: 90
38 Sr52 を起点とする崩壊系列の時間経過。図 5.3 を両対数スケールで表示しなおしたものである。
第 4 の例として, 第 1 の例としてとりあげた 137
55 Cs82 の場合を詳細に考えなおす。実は,
この崩壊は半減期 T 21 = 30.08 年で 5.4%の確率で β− 崩壊を起こし 137
56 Ba81 に崩壊し安定−137
核になるか, 94.6%の確率で β 崩壊を起こし 56 Ba81 の励起状態である 137m
56 Ba81 に崩壊
137
し, その後半減期 T 21 = 2.55 分で 56 Ba81 に崩壊し安定核になる 2 通りの経路を経る。
実際には多くの β 崩壊においてこのようなことが起きている。即ち,β 崩壊のほとん
どで, 親核の基底状態からの遷移が娘核の多くの励起状態に分岐して起こり, その直後に励起状態から基底状態への γ 崩壊が起きているのである。しかし,β 崩壊の半減期がγ崩壊の半減期よりはるかに長く, かつ, γ 崩壊の半減期程度の時間で起こる変動を無視してよい場合は, 親核が 100%直接に娘核の基底状態へ β 崩壊すると近似しても支障ないのである。本例はそのような近似的扱いの妥当性の数値的な確認としても役立つ。
137m
137
時刻 t における 137
55 Cs82 の個数を N1 (t) 個, 56 Ba81 の個数を N2 (t) 個, 56 Ba81 の個数をN3 (t) 個とする。初期値をそれぞれ N1 (0) = N0 , N2 (0) = 0, N3 (0) = 0 とする。また, 崩壊
定数を λ1 , λ2 , 同じく部分崩壊定数を λ1→2 , λ1→3 とすると, それぞれ

0.693
ln 2
≈
= 7.298 × 10−10 s−1
T 1/2 30.09y(5.46)

ln 2
0.693
≈
= 4.53 × 10−3 s−1
T 1/2 2.55m(5.47)

λ1→2 = λ1 × 0.946 = 6.904 × 10−10 s−1(5.48)

λ1→3 = λ1 × 0.054 = 3.94 × 10−11 s−1(5.49)

λ1 =

λ2 =

29

となる。
N1 (0) = N0 , N2 (0) = 0, N3 (0) = 0 の従う連立 1 階常微分方程式を考えると,



 
 N1 (t)   −λ1
0 0   N1 (t) 


 
d 
 N2 (t)  =  λ1→2 −λ2 0   N2 (t) 

 
 
dt 
N3 (t)
λ1→3 λ2 0
N3 (t)(5.50) N1 (t)  −λ10 0 

となる。 N2 (t)  = ⃗n(t),  λ1→2 −λ2 0  = Λ とおくと, (5.50) 式は




N3 (t)
λ1→3 λ2 0d⃗n(t) = Λ⃗n(t)dt(5.51)とも表すことができる。
行列 Λ は三角行列なので, その固有値は対角要素 λ1 , λ2 , 0 である。以下では固有値が縮退していない (即ちすべて値が異なる) と仮定して解を求める。
形式的な解は⃗n(t) = eΛt⃗n(0)(5.52)である。固有値に縮退がないので
Λ は相似変換で対角行列にできる。
 −λ1 0 0 
即ち, P−1 ΛP = D =  0 −λ2 0  を満たす行列 P が存在する。このとき,


0
0 0

 −λ t e 10 0 eΛt = PeDt P−1 = P  0 e−λ2 t 0  P−1


0
0 0

(5.53)となる。したがって, n⃗i (t)(1 ≤ i ≤ 3) は, e−λ1 t , e−λ2 t , 1 の線形結合となるので,

ni (t) = Ci1 e−λ1 t + Ci2 e−λ2 t + Ci3 (1 ≤ i ≤ 3)

(5.54)とおける。9 係数 Ci j (1 ≤ i ≤ 3, 1 ≤ j ≤ 3) は,

d⃗n = Λ⃗ndt(5.55)

から得られる 9 条件のうち, 詳しく考察すると独立なものが 6 条件あり, これと初期条
件の 3 条件とから定めることができる。
以上のようにして得られた解は

N1 (t) = N0 e−λ1 t
λ1→2 −λ2 t(e− e−λ1 t )N2 (t) = N0
λ1 − λ2
{
(
)
}
λ1→2
λ1→2 −λ2 t−λ1 t
N3 (t) = N0 1 − 1 −e−e
λ1 − λ2
λ1 − λ230(5.56)(5.57)(5.58)

である。
したがって, 以上の事柄を踏まえ (5.56) 式, (5.57) 式, (5.58) 式の時間経過をグラフで表したものが図 5.5 となる。しかし, 5.4 節と同じように図 5.5 では N3 (t) の値が小さす
ぎるため確認できない。したがって, 両軸を対数目盛で表したものを図 5.6 に示す。

N1
N2
N3

1

0.8

N

0.6

0.4

0.2

0
0

10

20

30

40
50
time[year]

60

70

80

90

137
図 5.5: 137
55 Cs82 からの崩壊の時間変化。赤色の実線は 55 Cs82 の個数（時刻 0 で 1 に規格化), 緑色の破線は 137m
56 Ba81 の個数（本グラフではほとんど水平軸に重なるため描かれていない）, 青色の点線は 137
56 Ba81 の個数を表す。

31

1

0.1

N1
N2
N3

0.01

N

0.001

0.0001

1e-0051e-0061e-0071e-0080.1

1

10time[year]

図 5.6: 137
55 Cs82 からの崩壊の時間変化。図 5.5 を両対数スケールで描き直したものである。

32

第6章
6.1

計算結果

数値計算プログラム作成の概略

前章で論じた 4 通りの具体例での解の求め方を一般化することで, 任意の崩壊系列の時間変化を表す連立常微分方程式の解法を作り上げることが可能であり, 実際に, その結果は, 大沼・森近の卒業論文 [1] の 18〜19 頁の 3.2.2 節の数式 (3.10)〜(3.16) に示されている。なお, 大沼・森近の卒業論文では, 半減期に偶然縮退がある場合には縮退をわずかに解くという近似法が使用されたが, 本研究の対象とする崩壊系列では偶然縮退はなかったので, この近似法は使用されていない。
大沼・森近はその一般解法を C 言語プログラムとして実装したが, 本論文では, それをプログラムの中核としてそのまま利用し, そこに新たに取り組もうとしている崩壊系
列上に存在する核種の半減期と崩壊の分岐比のデータを設定するサブルーチン (C 言語の関数) を付け加えて, 重い原子核のアルファ崩壊系列に適用できるようにした。付録には, そのようにして拡張した C 言語のプログラムの全ソースコードを公開した。
大沼・森近の卒業論文 [1] に掲載されたプログラムでは, 半減期・分岐比等のデータ
設定は, 質量数の値毎に存在する β 崩壊系列それぞれについて, (5.31) 式で定義された行列 Λ の行列要素に値を代入するサブルーチンを作り, それを呼び出すことで行われていた。本研究でのプログラム追加に当たっては, より本格的なプログラムに仕立て直すために, プログラムが入力データファイルを読み込むという形式に変更することも検討した。しかし, プログラミング作業のほとんどがデータファイルの入力作業になってしま
い, C 言語の使用経験を積みたいという卒業研究の副目的のひとつにあまり適さないという短所がともなった。これに対して, 大沼・森近の卒業論文のように, データを C 言
語プログラムのステートメントとして記述すれば, データの入力も, C 言語のプログラ
ミングとして行われることになるので, C 言語に触れている時間は遥かに長くなる。そこで, 本卒業研究でも, データを C 言語のステートメントとして記述する方法を踏襲することに決定した。
本研究でのデータ入力ステートメントはステートメント間のデータ整合性を, 部分的にだが, 自動的に保証するために C 言語のマクロ機能を使うことにした。このため, 大
森・森近の卒業論文と較べてやや抽象化している。しかし, それが C 言語のプログラムであることを意識することなくしては, 正しいデータ設定サブルーチンを作り上げることは決してできない。したがって, C 言語によるプログラミング経験を積むという副目
的は本卒業研究を通じてかなり果たすことができたという実感を抱いている。

33

6.2

α 崩壊系列の簡単な計算例

実際の α 崩壊系列について論じる前に, 簡単な例を用いて説明する。原子核 A,B,C,D,E
が,A→B→C→D→E の経路で崩壊する崩壊系列を考える。半減期は,T 21 (A) = 107 秒, T 12 (B) = 102 秒, T 12 (C) = 100 秒, T 12 (D) = 104 秒, T 12 (E) = ∞（安定核）と設定する。
これらの原子核の個数を NA ∼ NE とすると, 崩壊過程を表す微分方程式は下記のようになる。
dNA= −λA NA
(6.1)
dt
dNB
= λA NA − λB NB
(6.2)
dt
dNC
(6.3)
=
λB NB − λC NC
dt
dND=
λC NC − λD ND
(6.4)
dt
dNE
=
λD ND
(6.5)
dt
上式は定数係数連立 1 階線形常微分方程式なので一般解を求めることができる。
図 6.1 は NA ∼ NE の崩壊の時間変化を両対数スケールでプロットしたものである。但し, 時刻ゼロにおける起点の核の個数を１個に規格化した。図 6.2 は, 放射平衡とそれに至る前の過渡過程についての説明を図 6.1 に書き込んだものである。まず, 放射平衡が起きる前については, 黒色の矢印の部分が傾きが 1, オレンジ色の矢印の部分の傾きが 2,
水色の部分の矢印の傾きが 3 となっている。理由は以下の通りである。先程の方程式(6.2) において,λB NB が小さいうちは,

NB= λA NAdt(6.6)と考えることができ,λA NA を定数とみなすことができる。したがって,

NB = λA NA t

(6.7)

となり,t の 1 乗に比例することがわかる。同じように考えると,NC は λC NC が小さいう
ちは t の 2 乗に比例,ND は λD ND は t の 3 乗に比例していることがわかる。
次に, 紫色の矢印の区間は安定核である E の値を除いて値の変わらないプラトー領域となっており, このプラトー領域を含め以後で放射平衡 (secular equilibrium in radioactive
decay) が起きている。
放射平衡とは, 親核と娘核の量的な関係が時間的にほぼ一定の割合で推移することであり, 親核のほうが娘核より半減期が長いときに起こる状態である。
続いて, 放射能の量で表したものが図 6.3 になる。こちらでも, 放射能の量が同じ値で

34

A
B
C
D
E

T1/2(C)

T1/2(B)

T1/2(D)

T1/2(A)1

0.01

0.0001N

1e−0061e−0081e−0101e−0121e−014
1

100

10000
time[sec]

1e+0061e+0081e+010

図 6.1: 原子核 A∼E の時間変化。横軸は経過時間（秒）, 縦軸は核種の個数（時刻ゼロ
での A の個数を 1 に規格化した）を表す。

35

A
B
C
D
E

T1/2(C)

T1/2(B)

T1/2(D)

T1/2(A)1

0.01

Slope=1

0.00011e−006N

Slope=2
1e−0081e−010

Plateau region

1e−012

Slope=3

1e−014
1

100

10000
time[sec]

1e+0061e+0081e+010

図 6.2: 原子核 A∼E の時間変化。横軸は経過時間（秒）, 縦軸は核種の個数（時刻ゼロ
での A の個数を 1 に規格化した）を表す。

36

1e−007

A
B
C
D

1e−0081e−009D [Bq]

1e−0101e−0111e−0121e−0131e−0141e−015
1

100

10000
time[sec]

1e+0061e+0081e+010

図 6.3: 原子核 A∼E の時間変化。横軸は経過時間（秒）, 縦軸は放射能の量（ベクレル）
を表す。

37

減少しており, 放射平衡が起こっていることがわかる。
さて, このような長寿命核種を起点とする α 崩壊系列は, ダムが途中に存在する川の流れに例えることができる。図 6.4 のようにダム A∼D が存在する川があると考える。
ただし, 一本の川とし, 途中の合流はないものとし, ダム D の流出先は海（安定核）と考える。

図 6.4: 川およびダムを図示したもの。上流からダムを A∼D としている。
まず, 放射平衡が起こる前は, ダムに流入する水の量が流出する量よりも多いため, ダ
ムに水がたまっていく。プラトー領域のところでは, ダムに流入する水の量と流出する水の量が同じとなるため, 水位は変わらない。それ以後では流出する量が流入する量よりも多いため, ダム内の水は減っていく。但し, 減っていく割合はそれぞれ同じである。

6.3

超長寿命核種を起点とする放射平衡

ここから実際の α 崩壊系列について論じていく。自然に存在する α 崩解系列の源には, 系列上の他の核種より桁違いに長い半減期を持ったウランないしトリウム同位体が存在している。まず最初に, それらの長寿命核種を源とする崩壊系列の下流の時間変化を計算する。α 崩壊系列は質量数が 4 だけ異なる核種が連なって作られるのであるから, 4 章で説明した通り, トリウム系列 (A = 4n), ネプツニウム系列 (A = 4n + 1), ウラン
系列 (A = 4n + 2), アクチニウム系列 (A = 4n + 3) の４系列が存在する。
図 6.5 は半減期 140 億年の 232 Th を起点とするトリウム系列の崩壊の時間変化を両対数スケールでプロットしたものである。横軸は経過時間（秒）, 縦軸は崩壊系列上の核種の個数（時刻ゼロにおける起点の核の個数を１個に規格化した）を表している。
109 ∼1017 秒の間に, 安定核である終点の 208 Pb および核分裂生成物 (fission product, FP)

38

1

0.01

0.0001

1e-006

N

Th232
Ra228
Ac228
Th228
Ra224
Rn220
Po216
Pb212
Bi212
Po212
Pb210
Bi210
Po210
Tl208
Pb208
Hg206
Tl206
Pb206
FP

1e-0081e-0101e-012

1e-014
1

1000001e+010time[sec]

1e+0151e+020

図 6.5: α 崩壊のトリウム系列の時間変化。横軸は経過時間（秒）, 縦軸は核種の個数（時刻ゼロでの 232 Th の個数を 1 に規格化した）を表す。
以外の核種の個数はほぼ一定値を保っている。これは 6.4 節で説明したように放射平衡のプラトー領域を示している。
続いて図 6.6 は図 6.5 を放射能の量で両対数スケールでプロットしたものである。横軸は経過時間（秒）, 縦軸は崩壊系列における放射能の量を表している。こちらも 109
秒以降では同じ割合で減少していることがわかる。したがって、こちらでも放射平衡が起きていることが読み取れる。
図 6.7 は半減期 7 億年の 235 U を起点とするアクチニウム系列の崩壊の時間変化を両対数スケールでプロットしたものである。横軸は経過時間（秒）, 縦軸は崩壊系列上の核種の個数（時刻ゼロにおける起点の核の個数を１個に規格化した）を表している。
1012 ∼1015 秒の間に, 安定核である終点の 207 Pb,209 Bi, および核分裂生成物の総量は増加しているが, 他の核種の個数はほぼ一定値を保っている。これもやはり放射平衡のプラトー領域を示している。
続いて図 6.8 は図 6.7 を放射能の量で両対数スケールでプロットしたものである。横軸は経過時間（秒）, 縦軸は崩壊系列における放射能の量を表している。こちらも 1012秒以降では同じ割合で減少していることがわかる。したがって、こちらでも放射平衡が起きていることが読み取れる。
図 6.9 は半減期 45 億年の 238 U を起点とするウラン系列の崩壊の時間変化を両対数スケールでプロットしたものである。横軸は経過時間（秒）, 縦軸は崩壊系列上の核種の個数（時刻ゼロにおける起点の核の個数を１個に規格化した）を表している。1014 ∼1017
秒の間に, 安定核である終点の 209 Bi,208 Pb,206 Pb, および核分裂生成物の総量は増加しているが, 他の核種の個数はほぼ一定値を保っている。これもやはり放射平衡のプラトー

39

1e-016

1e-0181e-020

D[Bq]

Th232
Ra228
Ac228
Th228
Ra224
Rn220
Po216
Pb212
Bi212
Po212
Pb210
Bi210
Po210
Tl208
Hg206
Tl206
total

1e-022

1e-0241e-0261e-028
1

1000001e+010time[sec]

1e+0151e+020

図 6.6: α 崩壊のトリウム系列の時間変化。横軸は経過時間（秒）, 縦軸は崩壊系列における放射能の量を表している。全放射能量は total と表示されている。
1

0.01

0.0001

1e-006

N

U235
Th231
Pa231
Ra227
Ac227
Th227
Fr223
Ra223
At219
Rn219
Pb215
Bi215
Po215
At215
Pb211
Bi211
Po211
Pb209
Bi209
Hg207
Tl207
Pb207
FP

1e-0081e-0101e-012

1e-014
1

1000001e+010time[sec]

1e+0151e+020

図 6.7: アクチニウム崩壊系列の時間変化。横軸は経過時間（秒）, 縦軸は核種の個数（時刻ゼロでの 235 U の個数を 1 に規格化した）を表す。

40

1e-016

1e-0181e-020

D[Bq]

U235
Th231
Pa231
Ra227
Ac227
Th227
Fr223
Ra223
At219
Rn219
Pb215
Bi215
Po215
At215
Pb211
Bi211
Po211
Pb209
Hg207
Tl207
total

1e-022

1e-0241e-0261e-028
1

1000001e+010time[sec]

1e+0151e+020

図 6.8: アクチニウム崩壊系列の時間変化。横軸は経過時間（秒）, 縦軸は崩壊系列における放射能の量を表している。全放射能量は total と表示されている。

41

領域を示している。
1

0.01

0.0001

1e-006

N

U238
Th234
Pa234m
Pa234
U234
Th230
Ra226
Rn222
Po218
At218
Rn218
Pb214
Bi214
Po214
Pb212
Bi212
Po212
Hg210
Tl210
Pb210
Bi210
Po210
Pb209
Bi209
Tl208
Pb208
Hg206
Tl206
Pb206
FP

1e-0081e-0101e-012

1e-014
1

1000001e+010time[sec]

1e+0151e+020

図 6.9: ウラン崩壊系列の時間変化。横軸は経過時間（秒）, 縦軸は核種の個数（時刻ゼ
ロでの 238 U の個数を 1 に規格化した）を表す。
続いて図 6.10 は図 6.9 を放射能の量で両対数スケールでプロットしたものである。横軸は経過時間（秒）, 縦軸は崩壊系列における放射能の量を表している。こちらも 1014秒以降では同じ割合で減少していることがわかる。したがって、こちらでも放射平衡が起きていることが読み取れる。
図 6.11 は半減期 214 万年の 237 Np を起点とするネプツニウム系列の崩壊の時間変化を両対数スケールでプロットしたものである。横軸は経過時間（秒）, 縦軸は崩壊系列上の核種の個数（時刻ゼロにおける起点の核の個数を 1 個に規格化した）を表している。この場合にプラトー領域を示すのは, 106 〜 1013 秒の間の 233 Pa だけである。これはその直後の核種の 233 U が 16 万年という起点の核とさほど違わない半減期を持つためである。233 U を起点とすれば, その下流はプラトー領域を形成しうるが, 起点が 237 Npであるため, 1013 秒付近までは 233 U の個数が増え続けるので, 下流の核種の個数もそれに比例して増え続けるのである。しかし、1013 秒付近からは同じ割合で個数が減少しているので、放射平衡は起きている。
続いて図 6.12 は図 6.11 を放射能の量で両対数スケールでプロットしたものである。
横軸は経過時間（秒）, 縦軸は崩壊系列における放射能の量を表している。こちらも1013 秒以降では同じ割合で減少していることがわかる。したがって、こちらでも放射平衡が起きていることが読み取れる。

42

1e-016

1e-0181e-020

D[Bq]

U238
Th234
Pa234m
Pa234
U234
Th230
Ra226
Rn222
Po218
At218
Rn218
Pb214
Bi214
Po214
Pb212
Bi212
Po212
Hg210
Tl210
Pb210
Bi210
Po210
Pb209
Tl208
Hg206
Tl206
total

1e-022

1e-0241e-0261e-028
1

1000001e+010time[sec]

1e+0151e+020

図 6.10: ウラン崩壊系列の時間変化。横軸は経過時間（秒）, 縦軸は崩壊系列における放射能の量を表している。全放射能量は total と表示されている。

43

1

Np237
Pa233
U233
Th229
Ra225
Ac225
Fr221
Ra221
At217
Rn217
Bi213
Po213
Bi211
Po211
Tl209
Pb209
Bi209
Tl207
Pb207
Hg205
Tl205
FP

0.01

0.0001N

1e-0061e-0081e-0101e-012

1e-014
1

1000001e+010time[sec]

1e+0151e+020

図 6.11: ネプツニウム崩壊系列の時間変化。横軸は経過時間（秒）, 縦軸は核種の個数（時刻ゼロでの 237 Np の個数を 1 に規格化した）を表す。

44

1e−0141e−0161e−0181e−020

D[Bq]

Np237
Pa233
U233
Th229
Ra225
Ac225
Fr221
Ra221
At217
Rn217
Bi213
Po213
Bi211
Po211
Tl209
Pb209
Tl207
Hg205
total

1e−0221e−0241e−0261e−028
1

100

10000 1e+006 1e+008 1e+010 1e+012 1e+014 1e+016 1e+018
time[sec]

図 6.12: ネプツニウム崩壊系列の時間変化。横軸は経過時間（秒）, 縦軸は崩壊系列における放射能の量を表している。全放射能量は total と表示されている。

45

表 6.1: ウンウントリウムの崩壊系列。238 U 以降はウラン系列と同様であるため省略する。
元素 質量数
陽子数 中性子数
半減期
崩壊様式
Uut
278
113
165
0.24ms
α
Rg
274
111
163
6.4ms
α
Mt
270
109
161
2.0s
α
Bh
266
107
159
1.7s
α
Db
262
105
157
35s
α,67%
SF,33%
Lr
258
103
155
4.1s
α,95%
SF,5%
Md
254
101
153
28m
ϵ
Fm
254
100
154
3.240h
α,99.94%
SF,0.06%
Fm
250
98
152
13.08yα,99.92%
SF,0.08%
Fm
246
96
150
4706y
α,99.97%
SF,0.03%
Pu
242
94
148
3.240hαSF,5.5 × 10−4 %
U
238
92
146
4.468 × 109 yαSF, 5.5 × 10−5 %6.4

短寿命な超重核を起点とする α 崩壊系列

本研究では, 教科書にはのっていないような非常に重い短寿命な核種からの α 崩壊系
列についての計算を行ってみた。そのような崩壊系列は, 半減期や分岐のデータが完全には揃っておらず, 正確に崩壊の時間変化を計算できるものは少ない。我々は, データの揃った数少ない系列の中で, 2004 年に日本の理化学研究所で発見されたウンウント
リウム (Uut), 2002 年にロシアのドゥブナ合同原子核研究所で発見されたウンウンオク
チウム (Uuo) の崩壊について計算する。前節と同じように計算に必要なプログラムは付録に収録しておく。また,4.3.1 節と同じように,2 つの崩壊系列を表 6.1 及び表 6.2 に記す。
図 6.13 には, 278
113 Uut165 を起点とする崩壊系列の時間変化を両対数スケールでプロットしたものである。図からわかるとおり, この場合はほぼ放射平衡にはならないことがわ
かる。Uut はウラン系列の起点である 238 U へとつながっていくが, 238 U が桁違いに超寿命であるため, 計算をこの核で打ち切ってある。
続いて図 6.14 は図 6.13 を放射能の量で両対数スケールでプロットしたものである。
横軸は経過時間（秒）, 縦軸は崩壊系列における放射能の量を表している。こちらでもほぼ放射平衡は起きていないことが読み取れる。ここで, 減少していく部分で直線の傾きを取ると −1 になっている。これは, 半減期の短い核から長い核へと崩壊が進んでい

46

元素
Uuo
Lv
Fl

質量数
294
290
286

Cn

282

表 6.2: ウンウンオクチウムの崩壊系列。
陽子数 中性子数 半減期 崩壊様式
118
176
0.89ms
α
116
174
15ms
α
114
172
0.16s
α,40%
SF,60%
112
170
0.50msSF

1

Uut278
Rg274
Mt270
Bh266
Db262
Lr258
Md254
Fm254
Cf250
Cm246
Pu242
U238
FP

0.01

0.0001N

1e-0061e-0081e-0101e-012

1e-0141e-0101e-005

1

100000time[sec]

1e+0101e+015

図 6.13: 278 Uut を起点とする α 崩壊系列の時間変化。横軸は経過時間（秒）, 縦軸は核種の個数（時刻ゼロでの起点の核の個数を１個に規格化した）を表す。

47

くためで, これは大沼・森近の卒業論文 [1] で詳しく書いてある β 崩壊系列と同じことが起きている。なぜならば,278 Uut を起点とするそれぞれの核種の放射能の量をD = Ni λi(6.8)

とおくと, 式 (6.8) の値が一定になり,D は t の −1 乗に比例することになるためである。
また, 節で例えたようにダムがある川で例えると, 流入元のダムが空になり, 流入先のダムの水が溜りきった後で次のダムへと水が流出していく。このことの繰り返しであり,節の例えとは異なる川の流れになる。
1e+004

Uut278
Rg274
Mt270
Bh266
Db262
Lr258
Md254
Fm254
Cf250
Cm246
Pu242
total

1e+002
1e+0001e−002D[Bq]1e−0041e−0061e−0081e−0101e−0121e−0141e−005

1

100000time[sec]

1e+0101e+015

図 6.14: 278 Uut を起点とする α 崩壊系列の時間変化。横軸は経過時間（秒）, 縦軸は崩壊系列における放射能の量を表している。全放射能量は total と表示されている。
図 6.15 には, 294
118 Uuo176 を起点とする崩壊系列の時間変化を両対数スケールでプロットしたものである。図からわかるとおり, この場合も放射平衡にはならないことがわかる。この崩壊系列は全て核分裂へとつながって行く。核分裂生成物は, 大沼・森近の卒業論文で計算されたような β 崩壊系列で崩壊を続けることになる。
続いて図 6.16 は図 6.15 を放射能の量で両対数スケールでプロットしたものである。
横軸は経過時間（秒）, 縦軸は崩壊系列における放射能の量を表している。こちらでも放射平衡は起きていないことが読み取れる。こちらもデータは少ないが同じように、減少していく部分で直線の傾きを取ると −1 になっており、先程と同じようにベータ崩壊
系列と同じことが起きている。理由も同様である。

48

1

Uuo294
Lv290
Fl286
Cn282
FP

0.01

0.0001N

1e-0061e-0081e-0101e-012

1e-0141e-0101e-0081e-0060.0001time[sec]

0.011

図 6.15: 294 Uuo を起点とする α 崩壊系列の時間変化。横軸は経過時間（秒）, 縦軸は核種の個数（時刻ゼロでの起点の核の個数を１個に規格化した）を表す。

49

Uuo294
Lv290
Fl286
Cn282
total

1e+002

1e+000D[Bq]1e-002

1e-0041e-0061e-0081e-0101e-0101e-0081e-0060.0001time[sec]

0.011

図 6.16: 294 Uuo を起点とする α 崩壊系列の時間変化。横軸は経過時間（秒）, 縦軸は崩壊系列における放射能の量を表している。全放射能量は total と表示されている。

50

第7章

まとめ

本研究では原子核の α 崩壊系列について, その崩壊を決定する時間に関する一階の定数係数線形連立常微分方程式を解き, その時間変化を調べた。崩壊系列中の他の核種と比較して, 半減期が圧倒的に長い核を起点とする崩壊については, 放射平衡と呼ばれる現象が見られることを確認できた。プラトー領域を伴う放射平衡の例として, まずトリウム系列を含む 3 系列を示した。また, プラトー領域を伴わない放射平衡の例として, 同程度の半減期を持つ核種を 2 個含むネプツニウム系列を示した。次に, 半減期のごく短い超重核を起点とする崩壊を 2 例挙げ, 放射平衡を示さないことを見た。
本研究では遣り残した課題として, 崩壊系列の起点を他の核種に取った場合の計算も行ってみることで, 部分的な放射平衡が現れることがあるのかどうかをより詳細に論じるという課題が考えられる。以上, 今後本研究を引き継いで発展させてくれる卒業研究
生への課題としてここに記しておくものである。

51

参考文献
[1] 大沼正季, 森近厚平, 著,「235 U の全核分裂生成物の放射能の時間変化」, 福井大学工学部物理工学科卒業論文（2012）.
[2] Particle Physics Booklet, Particle Data Group (2012).
[3] J.K. Tuli, NUCLAR WALLET CARDS, http://www.nndc.bnl.gov/wallet/ (2011).
[4] B・ポッフ,K. リーツ,C. ショルツ,F. サッチャ, 著 柴田利明, 訳,「素粒子・原子核物理入門 改訂新版」, 丸善出版,（2012）.
[5] 藤家洋一,「原子力‐自然に学び, 自然を学ぶ‐」,ERC 出版,（2005）.
[6] 米国立核データセンター (National Nuclear Data Center), ウェブサイト
http://www.nndc.bnl.gov

52

url

謝辞
本論文を作成するにあたり, 田嶋直樹先生には終始丁寧なご指導をしていただいたことに感謝し, お礼申し上げます。また博士前期課程 1 年生の伊藤研人氏, 杉浦友章氏には研究上の様々な面でお世話をしていただいたことに深く感謝いたします。同じ研究室で卒業研究に取り組んだ杉本智彦氏, 團野良樹氏にも同様に感謝いたします。最後に,本研究に対してご意見をいただいた, 多くの物理工学科の先生方にもお礼申し上げ, 謝
辞の言葉とさせていただきます。

53

付録

1) 原子質量超過から原子核の束縛エネルギーを計算するプログラム#include <stdio.h>#include <math.h>main(){
int Z; //計算したい原子の原子番号int N; //中性子数
double Mx; //原子質量超過
double H = 938.7829; //(MeV):水素原子の質量
double n = 939.5652; //(MeV):中性子の質量
double u = 931.4939; //(MeV):原子質量単位
double a = 1.433e-5; //(MeV):原子中の電子の束縛エネルギーの計算式の係数
double E; //結合エネルギー
printf("Z を入力");
scanf("%d",&Z);printf("N を入力");
scanf("%d",&N);printf("M を入力");
scanf("%f",&M);printf("E=%f\n",M - Z *( H - u) - N *( n - u ) + a *( pow(Z,2.39) - Z );return 0;}2) 5.2 節で取り上げた放射性崩壊系列の解析解の値を計算するプログラム/*
原子核の放射性崩壊を表す
時間に関する１階の常微分方程式の解の値を表示するプログラム
2014 年 7 月 23 日 created from scratch by 松岡 弘和

54decay-ex1f に改良をくわえたもの
if 分岐を使い、年、日、時間、分、秒にそれぞれ対応。
*/#include <stdio.h>#include <math.h>#include <stdlib.h>#include <ctype.h>
int main(){double t; //時間
double halflife; //半減期
double lambda; //崩壊定数λdouble n1; //n1 原子核種１の個数（初期値＝１と規格化）
double n2; //n2 原子核種２の個数（初期値＝０と規格化）
double a; //作業変数 1
double t2; //作業変数 2
double halflife2; //作業変数 3(halflife で入力した年、日、時間、分を秒に変換)
char timeunit[16]; // 読み込んだ半減期の数値の時間の単位/*
printf("# λを入力 : ");scanf("%lf",&lambda);
printf("# lambda=%f\n",lambda);*/
printf("# T_(1/2)");
scanf("%lf",&halflife);
printf("# timeunit");scanf("%S",&timeunit);
printf("# halflife=%f unit=%s\n",halflife,timeunit);//
exit(0);if(toupper(timeunit[0])==’Y’){
halflife2=halflife * 365.25 * 86400;}
else if(toupper(timeunit[0])==’D’){
halflife2=halflife * 86400;}
else if(toupper(timeunit[0])==’H’){
halflife2=halflife * 3600;}
else if(toupper(timeunit[0])==’M’){
halflife2=halflife * 60;}
else if(toupper(timeunit[0])==’S’){
halflife2=halflife;}55

else{
printf("eroor\n");};printf("%f\n",halflife2);
t2=halflife2 / 500;
for(t=0;t<halflife2 * 3;t=t+t2){
a=-log(2) / halflife2 * t;n1=exp(a);n2=1 - n1;
printf("%e %f %f\n",t,n1,n2);}
return 0;}2’) プログラム 2 への入力データ
30.1
Y3) 5.3 節で取り上げた放射性崩壊系列の解析解の値を計算するプログラム/*
原子核の放射性崩壊を表す
時間に関する１階の常微分方程式の解の値を表示するプログラム
卒業論文の 5 章 3 節の例のような場合において使用できる。
2014 年 7 月 27 日 created by 松岡 弘和*/#include <stdio.h>#include <math.h>#include <stdlib.h>#include <ctype.h>

int main(){double t; //時間 (秒)double halflife; //半減期 (timeunit で指定される単位)
double n1; //n1 原子核種 1 の個数（初期値=1 と規格化）
double n2; //n2 原子核種 2 の個数（初期値=0 と規格化）

56

double n3; //n3 原子核種 3 の個数（初期値=0 と規格化）
double a; //作業変数 1
double b; //作業変数 2
double c; //作業変数 3
double p1; //分岐比 1
double p2; //分岐比 2
double halflife2; //半減期 (秒単位)
double t2; //halflife2 を 500 で割った値
char timeunit[16]; // halflife の時間の単位

printf("# T_(1/2)?:");
scanf("%lf",&halflife);
printf("# timeunit?:");
scanf("%s",timeunit);
printf("分岐比 1：");
scanf("%lf",&p1);
printf("分岐比 2：");
scanf("%lf",&p2);
printf("# halflife=%f

unit=%s\n",halflife,timeunit);

if(toupper((int)timeunit[0])==(int)’Y’){
halflife2=halflife * 365.25 * 86400;}
else if(toupper((int)timeunit[0])==(int)’D’){
halflife2=halflife * 86400;}
else if(toupper((int)timeunit[0])==(int)’H’){
halflife2=halflife * 3600;}
else if(toupper((int)timeunit[0])==(int)’M’){
halflife2=halflife * 60;}
else if(toupper((int)timeunit[0])==(int)’S’){
halflife2=halflife;}
else{
printf("eroor\n");};printf("# halflife=%f (sec)\n",halflife2);
printf("# time(sec), no. of nuclide 1, no. of nuclide 2,no. of nuclide 3\n");t2=halflife2 / 500;
for(t=0;t<halflife2 * 3;t=t+t2){
a=-log(2) / halflife2 * t;b=a * p1;
c=a * p2;n1=exp(a);57n2=b / a *(1-n1);n3=c / a *(1-n1);// 出力は、時刻（秒）、核 1 の個数（初期値=1）、核 2 の個数（初期値=0)、核 3 の個数 (初期値=0)
printf("%e %f %f %f\n",t,n1,n2,n3);}
return 0;}3’) プログラム 3 への入力データ
1.277e04
year
0.10720.89284) 5.4 節で取り上げた放射性崩壊系列の解析解の値を計算するプログラム/*
原子核の放射性崩壊を表す
例 3 のような場合で使用できる 2 段階で崩壊した場合のプログラム
2014 年 8 月 27 日 created from scratch by 松岡 弘和*/#include <stdio.h>#include <math.h>#include <stdlib.h>#include <ctype.h>double time_unit_converter(double halflife1,char *timeunit);
double time_unit_converter2(double halflife2,char *timeunit2);
int main(){double t; //時間
double halflife1; //半減期
double halflife2; //半減期
double n1; //n1 原子核種１の個数（初期値＝１と規格化）
double n2; //n2 原子核種２の個数（初期値＝０と規格化）
double n3; //n2 原子核種 3 の個数（初期値＝０と規格化）
double a; //作業変数 1
double b; //作業変数 2
double t2; //作業変数 3
double tuc; //半減期を表す関数 1
double tuc2; //半減期を表す関数 2
double hl1; //関数内における作業変数 1

58

double hl2; //関数内における作業変数 2char tu1[16]; //関数内における作業変数 3
char tu2[16]; //関数内における作業変数 4
printf("# T_(1/2)1");
scanf("%lf",&hl1);
printf("# timeunit");scanf("%s",tu1);
printf("# T_(1/2)2");
scanf("%lf",&hl2);
printf("# timeunit2");
scanf("%s",tu2);
printf("# halflife=%f unit=%s\n",hl1,tu1);
printf("# halflife=%f unit=%s\n",hl2,tu2);
tuc=time_unit_converter(hl1,tu1);
tuc2=time_unit_converter2(hl2,tu2);

printf("%f\n",tuc);for(t=3600;t<tuc * 3;t=t*1.03){
a=log(2) / tuc ;b=log(2) / tuc2 ;
n1=exp(-a * t);n2=a / (a - b) * (exp(-b * t) - exp(-a * t));n3=b / (a - b) * exp(-a * t) - (a / (a - b)) * exp(-b * t) + 1;
printf("%e %e %e %e\n",t,n1,n2,n3);}
return 0;}
double time_unit_converter(double halflife1,char *timeunit){
double halflife1a;
if(toupper((int)timeunit[0])==(int)’Y’){
halflife1a=halflife1 * 365.25 * 86400;}
else if(toupper((int)timeunit[0])==(int)’D’){
halflife1a=halflife1 * 86400;}
else if(toupper((int)timeunit[0])==(int)’H’){
halflife1a=halflife1 * 3600;}
else if(toupper((int)timeunit[0])==(int)’M’){
halflife1a=halflife1 * 60;59

}else if(toupper((int)timeunit[0])==(int)’S’){
halflife1a=halflife1;}
else{
printf("eroor\n");};return halflife1a;}double time_unit_converter2(double halflife2,char *timeunit2){
double halflife2a;if(toupper((int)timeunit2[0])==(int)’Y’){
halflife2a=halflife2 * 365.25 * 86400;}
else if(toupper((int)timeunit2[0])==(int)’D’){
halflife2a=halflife2 * 86400;}
else if(toupper((int)timeunit2[0])==(int)’H’){
halflife2a=halflife2 * 3600;}
else if(toupper((int)timeunit2[0])==(int)’M’){
halflife2a=halflife2 * 60;}
else if(toupper((int)timeunit2[0])==(int)’S’){
halflife2a=halflife2;}
else{
printf("eroor\n");};return halflife2a;}4’) プログラム 4 への入力データ
28.79
Y
64.10
H

5) 5.5 節で取り上げた放射性崩壊系列の解析解の値を計算するプログラム
60

/*
原子核の放射性崩壊を表す
2 段階で崩壊した場合のプログラム
2014 年 9 月 5 日 created from scratch by 松岡 弘和
if 分岐を使い、年、日、時間、分、秒にそれぞれ対応。
*/#include <stdio.h>#include <math.h>#include <stdlib.h>#include <ctype.h>double time_unit_converter(double halflife1,char *timeunit);
double time_unit_converter2(double halflife2,char *timeunit2);
int main(){double t; //時間
double halflife1; //半減期
double halflife2; //半減期
double n1; //n1 原子核種１の個数（初期値＝１と規格化）
double n2; //n2 原子核種２の個数（初期値＝０と規格化）
double n3; //n2 原子核種 3 の個数（初期値＝０と規格化）
double p1; //分岐比 1
double p2; //分岐比 2
double a; //作業変数 1
double b; //作業変数 2
double c; //作業変数 3
double d; //作業変数 4
double t2; //作業変数 5
double sum; //半減期を表す関数 1
double sum2; //半減期を表す関数 2
double hl1; //double hl2; //char tu1[16]; //char tu2[16]; //printf("# T_(1/2)1");
scanf("%lf",&hl1);
printf("# timeunit");scanf("%S",&tu1);
printf("p1");
scanf("%lf",&p1);
printf("# T_(1/2)2");
scanf("%lf",&hl2);
printf("# timeunit2");
scanf("%S",&tu2);
printf("p2");
scanf("%lf",&p2);

61

printf("# halflife=%f
printf("# halflife=%f

unit=%s
unit=%sp1=%f\n",hl1,tu1,p1);
p2=%f\n",hl2,tu2,p2);

sum=time_unit_converter(hl1,tu1);
sum2=time_unit_converter2(hl2,tu2);

printf("%f\n",sum);t2=sum / 500;
for(t=0;t<sum * 3;t=t+t2){a=log(2) / sum ;b=log(2) / sum2 ;c=a * p1;
d=b * p2;n1=exp(-a * t);n2=c / (a - b) * (exp(-b * t) - exp(-a * t));n3=1 - (1 - c / (a - b) ) * exp(-a * t) - c / (a - b) * exp(-b * t);printf("%e %f %f %f\n",t,n1,n2,n3);}
return 0;}
double time_unit_converter(double halflife1,char *timeunit){
double halflife1a;
if(toupper(timeunit[0])==’Y’){
halflife1a=halflife1 * 365.25 * 86400;}
else if(toupper(timeunit[0])==’D’){
halflife1a=halflife1 * 86400;}
else if(toupper(timeunit[0])==’H’){
halflife1a=halflife1 * 3600;}
else if(toupper(timeunit[0])==’M’){
halflife1a=halflife1 * 60;}
else if(toupper(timeunit[0])==’S’){
halflife1a=halflife1;}
else{
printf("eroor\n");};return halflife1a;}62

double time_unit_converter2(double halflife2,char *timeunit2){
double halflife2a;if(toupper(timeunit2[0])==’Y’){
halflife2a=halflife2 * 365.25 * 86400;}
else if(toupper(timeunit2[0])==’D’){
halflife2a=halflife2 * 86400;}
else if(toupper(timeunit2[0])==’H’){
halflife2a=halflife2 * 3600;}
else if(toupper(timeunit2[0])==’M’){
halflife2a=halflife2 * 60;}
else if(toupper(timeunit2[0])==’S’){
halflife2a=halflife2;}
else{
printf("eroor\n");};return halflife2a;}5’) プログラム 5 への入力データ/*
原子核の放射性崩壊を表す
2 段階で崩壊した場合のプログラム
2014 年 9 月 5 日 created from scratch by 松岡 弘和
if 分岐を使い、年、日、時間、分、秒にそれぞれ対応。
*/#include <stdio.h>#include <math.h>#include <stdlib.h>#include <ctype.h>double time_unit_converter(double halflife1,char *timeunit);
double time_unit_converter2(double halflife2,char *timeunit2);
int main(){double t; //時間

63

double halflife1; //半減期
double halflife2; //半減期
double n1; //n1 原子核種１の個数（初期値＝１と規格化）
double n2; //n2 原子核種２の個数（初期値＝０と規格化）
double n3; //n2 原子核種 3 の個数（初期値＝０と規格化）
double p1; //分岐比 1
double p2; //分岐比 2
double a; //作業変数 1
double b; //作業変数 2
double c; //作業変数 3
double d; //作業変数 4
double t2; //作業変数 5
double sum; //半減期を表す関数 1
double sum2; //半減期を表す関数 2
double hl1; //double hl2; //char tu1[16]; //char tu2[16]; //printf("# T_(1/2)1");
scanf("%lf",&hl1);
printf("# timeunit");scanf("%S",&tu1);
printf("p1");
scanf("%lf",&p1);
printf("# T_(1/2)2");
scanf("%lf",&hl2);
printf("# timeunit2");
scanf("%S",&tu2);
printf("p2");
scanf("%lf",&p2);
printf("# halflife=%f unit=%s
printf("# halflife=%f unit=%s

p1=%f\n",hl1,tu1,p1);
p2=%f\n",hl2,tu2,p2);

sum=time_unit_converter(hl1,tu1);
sum2=time_unit_converter2(hl2,tu2);

printf("%f\n",sum);t2=sum / 500;
for(t=0;t<sum * 3;t=t+t2){a=log(2) / sum ;b=log(2) / sum2 ;c=a * p1;
d=b * p2;n1=exp(-a * t);n2=c / (a - b) * (exp(-b * t) - exp(-a * t));n3=1 - (1 - c / (a - b) ) * exp(-a * t) - c / (a - b) * exp(-b * t);64

printf("%e %f %f %f\n",t,n1,n2,n3);}
return 0;}
double time_unit_converter(double halflife1,char *timeunit){
double halflife1a;
if(toupper(timeunit[0])==’Y’){
halflife1a=halflife1 * 365.25 * 86400;}
else if(toupper(timeunit[0])==’D’){
halflife1a=halflife1 * 86400;}
else if(toupper(timeunit[0])==’H’){
halflife1a=halflife1 * 3600;}
else if(toupper(timeunit[0])==’M’){
halflife1a=halflife1 * 60;}
else if(toupper(timeunit[0])==’S’){
halflife1a=halflife1;}
else{
printf("eroor\n");};return halflife1a;}double time_unit_converter2(double halflife2,char *timeunit2){
double halflife2a;if(toupper(timeunit2[0])==’Y’){
halflife2a=halflife2 * 365.25 * 86400;}
else if(toupper(timeunit2[0])==’D’){
halflife2a=halflife2 * 86400;}
else if(toupper(timeunit2[0])==’H’){
halflife2a=halflife2 * 3600;}
else if(toupper(timeunit2[0])==’M’){
halflife2a=halflife2 * 60;}65

else if(toupper(timeunit2[0])==’S’){
halflife2a=halflife2;}
else{
printf("eroor\n");};return halflife2a;}6) 一般の放射性崩壊系列を計算するためのファイル群

6-1) Makefile.PHONY: all clean reset
all: decay_chain_ac.exe
decay_chain_ac.exe: decay_chain_ac.o set_lambda.ogcc -o decay_chain_ac.exe decay_chain_ac.o set_lambda.o -lm
decay_chain_ac.o: decay_chain_ac.c
gcc -g -Wall -c decay_chain_ac.c
set_lambda.o: set_lambda.c
gcc -g -Wall -c set_lambda.cclean:
rm -f decay_chain_ac.o set_lambda.oreset:cleanrm -f decay_chain_ac.exe6-2) Header file/*
decay_chain_ac.h : to be included by decay_chain_ac.c and set_lambda.c
2014/9/5-10 created*/#define M 300 // maximum number of nuclides considered in the calculation#define BIGNUM 1.0e+300 // a number used as infinity#define BIGNUM0 1.0e+299 // a number used as maximum finite numberenum{sec,min,hour,day,year};#ifdef MAIN_SOURCE_FILE
double lambda[M][M]; // decay-rate matrix
double halflife[M]; // array of halflives
char *nuclide[M]; // array of pointers to character strings for nuclides
int nn; // the number of nuclides considered in the decay calculation66

#elseextern double lambda[M][M];extern char *nuclide[M];extern double halflife[M];extern int nn;#endifdouble time_in_sec(double time,int timeunit);
double halflife_to_decayrate(double halflife);6-3) The main source file of the C program/*
decay_chain_ac.c : decay chain of Actinides
2014/9/2-10 Version 1 created by utilizing decay_chain_fp.c (fp=fission products)*/#include <stdio.h>#include <math.h>#include <stdlib.h>#include <string.h> // for strncmp()#define MAIN_SOURCE_FILE#include "decay_chain_ac.h"double n[M]; // population of each nuclei
double n0[M]; // initial population of nuclei
double c[M][M]; // coefficients used to express the solutions
int alpha_chain(char *init_nucl, double tmin, double tmax, int ntime, FILE *fp_evo);
int set_initial_state(char *init_nucl);
int solve_c();double n_of_t(int i,double time);
int time_evolution(double tmin, double tmax, int ntime, FILE *fp_evo);
int make_script(FILE *fp);
int make_gnuplot_script(FILE *fp);
int check_sum_n(double tmin,double tmax,int ntime);
int check_derivative(double tmin,double tmax,int ntime);
int set_lambda();int main(){double tmin=1.0e-3; // 1ms ; -1.0 for automatic determination;
double tmax=1.0e+16; // 100years ; -1.0 for automatic determination;
int ntime=2000; // -1 for default value = 100;
FILE *fp_evo=NULL; // file to write time evolution
FILE *fp_scr=NULL; // file to write a shell script to divide time evolution data67

fp_evo=fopen("time_evolution.dat","wt");
if(fp_evo == NULL){
fprintf(stderr,"failed to open time evolution output file\n");exit(1);}//alpha_chain("Uuo294",tmin=1e-005,tmax=1e+40,ntime,fp_evo);//alpha_chain("Uut278",tmin=1e-06,tmax=1e+16,ntime,fp_evo);//alpha_chain("U238",tmin=1e-005,tmax=1e+40,ntime,fp_evo);//alpha_chain("Np237",tmin=1e-005,tmax=1e+40,ntime,fp_evo);//alpha_chain("U235",tmin=1e-005,tmax=1e+30,ntime,fp_evo);//alpha_chain("Th232",tmin=1e-005,tmax=1e+40,ntime,fp_evo);
if(fclose(fp_evo)!=0){
fprintf(stderr,"failed to close evolution output file\n");}fp_scr=fopen("time_evolution.sh","wt");
if(fp_scr == NULL){
fprintf(stderr,"failed to open shell-script to grep time_evolution.dat\n");exit(1);}
make_script(fp_scr);
if(fclose(fp_scr)!=0){
fprintf(stderr,"failed to close shell-script to grep time_evolution.dat\n");}fp_scr=fopen("time_evolution.plt","wt");
if(fp_scr == NULL){
fprintf(stderr,"failed to open gnuplot script to plot time_evolution data\n");exit(1);}
make_gnuplot_script(fp_scr);
if(fclose(fp_scr)!=0){
fprintf(stderr,"failed to close gnuplot script to plot time_evolution data\n");}return 0;}
int alpha_chain(char *init_nucl, double tmin, double tmax, int ntime, FILE *fp_evo){
int i,j,k,l,ng,degeneracy;
set_lambda();set_initial_state(init_nucl);
printf("initial population = 1 for %s\n",init_nucl);
for(j=0;j<nn;j++){printf(" For nuclide %s, n(t=0)=%6.3f\n",nuclide[j],n0[j]);for(i=0;i<j;i++){ // check whether there is the same nuclei
if(strncmp(nuclide[i],nuclide[j],16)==0){
printf("Error: Entries no. %d and no. %d are the same\n",i,j);

68

}
}
}

printf("\n lambda[0..%d][0..%d]=\n",nn-1,nn-1);ng=0; // increased from zero if matrix lambda is not lower triangulardegeneracy=0; // changed to >0 if diagonal elements are equalfor(k=0;k<nn;k++){for(l=0;l<nn;l++){//printf(" %15.7e",lambda[k][l]);printf(" %s",lambda[k][l] == 0.0 ? "." : "*" );if(l>k && lambda[k][l]!=0.0) {ng++;
fprintf(stderr,"\n error non zero lambda for %d %s %d %s\n",k,nuclide[k],l,nuclide[if(l!=k &&fabs(lambda[k][k]-lambda[l][l])/(fabs(lambda[k][k])+fabs(lambda[l][l])+1.0e-100)<1.0e-13) degeneracy++;}
printf("\n");}if(ng>0){
printf("lambda is not lower triangular %d\n",ng);exit(1);}if(degeneracy>0){
printf("eigen values of lambda are degenerated %d/2\n",degeneracy);exit(1);}
printf("\n n0[0..%d]=\n",nn-1);
for(k=0;k<nn;k++){printf(" %12.5e",n0[k]);}printf("\n");solve_c();
printf("\n c[0..%d][0..%d]=\n",nn-1,nn-1);
for(k=0;k<nn;k++){for(l=0;l<nn;l++){printf(" %15.7e",c[k][l]);}printf("\n");}if(tmin < 0.0){
tmin=BIGNUM;
for(i=0;i<nn;i++){if(halflife[i]>BIGNUM0) continue;69if(tmin>halflife[i]) tmin=halflife[i];}tmin*=0.01;}
if(tmax < 0.0){
tmax=0.0;
for(i=0;i<nn;i++){if(halflife[i]<BIGNUM0) tmax+=halflife[i];}tmax*=30.0;}
if(tmin>tmax || tmin <= 0.0 || tmax > BIGNUM0){
fprintf(stderr,"error: tmin=%e tmax=%e\n",tmin,tmax);exit(0);}if(ntime<0) ntime=100;
time_evolution(tmin,tmax,ntime,fp_evo);
check_sum_n(tmin,tmax,ntime);
check_derivative(tmin,tmax,ntime);
return 0;}
int set_initial_state(char *init_nucl){
int i,c=0;
for(i=0;i<nn;i++){if(strncmp(init_nucl,nuclide[i],16)==0){n0[i]=1.0;
c++;
}
else{
n0[i]=0.0;}}
if(c==1){
return 0;}
else{
fprintf(stderr,"set_initial_state: error: %d %s\n",c,init_nucl);
for(i=0;i<nn;i++) fprintf(stderr,"%d %s\n",i,nuclide[i]);exit(1);}}
int solve_c(){
int i,j,k;double cjk; //c_[j][k]double laikckj; //lambda_[i][k]*c_[k][j]70

for(j=0;j<nn;j++){if(j==0){c[j][j]=n0[0];}else{
cjk=0;
for(k=0;k<=j-1;k++){cjk+=c[j][k];}c[j][j]=n0[j]-cjk;}
if(j<nn-1){for(i=j+1;i<nn;i++){laikckj=0;
for(k=j;k<=i-1;k++){laikckj+=(lambda[i][k]*c[k][j]);}c[i][j]=(laikckj)/(-lambda[i][i]+lambda[j][j]);}}}return 0;}
double n_of_t(int i,double time){
int j;
double s=0;
for(j=0;j<=i;j++){s+=c[i][j]*exp(lambda[j][j]*time);}
return s;}
int time_evolution(double tmin, double tmax, int ntime, FILE *fp_evo){
int i,it;double r,rr,t,bq;
r=pow(tmax/tmin,1.0/(ntime-1));
printf("\n# time n[0..%d]\n",nn-1);
rr=1;
for(it=0; it<=ntime; it++){if(it==0){
t=0;
}
else{
t=tmin*rr;
rr*=r;}
for(i=0;i<nn;i++){71

n[i]=n_of_t(i,t);bq=n[i]*fabs(lambda[i][i]); // becquerel= decays per second
if(bq<1.0e-99) bq=0.0;
if(fabs(n[i])<1.0e-99) n[i]=0.0;
if(fp_evo != NULL){
fprintf(fp_evo,"%11.5e %12.5e %12.5e %d %s \n",t,n[i],bq,it,nuclide[i]);}}}return 0;}
int make_script(FILE *fp){
int i;
if(fp==NULL) return 0;for(i=0;i<nn;i++){fprintf(fp,"grep %s time_evolution.dat > %s.dat\n",nuclide[i],nuclide[i]);}return 0;}
int make_gnuplot_script(FILE *fp){
int i;
if(fp==NULL) return 0;
fprintf(fp,"set logscale xy\n i=1\n j=2\n plot [1e-2:1e+32][1e-16:1.5]\\\n");for(i=0;i<nn-1;i++){fprintf(fp,"\"%s.dat\" u i:j t \"%s\" w l,\\\n",nuclide[i],nuclide[i]);}fprintf(fp,"\"%s.dat\" u i:j t \"%s\" w l\n",nuclide[nn-1],nuclide[nn-1]);
return 0;}
int check_sum_n(double tmin,double tmax,int ntime){
double s=0.0,s0=0.0,ds=0.0,maxerror=0.0,time=0.0;
double r,rr;
int i,itime=0;
for(i=0;i<=nn-1;i++){s0+=n0[i];}r=pow(tmax/tmin,1.0/(ntime-1));
rr=1;
for(itime=1;itime<=ntime;itime++){time=tmin*rr;
rr*=r;
s=0.0;
for(i=0;i<nn;i++){n[i]=n_of_t(i,time);72

if(n[i]>=0){s+=n[i];}else if(n[i]<-1.0e-16){ // -1.0e-20 makes false err message
printf("check_sum_n:error:negative n[%d]=%15.7e at t=%d\n",i,n[i],itime);}}
ds=fabs(s-s0);
if(ds>1.0e-9){
printf("check_sum_n:warning:large error %.5e at time %.5e\n",ds,time);}
if(maxerror<ds) maxerror=ds;}
printf("\nSigma n’s max error =%.7e, max relative error=%.7e\n",maxerror,maxerror/(s0+1.0e-99));
return 0;}
int check_derivative(double tmin,double tmax,int ntime){
double sumn,dn,ds,dt,time,dtopt,r,rr;
double relative_error, max_relative_error=0.0, maxerror=0.0;
double ds1,ds2,ds3,ds4,ds5,dsmin;
int i,j,k,itime;
const int printlevel = 0;
r=pow(tmax/tmin,1.0/(ntime-1));
rr=1;
for(itime=0;itime<=ntime;itime++){ //ntime_partition
if(itime==0){time=0;}else{
time=tmin*rr;
rr*=r;}
if(printlevel>0) printf("%11.5e",time);
for(i=0;i<nn;i++){n[i]=n_of_t(i,time);
sumn=0.0;
for(j=0;j<nn;j++){sumn+=lambda[i][j]*n[j];}dsmin=ds1=ds2=ds3=ds4=ds5=BIGNUM;
dt=1.0e-12;
for(k=1;k<=200;k++){ // search dt which minimizes the error
dn=(- n_of_t(i,time+dt*2)+8*n_of_t(i,time+dt)-8*n_of_t(i,time-dt) + n_of_t(i,time-dt*2))/(12*dt);
ds5=ds4; ds4=ds3; ds3=ds2; ds2=ds1; ds1=fabs(dn-sumn);
ds=(ds1+ds2+ds3+ds4+ds5)/5.0;73

if(dsmin > ds) {dsmin=ds; dtopt=dt;}
dt=dt*1.77827941; // 10ˆ(1/4)=1.778
if(dt>tmax) break;}ds=dsmin;if(maxerror<ds) maxerror=ds;
relative_error=ds/(fabs(sumn)+1.0e-50);
if(printlevel==1) printf(" %12.5e",ds);if(printlevel==2) printf(" %12.5e",relative_error);if(printlevel==3) printf(" %12.5e %12.5e",ds,dtopt);
if(printlevel>3) printf(" %12.5e %12.5e",relative_error,dtopt);
if(max_relative_error<relative_error) max_relative_error=relative_error;}
if(printlevel>0)printf("\n");}printf("\ndn/dt’s max error=%.7e, max relative error=%.7e\n",maxerror,max_relative_error);
return 0;}//-------------------------------------------------------------------double time_in_sec(double time,int timeunit){
switch(timeunit){case sec :return time;break;case min :return time*60.0;break;case hour :return time*3600.0;break;case day :return time*86400.0;break;case year :return time*86400.0*365.0;break;default :fprintf(stderr,"time_in_sec:unknown time unit\n");exit(1);}}
double halflife_to_decayrate(double halflife){ // T_{1/2} -> lambda converter
const double ln2=log(2.0);
static int zero_counter=0;
if(halflife>BIGNUM0){
return (++zero_counter +10)*1.0e-31; // to avoid degeneracy of zero}
return (ln2/(halflife));}//--------------------------------------------------------------------

6-4) Data routine for the example chain/*74set_lambda.c : set decay rate values for decay chains of ex
2015/1/30 created*/#include <stdio.h>#include <stdlib.h>// for exit()#include "decay_chain_ac.h"#define DIAGELEM(NC,HL,TU) {nuclide[NC]=#NC;\
halflife[NC]=time_in_sec(HL,TU);\
lambda[NC][NC]=-halflife_to_decayrate(halflife[NC]);}// NC:nuclide, HL:halflife, TU:time unit
#define OFFDIAG(P,D,BR); lambda[D][P]=-lambda[P][P]*(BR)// P:parent nuclide, D:daughter nuclide, BR:branching ratio
int set_lambda(){enum{A,
B,
C,
D,
E,
END};
int i,j;nn=END; // the number of nuclides, nn is a global variableif(M<nn) exit(1);for(i=0;i<nn;i++){for(j=0;j<nn;j++)lambda[i][j]=0;}// diagonal elements of decay-rate matrix
DIAGELEM(A, 10e7
,sec );
DIAGELEM(B, 10e2
,sec );
DIAGELEM(C, 1.0,sec );
DIAGELEM(D, 10e4
,sec );
DIAGELEM(E, BIGNUM,sec );// off-diagonal elements of decay-rate matrix
OFFDIAG(A, B, 1.0 );OFFDIAG(B, C, 1.0 );OFFDIAG(C, D, 1.0 );OFFDIAG(D, E, 1.0 );return 0;}

6-5) Data routine for the Thorium chain/*
set_lambda.c : set decay rate values for decay chains of Actinides
2014/9/8-10 created75*/#include <stdio.h>#include <stdlib.h>// for exit()#include "decay_chain_ac.h"#define DIAGELEM(NC,HL,TU) {nuclide[NC]=#NC;\
halflife[NC]=time_in_sec(HL,TU);\
lambda[NC][NC]=-halflife_to_decayrate(halflife[NC]);}// NC:nuclide, HL:halflife, TU:time unit
#define OFFDIAG(P,D,BR) lambda[D][P]=-lambda[P][P]*(BR)// P:parent nuclide, D:daughter nuclide, BR:branching ratio
int set_lambda(){enum{Th232,
Ra228,
Ac228,
Th228,
Ra224,
Rn220,
Po216,
Pb212,
Bi212,
Po212,
Pb210,
Bi210,
Po210,
Tl208,
Pb208,
Hg206,
Tl206,
Pb206,FP,
END};
int i,j;double r; // branching ratio
nn=END; // the number of nuclides, nn is a global variableif(M<nn) exit(1);for(i=0;i<nn;i++){for(j=0;j<nn;j++)lambda[i][j]=0;}// diagonal elements of decay-rate matrix//
nuclides in A=4n sequence
DIAGELEM(Th232, 1.40e10,year );
DIAGELEM(Ra228, 5.75
,year );
DIAGELEM(Ac228, 6.15
,hour );
DIAGELEM(Th228, 1.9116,year );
DIAGELEM(Ra224, 3.6319,day );
DIAGELEM(Rn220, 55.6
,sec );

76

DIAGELEM(Po216, 0.145,sec );
DIAGELEM(Pb212, 10.64,hour);
DIAGELEM(Bi212, 60.55,min );
DIAGELEM(Po212, 0.299e-6 ,sec );DIAGELEM(Pb210, 22.20 ,year );DIAGELEM(Bi210, 5.012 ,day );DIAGELEM(Po210, 138.376 ,day );DIAGELEM(Tl208, 3.053,min );
DIAGELEM(Pb208, BIGNUM,sec );
DIAGELEM(Hg206, 8.32
,min );
DIAGELEM(Tl206, 4.202,min );
DIAGELEM(Pb206, BIGNUM,min );
DIAGELEM(FP,BIGNUM, min );//nuclides in A=4n+1 sequence//

nuclides in A=4n+2 sequence//

nuclides in A=4n+3 sequence// off-diagonal elements of decay-rate matrix//
decays from nuclides in A=4n sequence
r=1.1e-11;
OFFDIAG(Th232, Ra228, 1-r );OFFDIAG(Th232, FP,r );OFFDIAG(Ra228, Ac228, 1.0 );OFFDIAG(Ac228, Th228, 1.0 );
r=1e-13;
OFFDIAG(Th228, Ra224, 1-r );OFFDIAG(Th228, Pb208, r );r=4.0e-11;
OFFDIAG(Ra224, Rn220, 1-r );OFFDIAG(Ra224, Pb210, r );OFFDIAG(Rn220, Po216, 1.0 );OFFDIAG(Po216, Pb212, 1.0 );OFFDIAG(Pb212, Bi212, 1.0 );r=0.3594e-2;
OFFDIAG(Bi212, Po212, 1-r );OFFDIAG(Bi212, Tl208, r);
OFFDIAG(Po212, Tl208, 1.0 );r=1.9e-8;
OFFDIAG(Pb210, Bi210, 1-r );OFFDIAG(Pb210, Hg206, r );r=1.3e-6;
OFFDIAG(Bi210, Po210, 1-r );OFFDIAG(Bi210, Tl206, r );OFFDIAG(Tl208, Pb208, 1.0 );OFFDIAG(Hg206, Tl206, 1.0 );OFFDIAG(Tl206, Pb206, 1.0 );//decays from nuclides in A=4n+1 sequence

77

//

decays from nuclides in A=4n+2 sequence//

decays from nuclides in A=4n+3 sequence
return 0;}

6-6) Data routine for the Neptunium chain/*
set_lambda.c : set decay rate values for decay chains of Actinides
2014/9/8-10 created*/#include <stdio.h>#include <stdlib.h>// for exit()#include "decay_chain_ac.h"#define DIAGELEM(NC,HL,TU) {nuclide[NC]=#NC;\
halflife[NC]=time_in_sec(HL,TU);\
lambda[NC][NC]=-halflife_to_decayrate(halflife[NC]);}// NC:nuclide, HL:halflife, TU:time unit
#define OFFDIAG(P,D,BR); lambda[D][P]=-lambda[P][P]*(BR)// P:parent nuclide, D:daughter nuclide, BR:branching ratio
int set_lambda(){enum{Np237,
Pa233,
U233,
Th229,
Ra225,
Ac225,
Fr221,
Ra221,
At217,
Rn217,
Bi213,
Po213,
Bi211,
Po211,
Tl209,
Pb209,
Bi209,
Tl207,
Pb207,

78

Hg205,
Tl205,FP,
END};
int i,j;double r,r2,r3; // branching ratio
nn=END; // the number of nuclides, nn is a global variableif(M<nn) exit(1);for(i=0;i<nn;i++){for(j=0;j<nn;j++)lambda[i][j]=0;}// diagonal elements of decay-rate matrix//
nuclides in A=4n sequence//

nuclides in A=4n+1 sequence
DIAGELEM(Np237, 2.144e6,year );
DIAGELEM(Pa233, 26.975,day );
DIAGELEM(U233, 1.592e5,year );
DIAGELEM(Th229, 7932
,year );
DIAGELEM(Ra225, 14.9
,day );
DIAGELEM(Ac225, 10.0
,day );
DIAGELEM(Fr221, 4.9,min );
DIAGELEM(Ra221, 28,sec );
DIAGELEM(At217, 32.3e-3,sec );
DIAGELEM(Rn217, 0.54e-3,sec );
DIAGELEM(Bi213, 45.59,min );
DIAGELEM(Po213, 3.72e-6,sec );
DIAGELEM(Bi211, 2.14
,min );
DIAGELEM(Po211, 0.516,sec );
DIAGELEM(Tl209, 2.161,min );
DIAGELEM(Pb209, 3.253,hour );
DIAGELEM(Bi209, BIGNUM,sec );
DIAGELEM(Tl207, 4.77
,min );
DIAGELEM(Pb207, BIGNUM,min );
DIAGELEM(Hg205, 5.14
,min );
DIAGELEM(Tl205, BIGNUM,min );
DIAGELEM(FP,BIGNUM,min );//
nuclides in A=4n+2 sequence//

nuclides in A=4n+3 sequence// off-diagonal elements of decay-rate matrix//
decays from nuclides in A=4n sequence//

decays from nuclides in A=4n+1 sequence
r=2.0e-12;
OFFDIAG(Np237, Pa233, 1-r );OFFDIAG(Np237, FP,r );OFFDIAG(Pa233, U233, 1.0 );
r=9e-12;r2=6e-13;79r3=1e-15;
OFFDIAG( U233, Th229, 1-r-r2-r3 );OFFDIAG( U233, Pb209, r );OFFDIAG( U233, FP,r2 );OFFDIAG( U233, Hg205, r3 );OFFDIAG(Th229, Ra225, 1.0 );OFFDIAG(Ra225, Ac225, 1.0 );
r=4e-14;
OFFDIAG(Ac225, Fr221, 1-r );OFFDIAG(Ac225, Bi211, r );r=0.1e-2;
OFFDIAG(Fr221, At217, 1-r );OFFDIAG(Fr221, Ra221, r );r=1e-14;
OFFDIAG(Ra221, Rn217, 1-r );OFFDIAG(Ra221, Pb207, r );r=7.0e-5;
OFFDIAG(At217, Bi213, 1-r );OFFDIAG(At217, Rn217, r );OFFDIAG(Rn217, Po213, 1.0 );r=2.20e-2;
OFFDIAG(Bi213, Po213, 1-r );OFFDIAG(Bi213, Tl209, r);
OFFDIAG(Po213, Pb209, 1.0 );r=0.28e-2;
OFFDIAG(Bi211, Tl207, 1-r );OFFDIAG(Bi211, Po211, r );OFFDIAG(Po211, Pb207, 1.0 );OFFDIAG(Tl209, Pb209, 1.0 );OFFDIAG(Pb209, Bi209, 1.0 );OFFDIAG(Tl207, Pb207, 1.0 );OFFDIAG(Hg205, Tl205, 1.0 );//decays from nuclides in A=4n+2 sequence//

decays from nuclides in A=4n+3 sequence
return 0;}

6-7) Data routine for the Uranium chain/*
set_lambda.c : set decay rate values for decay chains of Actinides
2014/9/8-10 created*/#include <stdio.h>#include <stdlib.h>// for exit()

80

#include "decay_chain_ac.h"#define DIAGELEM(NC,HL,TU) {nuclide[NC]=#NC;\
halflife[NC]=time_in_sec(HL,TU);\
lambda[NC][NC]=-halflife_to_decayrate(halflife[NC]);}// NC:nuclide, HL:halflife, TU:time unit
#define OFFDIAG(P,D,BR); lambda[D][P]=-lambda[P][P]*(BR)// P:parent nuclide, D:daughter nuclide, BR:branching ratio
int set_lambda(){enum{U238,
Th234,
Pa234m,Pa234,
U234,
Th230,
Ra226,
Rn222,
Po218,
At218,
Rn218,
Pb214,
Bi214,
Po214,
Pb212,
Bi212,
Po212,
Hg210,
Tl210,
Pb210,
Bi210,
Po210,
Pb209,
Bi209,
Tl208,
Pb208,
Hg206,
Tl206,
Pb206,FP,
END};
int i,j;double r,r2,r3; // branching ratio
nn=END; // the number of nuclides, nn is a global variableif(M<nn) exit(1);for(i=0;i<nn;i++){for(j=0;j<nn;j++)lambda[i][j]=0;}// diagonal elements of decay-rate matrix

81

//

nuclides in A=4n sequence//

nuclides in A=4n+1 sequence//

nuclides in A=4n+2 sequence
DIAGELEM(U238, 4.468e9,year );
DIAGELEM(Th234, 24.10,day );
DIAGELEM(Pa234m, 1.159,min );
DIAGELEM(Pa234, 6.70
,hour );
DIAGELEM(U234, 2.455e5,year );
DIAGELEM(Th230, 7.54e4,year );
DIAGELEM(Ra226, 1600
,year );
DIAGELEM(Rn222, 3.8235,day );
DIAGELEM(Po218, 3.098,min );
DIAGELEM(At218, 1.5,sec );
DIAGELEM(Rn218, 35e-3,sec );
DIAGELEM(Pb214, 26.8
,min );
DIAGELEM(Bi214, 19.9
,min );
DIAGELEM(Po214, 164.3e-6,min );
DIAGELEM(Pb212, 10.64,hour);
DIAGELEM(Bi212, 60.55,min);
DIAGELEM(Po212, 0.299e-6,sec);
DIAGELEM(Hg210, 300e-9,sec );
DIAGELEM(Tl210, 1.30
,min );
DIAGELEM(Pb210, 22.20,year );
DIAGELEM(Bi210, 5.012,day );
DIAGELEM(Po210, 138.376,day );
DIAGELEM(Pb209, 3.253,hour );
DIAGELEM(Bi209, BIGNUM,sec );
DIAGELEM(Tl208, 3.053,min);
DIAGELEM(Pb208, BIGNUM,sec );
DIAGELEM(Hg206, 8.32
,min );
DIAGELEM(Tl206, 4.202,min );
DIAGELEM(Pb206, BIGNUM,sec );
DIAGELEM(FP,BIGNUM,sec)//
nuclides in A=4n+3 sequence// off-diagonal elements of decay-rate matrix//
decays from nuclides in A=4n sequence//

decays from nuclides in A=4n+1 sequence//

decays from nuclides in A=4n+2 sequence
r=5.5e-7
OFFDIAG(U238, Th234, 1-r );OFFDIAG(U238, FP, r );OFFDIAG(Th234, Pa234m, 1.0 );r=0.16e-2
OFFDIAG(Pa234m, U234, 1-r );OFFDIAG(Pa234m, Pa234, r );82

OFFDIAG(Pa234, U234, 1.0 );r=1.6e-11;
r2=1e-13;
r3=9e-14;
OFFDIAG(U234, Th230, 1-r-r2-r3 );OFFDIAG(U234, FP, r );//SF
OFFDIAG(U234, Hg210, r2 );//Mg
OFFDIAG(U234, Pb214, r3 );//Ne
r=6e-13;
r2=4e-14;
OFFDIAG(Th230, Ra226, 1-r-r2 );OFFDIAG(Th230, Hg206, r);//24Ne
OFFDIAG(Th230, FP , r2 );
r=3.2e-11;
OFFDIAG(Ra226, Rn222, 1-r );OFFDIAG(Ra226, Pb212, r ); //14C
OFFDIAG(Rn222, Po218, 1.0 );r=0.02e-2;
OFFDIAG(Po218, Pb214, 1-r );OFFDIAG(Po218, At218, r );r=0.10e-2;
OFFDIAG(At218, Bi214, 1-r );OFFDIAG(At218, Rn218, r );OFFDIAG(Rn218, Po214, 1.0 );OFFDIAG(Pb214, Bi214, 1.0 );r=0.02e-2;
OFFDIAG(Bi214, Po214, 1-r );OFFDIAG(Bi214, Tl210, r );OFFDIAG(Po214, Pb210, 1.0 );OFFDIAG(Pb212, Bi212, 1.0 );r=35.94e-2
OFFDIAG(Bi212, Po212, 1-r );OFFDIAG(Bi212, Tl208, r );OFFDIAG(Po212, Pb208, 1.0 );OFFDIAG(Hg210, Tl210, 1.0 );r=7.0e-5;
OFFDIAG(Tl210, Pb210, 1-r );OFFDIAG(Tl210, Pb209, r ); //β nr=1.9e-8;
OFFDIAG(Pb210, Bi210, 1-r );OFFDIAG(Pb210, Hg206, r );r=1.3e-6;
OFFDIAG(Bi210, Po210, 1-r );OFFDIAG(Bi210, Tl206, r );OFFDIAG(Po210, Pb206, r );OFFDIAG(Pb209, Bi209, 1.0 );OFFDIAG(Tl208, Pb208, 1.0 );OFFDIAG(Hg206, Tl206, 1.0 );OFFDIAG(Tl206, Pb206, 1.0 );//decays from nuclides in A=4n+3 sequence83

return 0;}6-8) Data routine for the Actinide chain/*
set_lambda.c : set decay rate values for decay chains of Actinides
2014/9/8-10 created*/#include <stdio.h>#include <stdlib.h>// for exit()#include "decay_chain_ac.h"#define DIAGELEM(NC,HL,TU) {nuclide[NC]=#NC;\
halflife[NC]=time_in_sec(HL,TU);\
lambda[NC][NC]=-halflife_to_decayrate(halflife[NC]);}// NC:nuclide, HL:halflife, TU:time unit
#define OFFDIAG(P,D,BR); lambda[D][P]=-lambda[P][P]*(BR)// P:parent nuclide, D:daughter nuclide, BR:branching ratio
int set_lambda(){enum{U235,
Th231,
Pa231,
Ra227,
Ac227,
Th227,
Fr223,
Ra223,
At219,
Rn219,
Pb215,
Bi215, //Po215, //At215,
Pb211,
Bi211,
Po211,
Pb209,
Bi209,
Hg207,
Tl207,
Pb207,
FP,84

END};
int i,j;double r,r2,r3; // branching ratio
nn=END; // the number of nuclides, nn is a global variableif(M<nn) exit(1);for(i=0;i<nn;i++){for(j=0;j<nn;j++)lambda[i][j]=0;}// diagonal elements of decay-rate matrix//
nuclides in A=4n sequence//

nuclides in A=4n+1 sequence//

nuclides in A=4n+2 sequence//

nuclides in A=4n+3 sequence
DIAGELEM(U235, 7.038e8,year );
DIAGELEM(Th231, 25.52,hour );
DIAGELEM(Pa231, 3.276e4,year );
DIAGELEM(Ra227, 42.2
,min );
DIAGELEM(Ac227, 21.772,year );
DIAGELEM(Th227, 18.68,day );
DIAGELEM(Fr223, 22.00,min );
DIAGELEM(Ra223, 11.43,day );
DIAGELEM(At219, 56,sec );
DIAGELEM(Rn219, 3.96
,sec );
DIAGELEM(Pb215, 147,sec );
DIAGELEM(Bi215, 7.6,min );
DIAGELEM(Po215, 1.781e-3,sec );
DIAGELEM(At215, 0.10e-3,sec );
DIAGELEM(Pb211, 36.1
,min );
DIAGELEM(Bi211, 2.14
,min );
DIAGELEM(Po211, 0.516,sec );
DIAGELEM(Pb209, 3.253,hour );
DIAGELEM(Bi209, BIGNUM,sec );
DIAGELEM(Hg207, 2.9,min );
DIAGELEM(Tl207, 4.77
,min );
DIAGELEM(Pb207, BIGNUM,sec );
DIAGELEM(FP, BIGNUM,sec );// off-diagonal elements of decay-rate matrix//
decays from nuclides in A=4n sequence//

decays from nuclides in A=4n+1 sequence//

decays from nuclides in A=4n+2 sequence//

decays from nuclides in A=4n+3 sequence
r=7.0e-11;
r2=8.0e-12;
r3=8.0e-12;
OFFDIAG(U235, Th231, 1-r-r2-r3 );85

OFFDIAG(U235, FP, r );OFFDIAG(U235, Hg207, r2 );OFFDIAG(U235, Pb215, r3 );r=4.0e-13;
OFFDIAG(Th231, Pa231, 1-r );OFFDIAG(Th231, Ra227, r );r=2e-13;
OFFDIAG(Pa231, Ac227, 1-r );OFFDIAG(Pa231, FP, r );OFFDIAG(Ra227, Ac227, 1.0 );r=1.38e-2;
OFFDIAG(Ac227, Th227, 1-r );OFFDIAG(Ac227, Fr223, r );OFFDIAG(Th227, Ra223, 1.0 );r=6.0e-5;
OFFDIAG(Fr223, Ra223, 1-r );OFFDIAG(Fr223, At219, r );r=8.9e-10;
OFFDIAG(Ra223, Po215, 1-r );OFFDIAG(Ra223, Pb209, r );r=3.00e-2;
OFFDIAG(At219, Bi215, 1-r );OFFDIAG(At219, Rn219, r );OFFDIAG(Rn219, Po215, 1.0 );OFFDIAG(Pb215, Bi215, 1.0 );OFFDIAG(Bi215, Po215, 1.0 );r=2.3e-6;
OFFDIAG(Po215, Pb211, 1-r );OFFDIAG(Po215, At215, r );OFFDIAG(At215, Bi211, 1.0 );OFFDIAG(Pb211, Bi211, 1.0 );r=0.28e-2;
OFFDIAG(Bi211, Tl207, 1-r );OFFDIAG(Bi211, Po211, r );OFFDIAG(Po211, Pb207, 1.0 );OFFDIAG(Pb209, Bi209, 1.0 );OFFDIAG(Hg207, Tl207, 1.0 );OFFDIAG(Tl207, Pb207, 1.0 );return 0;}

6-9) Data routine for the chain from Uut/*
set_lambda.c : set decay rate values for decay chains of Actinides
2014/9/8-10 created*/86

#include <stdio.h>#include <stdlib.h>// for exit()#include "decay_chain_ac.h"#define DIAGELEM(NC,HL,TU) {nuclide[NC]=#NC;\
halflife[NC]=time_in_sec(HL,TU);\
lambda[NC][NC]=-halflife_to_decayrate(halflife[NC]);}// NC:nuclide, HL:halflife, TU:time unit
#define OFFDIAG(P,D,BR) lambda[D][P]=-lambda[P][P]*(BR)// P:parent nuclide, D:daughter nuclide, BR:branching ratio
int set_lambda(){enum{Uut278,
Rg274,
Mt270,
Bh266,
Db262,
Lr258,
Md254,
Fm254,
Cf250,
Cm246,
Pu242,
U238,
FP,
END};
int i,j;double r; // branching ratio
nn=END; // the number of nuclides, nn is a global variableif(M<nn) exit(1);for(i=0;i<nn;i++){for(j=0;j<nn;j++)lambda[i][j]=0;}// diagonal elements of decay-rate matrix
DIAGELEM(Uut278, 0.24e-3,sec );
DIAGELEM(Rg274, 6.4e-3,sec );
DIAGELEM(Mt270, 2.0e-3,sec );
DIAGELEM(Bh266, 1.7,sec );
DIAGELEM(Db262, 35,sec );
DIAGELEM(Lr258, 4.1,sec );
DIAGELEM(Md254, 28,min);
DIAGELEM(Fm254, 3.240,hour);
DIAGELEM(Cf250, 13.08,year);
DIAGELEM(Cm246, 4706
,year);
DIAGELEM(Pu242, 3.75e5,year);
DIAGELEM(U238, BIGNUM,sec);
DIAGELEM(FP, BIGNUM,sec );

87

// off-diagonal elements of decay-rate matrix
OFFDIAG(Uut278, Rg274, 1.0 );OFFDIAG(Rg274, Mt270, 1.0 );OFFDIAG(Mt270, Bh266, 1.0 );OFFDIAG(Bh266, Db262, 1.0 );r=0.33;
OFFDIAG(Db262, Lr258, 1-r );OFFDIAG(Db262, FP,r);
r=0.05;
OFFDIAG(Lr258, Md254, 1-r );OFFDIAG(Lr258, FP, r );OFFDIAG(Md254, Fm254, 1.0 );r=0.06e-2;
OFFDIAG(Fm254, Cf250, 1-r );OFFDIAG(Fm254, FP, r );r=0.08e-2;
OFFDIAG(Cf250, Cm246, 1-r );OFFDIAG(Cf250, FP, r );r=0.03e-2;
OFFDIAG(Cm246, Pu242, 1-r );OFFDIAG(Cm246, FP, r );r=5.5e-6;
OFFDIAG(Pu242, U238, 1-r );OFFDIAG(Pu242, FP, r );return 0;}6-10) Data routine for the chain from Uuo/*
set_lambda.c : set decay rate values for decay chains of Actinides
2014/9/8-10 created*/#include <stdio.h>#include <stdlib.h>// for exit()#include "decay_chain_ac.h"#define DIAGELEM(NC,HL,TU) {nuclide[NC]=#NC;\
halflife[NC]=time_in_sec(HL,TU);\
lambda[NC][NC]=-halflife_to_decayrate(halflife[NC]);}// NC:nuclide, HL:halflife, TU:time unit
#define OFFDIAG(P,D,BR); lambda[D][P]=-lambda[P][P]*(BR)// P:parent nuclide, D:daughter nuclide, BR:branching ratio
int set_lambda(){88

enum{Uuo294,
Lv290,
Fl286,
Cn282,FP,
END};
int i,j;double r; // branching ratio
nn=END; // the number of nuclides, nn is a global variableif(M<nn) exit(1);for(i=0;i<nn;i++){for(j=0;j<nn;j++)lambda[i][j]=0;}// diagonal elements of decay-rate matrix
DIAGELEM(Uuo294, 0.89e-3,sec );
DIAGELEM(Lv290, 15e-3,sec );
DIAGELEM(Fl286, 0.16
,sec );
DIAGELEM(Cn282, 0.50e-3,sec );
DIAGELEM(FP, BIGNUM,sec );// off-diagonal elements of decay-rate matrix
OFFDIAG(Uuo294, Lv290, 1.0 );OFFDIAG(Lv290, Fl286,
1.0 );
r=0.4
OFFDIAG(Fl286, FP, 1-r );OFFDIAG(Fl286, Cn282, r);
OFFDIAG(Cn282, FP,
1.0);
return 0;}7) 放射性崩壊系列の総放射能量をベクレル単位で計算するプログラム/*
総放射線量を計算するプログラム。
2014/10/29 created Matsuoka Hirokazu*/#include <stdio.h>#include <math.h>#include <stdlib.h>#include <string.h> // for strncmp()
int main(){
FILE *fp; // file to write time evolution
FILE *fp2;
char buffer[256];89

double t,n,d;
int i;
int i0=-1;
double dsum=0.0;
double t0=-1.0;
fp=fopen("time_evolution.dat","rt");
fp2=fopen("totBq.dat","wt");
while(fgets(buffer,256,fp)!=NULL){//printf("%s",buffer);sscanf(buffer,"%le %le %le %d",&t,&n,&d,&i);if(i!=i0){
if(i0>=0)printf("%e %e :totbq\n",t0,dsum);if(i0>=0)fprintf(fp2,"%e %e\n",t0,dsum);
i0=i;
t0=t;
dsum=0.0;}
dsum+=d;//printf("%e %e %e %d\n",t,n,d,i);}
printf("%e %e :totbq\n",t0,dsum);
fprintf(fp2,"%e %e\n",t0,dsum);close(fp);close(fp2);}

90