確率論的解釈

　ミクロな粒子を用いた二重スリット実験で、干渉縞が現れるということは、

• １つの粒子が一方のスリットを通ったのと同時に

• もう一方のスリットを通ったことを意味する

このミクロな粒子は何らかの空間的な広がりを持つことを示している。

このような実験結果を解釈するアイデアとして、以下のような量子力学のコペンハーゲン解釈が生まれました。

• 観測前には波動関数が空間的な広がりをもちシュレディンガー方程式に従う

• 観測時点では一点に収束している

• 検出確率が波動関数の二乗に比例する

　ここでは確率をつかって出来事を解釈する問題について取り上げていきます。

3 枚のコインの裏表

　３つのコインを投げて、表 (H:heads) と裏 (T:tails) の枚数を数えてみます。

Coin 1	Coin 2	Coin 3
H	H	H
H	H	T
H	T	H
T	H	H
H	T	T
T	T	H
T	H	T
T	H	T

　すべてのコインは見分けがつかないこととします。このとき H の個数を考えます。

　上記の表から数を数えると、

• 3つHがある場合は、HHH

• 2つHがある場合は、HHT, HTH, THH

• 1つHがある場合は、HTT, TTH, THT

• Hがない場合は、TTT

となっていることがわかります。

　ここで特定の数の H 状態をまとめてマクロ状態（macrostate）、またそれぞれの状態をミクロ状態（microstate）と呼びます。マクロ状態は一般的に観測することが容易なもの、ミクロ状態は観測が難しいものに対応させて考えるとよいでしょう。

　各ミクロ状態は等確率で存在しますが、マクロ状態は存在確率が異なります。

　上の例を見ると、3つが H のマクロ状態よりも、2つが H のマクロ状態の方が存在確率が高いことがわかります。

　ここで 多重度 \(\Omega\) (multiplicity)という量を導入します。多重度とは、特定のマクロ状態におけるミクロ状態の数です。

　上のコインゲームで、 n 回表が出た場合の数として \(\Omega(n)\) を定義します。

　コインの総数を N とすると、多重度について次の式を導くことができます。

\[\Omega(N,n) == \frac{N!}{n!\cdot(N-n)!} = \binom{N}{n}\]

　このような多重度を計算するプログラムを作成してください。正しい実行例は以下のようになります。

type a value for N: 3
when n= 0 ;  the multiplicity is 1.0
when n= 1 ;  the multiplicity is 3.0
when n= 2 ;  the multiplicity is 3.0
when n= 3 ;  the multiplicity is 1.0

N を増やすとどうなるでしょうか。5、10、20 などでも N を試してみましょう。

多重度 \(\Omega\) (multiplicity)の関数を定義する

import matplotlib.ticker as mtick
import matplotlib.pyplot as plt


#Multiplicity calculation
def factor(N):
    if N == 0 or N == 1:
        return 1
    else:
        total = 1
        for i in range(2,N+1):
            total *= i
        return total

def Multiplicity(N,n):
    return factor(N)/factor(n)/factor(N-n)

入力されたNに対して多重度を計算する

N=int(input("type a positive integer value for N: "))
N_series = range(N+1)
Omega = []

for n in N_series:
    print("when n=", n, ";  the multiplicity is", Multiplicity(N,n))

type a positive integer value for N: 3
when n= 0 ;  the multiplicity is 1.0
when n= 1 ;  the multiplicity is 3.0
when n= 2 ;  the multiplicity is 3.0
when n= 3 ;  the multiplicity is 1.0

確率をグラフにして考える

　次にNによる違いを詳しくみていきましょう。以下では曲線をプロットするための関数Plot_Omegaを定義しています。

def Plot_Omega(N):

    N_series = range(N+1)
    Omega = []

    for n in N_series:
        Omega.append(Multiplicity(N,n))

    plt.plot(N_series, Omega)
    plt.title('N='+str(N))
    plt.xlabel('n')
    plt.ylabel('$\Omega(n)$')

    plt.show()

Plot_Omega関数

これで任意の N の結果をより便利に呼び出すことができます。

Plot_Omega(50)

数を増やしましょう

Plot_Omega(100)

Nが大きい場合：統計物理学

　このように、Nに非常に大きな数値を代入すると、ランダム性に由来する揺らぎは測定不能なほど小さくなります。

　この場合に観測される可能性が最も高いマクロ状態は、非常に限られた局所状態になります。

　実在の物質において物理的特性を測定した場合、測定結果はほとんどばらつきのない値が得られます。

　これは、アボガドロ数＝6.02\(\times 10^{23}\)オーダーの多数の原子や分子で構成された大きな N の効果によるものです。

　このような状態を記述する物理学を統計物理学といいます。

さらに数を増やしましょう

Plot_Omega(1000)

---------------------------------------------------------------------------
OverflowError                             Traceback (most recent call last)
Cell In[5], line 1
----> 1 Plot_Omega(1000)

Cell In[2], line 25, in Plot_Omega(N)
     22 Omega = []
     24 for n in N_series:
---> 25     Omega.append(Multiplicity(N,n))
     27 plt.plot(N_series, Omega)
     28 plt.title('N='+str(N))

Cell In[2], line 16, in Multiplicity(N, n)
     15 def Multiplicity(N,n):
---> 16     return factor(N)/factor(n)/factor(N-n)

OverflowError: integer division result too large for a float

Diagnostics（診断）

プログラムは除算について問題があるようです。

N1/N2

実は、2 つの大きな整数の間で除算を行うことはできません。そこで、2 つの数値を除算する代わりに、次のようにします。 log() 関数を使用します。

log(N1/N2) = log(N1) - log(N2)

最終的に、次の式が得られます。

exp(log(N1) - log(N2))

import matplotlib.pyplot as plt
from math import log, exp

#Multiplicity calculation
def factor(N):
    N = int(N)
    total = 1
    if N<0:
        print('Error, N must be greater than 0')
    elif N>1:
        for i in range(2,N+1):
            total *= i
    return total

def Multiplicity(N,n):
    #return factor(N)/factor(n)/factor(N-n)
    return exp(log(factor(N)) - log(factor(n)) - log(factor(N-n))) #change it to log funtion

def Plot_Omega2(N):

    N_series = range(N+1)
    Omega = []

    for n in N_series:
        Omega.append(Multiplicity(N,n))

    plt.plot(N_series, Omega)
    plt.title('N='+str(N))
    plt.xlabel('n')
    plt.ylabel('$\Omega(n)$')

    plt.show()

Plot_Omega2(10)

Plot_Omega2(100)

Plot_Omega2(1000)

課題１：固体のアインシュタインモデル

　固体結晶は原子で構成されています。 アインシュタインは、各原子を調和振動子と考えるモデルを提案しました。

　このモデルでは、結晶格子内の隣接する原子はバネによって繋がれていると考えます。

　各振動子で、エネルギーは量子化された状態 (0,1,2,…..) をとります。

　振動子が 2 つしかないと仮定すると、エネルギーは、

Oscillator	1	2
Energy	0	0
	1	0
	0	1
	1	1
	2	0
	0	2

つまり、全エネルギー \(q\) が0のとき, \(\Omega\)=1、全エネルギー \(q\) が1のとき, \(\Omega\)=2、全エネルギー \(q\) が2のとき, \(\Omega\)=3

振動子の数をNとしたときの一般式は、

\[\Omega(N, q) = \binom{q+N-1}{q} = \frac{(q+N-1)!}{q!(N-1)!}\]

この式を用いて、任意の (\(N\), \(q\)) について \(\Omega\) を計算するプログラムを作成してください。

Plot_Omega3(100,100)

課題２：相互作用する 2 つのアインシュタイン固体

　ここで、相互作用する 2 つのアインシュタイン固体について考えてみましょう。\(A\) と \(B\) が弱く結合している (ちょうど理想気体のように) と仮定すると、それぞれの固体のエネルギー \(q_A\) と \(q_B\) はゆっくりと変化します。

　この仮定の下では、エネルギーの総数 \(q_\text{total}\) は単純に \(q_A\) と \(q_B\) の合計になります。簡単のために、\(q_\text{total}\) を固定しましょう。任意の \(q_A\) の \(\Omega\) はどうなるでしょうか。

\(A\) の振動子の数を \(N_A\) として \(A\) だけを数えると、

\[\Omega(A) = \binom{q_A+N_A-1}{q_A} = \frac{(q_A+N_A-1)!}{q_A!(N_A-1)!}\]

同様に \(B\) の振動子の数を \(N_B\) として、

\[\Omega(B) = \binom{q_B+N_B-1}{q_B} = \frac{(q_B+N_B-1)!}{q_B!(N_B-1)!}\]

これを用いると総数は以下のようになります。

\[\Omega(\text{total}) = \Omega(A) \Omega(B)\]

\(q_A + q_B = 2, N_A = N_B = 3\)の場合、結果は以下のような表にまとめられます。

q_A	\(\Omega\)(A)	q_B	\(\Omega\)(B)	\(\Omega\)(total)
0	1	2	6	6
1	3	1	3	9
2	6	0	1	6

この \(\Omega(total)\) を \(q_A\) に対してプロットすると以下のようになります。

import matplotlib.pyplot as plt

qa = [0, 1, 2]
omega = [6, 9, 6]

plt.plot(qa, omega, '.r')
plt.xlabel('$q_A$')
plt.ylabel('$\Omega(total)$')

Text(0, 0.5, '$\\Omega(total)$')

任意の\(N_A\)、\(N_B\)、\(q_\text{total}\)について、\(\Omega(\text{total})\) を \(q_A\) の関数として計算するコードを書いてください。

次に以下の２つの場合のグラフを作成し、\(N\) が増加したときの傾向を調べてください。

1, \(N_A\)=300、\(N_B\)=200、および \(q_\text{total}\)=100

2, \(N_A\)=3000、\(N_B\)=2000、および \(q_\text{total}\)=100

Plot_Omega(300,200,100)

Plot_Omega(3000,2000,100)

Plot_Omega(300,300,100)