中性子の基本的性質¶
もういろいろとめんどくさくなってきたからまとめておく． 計算ミスがあったら教えてください． あと原則，冷中性子を対象とした計算であることに留意すること．
とりあえず基本定数をまとめると
In [1]:
from math import pow, sqrt, pi, exp
from scipy import constants as phys
import handcalcs.render
In [2]:
%%render
#Parameters
e = phys.e*1e19 # $\times 10^{19}$ [c]
c = phys.c*1e-8 # $\times 10^{-8}$ [m/s]
m = phys.m_n*phys.c**2/phys.e*1e-6 # MeV
h = phys.h/phys.e *1e15 # $\times 10^{-15}$ []
k_B = phys.k * 1e23 # $\times 10^{-23}$ []
N_A = phys.N_A*1e-23 # $\times 10^{23}$ [$\rm /m^3$]
\[
\begin{aligned}
e &= 1.602 \; \;\textrm{($\times 10^{19}$ [c])}
&c &= 2.998 \; \;\textrm{($\times 10^{-8}$ [m/s])}
&m &= 939.565 \; \;\textrm{(MeV)}
\\[10pt]
h &= 4.136 \; \;\textrm{($\times 10^{-15}$ [])}
&k_{B} &= 1.381 \; \;\textrm{($\times 10^{-23}$ [])}
&N_{A} &= 6.022 \; \;\textrm{($\times 10^{23}$ [$\rm /m^3$])}
\\[10pt]
\end{aligned}
\]
In [3]:
# unify dimensions
e = phys.e # c
h = phys.h/e # eV
hbar = phys.hbar/e # eV
c = phys.c # m/s
m = phys.m_n*c**2/e # eV/c/c
k_B = phys.k #
N_A = phys.N_A
波長，エネルギー，温度，速度の換算¶
ある運動エネルギー $E$ を持った中性子は $$ E = \frac{\hbar^2 k^2}{2m} = \frac{h^2}{2m\lambda^2} = \frac{1}{2}mv^2 = k_BT $$ で書けるから,ある波長$\lambda$の時のエネルギーと温度は
In [4]:
%%render
lamb = 0.18e-9 # m
E = h**2/(2.*(m/c**2)*lamb**2)*1e3 #meV
T = E*1e-3/(k_B/e) #K
\[
\begin{aligned}
\lambda &= 0.000 \; \;\textrm{(m)}
\\[10pt]
E &= \frac{ \left( h \right) ^{ 2 } }{ 2. \cdot \left( \frac{ m }{ \left( c \right) ^{ 2 } } \right) \cdot \left( \lambda \right) ^{ 2 } } \cdot 1 \times 10 ^ {3} \\&= \frac{ \left( 0.000 \right) ^{ 2 } }{ 2. \cdot \left( \frac{ 939565420.520 }{ \left( 299792458.000 \right) ^{ 2 } } \right) \cdot \left( 0.000 \right) ^{ 2 } } \cdot 1 \times 10 ^ {3} \\&= 25.248 \; \;\textrm{(meV)}\\[10pt]
\\[10pt]
T &= E \cdot \frac{ 1 \times 10 ^ {-3} }{ \frac{ k_{B} }{ e } } = 25.248 \cdot \frac{ 1 \times 10 ^ {-3} }{ \frac{ 0.000 }{ 0.000 } } &= 292.993 \; \;\textrm{(K)}
\end{aligned}
\]
ある波長の中性子の速度$v$は $$ \lambda = \frac{h}{mv} $$ だから
In [5]:
import sympy
sympy.init_printing()
sympy.var('lamb, v')
sympy.Eq(lamb, h/(m/c/c)*1e9/v)
Out[5]:
$\displaystyle lamb = \frac{395.603401207146}{v}$
となる．
しれっと中性子エネルギーと温度の換算をしたが，実際のところ中性子の温度分布はボルツマン分布に従った幅のある分布を持っているため，温度とエネルギーは1対1で決めることはできない． 上記計算で求めた温度は平均値に相当．
TOF法¶
パルス中性子源で用いられるtof(time of fright)と対応させるには，小学生の時習った$v = L/\rm TOF$を使えばいいから
In [6]:
sympy.var('lamb, TOF, L')
sympy.Eq(lamb, (h/(m/c/c)*1e9 /L)*TOF)
Out[6]:
$\displaystyle lamb = \frac{395.603401207146 TOF}{L}$
ちなみにこの395.6という数字，よく使うから覚えといた方がいいかもしれない．
余談だが波長を表す係数にlambdaを用いずlambと書く理由は、pythonによってlambdaという関数が定義されているから変数名として扱えない、といった背景がある。
頑張ったら使えるんだろうけど、面倒くさいのでご愛嬌ということで
数密度の計算¶
数密度は，密度$\rho$, 原子量$M$, アボガドロ数$N_A$を用いて $$ N = \frac{\rho}{M}N_A $$ で書ける．
Niの場合，
In [7]:
%%render
#Parameters
rho_Ni = 8.908 #g/cm^3
M_Ni = 58.6934 #g/mol
N_Ni = rho_Ni/M_Ni*N_A*1e6*1e-28 #1e28/m^3
\[
\begin{aligned}
\rho_{Ni} &= 8.908 \; \;\textrm{(g/cm^3)}
&M_{Ni} &= 58.693 \; \;\textrm{(g/mol)}
&N_{Ni} &= 9.140 \; \;\textrm{(1e28/m^3)}
\\[10pt]
\end{aligned}
\]
有効ポテンシャルは $$ U = \frac{2\pi\hbar^2}{m}bN $$ であるから，
In [8]:
%%render
# Parameters
b_Ni = 10.3 #fm
U_Ni = (2*pi*hbar**2)/(m/c**2) * (b_Ni*1e-15)*N_Ni * 1e9 #neV
\[
\begin{aligned}
b_{Ni} &= 10.300 \; \;\textrm{(fm)}
&U_{Ni} &= 0.000 \; \;\textrm{(neV)}
\end{aligned}
\]
化合物の場合¶
化合物の場合、有効ポテンシャルは $$ U = \frac{2\pi\hbar^2}{m}\sum_i b_iN_i $$
と書ける。
つまり各元素の数密度がわかればいい。
化合物の数密度は
$$
N_i = \frac{\rho_\mathrm{tot}}{M_i}N_A\times \mathrm{wt\%}_i
$$
化合物の場合、物質量と散乱長は組成式$A_xB_yC_z$から
$$
M_\mathrm{tot} = xM_A + yM_B + zM_C
$$ $$
b_\mathrm{tot} = xb_A + yb_B + zb_C
$$
と書ける。多分高校の化学の教科書に載ってる
この時の有効ポテンシャルは
$$
U = \frac{2\pi\hbar^2}{m}b_\mathrm{tot} \frac{\rho_\mathrm{tot}}{M_\mathrm{tot}}N_A
$$
適当に$\mathrm{SiO_2}$の場合を考える
In [9]:
#rho_Si = 2.329 #g/cm^3
M_Si = 28.0855 # g/mol
#rho_O = 1.429 #g/L
M_O = 15.9994 # g/mol
M_SiO2 = M_O*2 + M_Si
rho_SiO2 = 2.648 #g/cm^3
#rho_SiO2 = 2.2 #g/cm^3
N_SiO2 = rho_SiO2/M_SiO2*N_A*1e6
print("N_SiO2 =", N_SiO2,"kg/mol")
N_SiO2 = 2.654042525664775e+28 kg/mol
In [10]:
%%render
N_SiO2 #/m^3
b_Si = 4.1491e-15 #m
b_O = 5.803e-15 #m
b_SiO2 = (b_Si*1 + b_O*2)
U_SiO2 = (2*pi*hbar**2)/(m/c**2) * b_SiO2*N_SiO2 #eV
\[
\begin{aligned}
N_{SiO2} &= 26540425256647748047586983936.000 \; \;\textrm{(/m^3)}
\\[10pt]
b_{Si} &= 0.000 \; \;\textrm{(m)}
\\[10pt]
b_{O} &= 0.000 \; \;\textrm{(m)}
\\[10pt]
b_{SiO2} &= 0.000 \;
\\[10pt]
U_{SiO2} &= \frac{ 2 \cdot \pi \cdot \left( \mathrm{hbar} \right) ^{ 2 } }{ \frac{ m }{ \left( c \right) ^{ 2 } } } \cdot b_{SiO2} \cdot N_{SiO2} \\&= \frac{ 2 \cdot 3.142 \cdot \left( 0.000 \right) ^{ 2 } }{ \frac{ 939565420.520 }{ \left( 299792458.000 \right) ^{ 2 } } } \cdot 0.000 \cdot 26540425256647748047586983936.000 \\&= 0.000 \; \;\textrm{(eV)}\\[10pt]
\end{aligned}
\]
In [11]:
(2*pi*hbar**2)/(m/c**2)*(b_Si*rho_SiO2/M_Si*N_A*2/3. + b_O*rho_SiO2/M_O*N_A*1/3.)*1e15 #neV
Out[11]:
$\displaystyle 91.0976114896725$
In [12]:
%%render
E = 0.025 #eV
delta_n = 1 - sqrt(1 - U_SiO2/E)
\[
\begin{aligned}
E &= 0.025 \; \;\textrm{(eV)}
\\[10pt]
\delta_{n} &= 1 - \sqrt { 1 - \frac{ U_{SiO2} }{ E } } = 1 - \sqrt { 1 - \frac{ 0.000 }{ 0.025 } } &= 0.000
\end{aligned}
\]
波長を用いて書き直すと
In [13]:
%%render
U = U_SiO2 #eV
h #eV*s
m # eV
lam = 1e-9 #m
n = sqrt(1-U/(h**2/2/(m/c/c)/lam**2))
\[
\begin{aligned}
U &= 0.000 \; \;\textrm{(eV)}
\\[10pt]
h &= 0.000 \; \;\textrm{(eV*s)}
\\[10pt]
m &= 939565420.520 \; \;\textrm{(eV)}
\\[10pt]
\mathrm{lam} &= 0.000 \; \;\textrm{(m)}
\\[10pt]
n &= \sqrt { 1 - \frac{ U }{ \frac{ \left( h \right) ^{ 2 } }{ 2 } \cdot \frac{1} { \frac{ m }{ c } \cdot \frac{1} { c } } \cdot \frac{1} { \left( \mathrm{lam} \right) ^{ 2 } } } } \\&= \sqrt { 1 - \frac{ 0.000 }{ \frac{ \left( 0.000 \right) ^{ 2 } }{ 2 } \cdot \frac{1} { \frac{ 939565420.520 }{ 299792458.000 } \cdot \frac{1} { 299792458.000 } } \cdot \frac{1} { \left( 0.000 \right) ^{ 2 } } } } \\&= 1.000 \\[10pt]
\end{aligned}
\]
In [14]:
print(n)
0.9999334476148128
屈折率はほとんど1に近い値を取るため、しばしば近似した式が用いられる。 $$ \sqrt{1+x} \simeq 1+\frac{x}{2} $$ であるから、 $$ n \simeq 1 - \frac{V}{2E} = 1 - \frac{b_c N \lambda^2}{2\pi} = 1 - \frac{m U_\mathrm{Fermi} \lambda^2}{h^2} $$ と書ける。
反射臨界角¶
中性子は有効ポテンシャルによって光学的な反射を行う．
その反射角は $$ \theta_c = \sqrt{\frac{bN}{\pi}}\lambda $$ Niの場合，$\theta_c$は
In [15]:
%%render
b_Ni = 10.3e-15 #m
N_Ni #n/m^3
lam = 1e-9 #m
theta_c = sqrt(b_Ni*N_Ni/pi)*lam # rad/nm
\[
\begin{aligned}
b_{Ni} &= 0.000 \; \;\textrm{(m)}
\\[10pt]
N_{Ni} &= 9.140 \; \;\textrm{(n/m^3)}
\\[10pt]
\mathrm{lam} &= 0.000 \; \;\textrm{(m)}
\\[10pt]
\theta_{c} &= \sqrt { b_{Ni} \cdot \frac{ N_{Ni} }{ \pi } } \cdot \mathrm{lam} = \sqrt { 0.000 \cdot \frac{ 9.140 }{ 3.142 } } \cdot 0.000 &= 0.000 \; \;\textrm{(rad/nm)}
\end{aligned}
\]
また、この時0.5eV中性子の反射臨界角は
In [16]:
%%render
E = 0.5 #eV
lam = h*sqrt(1/2/m*c**2/E)*1e9 #nm
theta = theta_c*lam*1e3 #mrad
\[
\begin{aligned}
E &= 0.500 \; \;\textrm{(eV)}
\\[10pt]
\mathrm{lam} &= h \cdot \sqrt { \frac{ 1 }{ 2 } \cdot \frac{1} { m } \cdot \frac{ \left( c \right) ^{ 2 } }{ E } } \cdot 1 \times 10 ^ {9} \\&= 0.000 \cdot \sqrt { \frac{ 1 }{ 2 } \cdot \frac{1} { 939565420.520 } \cdot \frac{ \left( 299792458.000 \right) ^{ 2 } }{ 0.500 } } \cdot 1 \times 10 ^ {9} \\&= 0.040 \; \;\textrm{(nm)}\\[10pt]
\\[10pt]
\theta &= \theta_{c} \cdot \mathrm{lam} \cdot 1 \times 10 ^ {3} = 0.000 \cdot 0.040 \cdot 1 \times 10 ^ {3} &= 0.000 \; \;\textrm{(mrad)}
\end{aligned}
\]
低エネルギー中性子の断面積¶
中性子の全散乱断面積$\sigma$は中性子散乱長$b$から得られる．
$$ \sigma = 4\pi b^2 $$
中性子散乱長は干渉性散乱長(coherent scat. length)$b_c$と、非干渉性散乱長(incoherent scat. length)$b_i$の二つに分離でき、
$$ b^2 = (b_c^2+b_i^2) $$
と書ける。 ここで干渉性散乱は各原子の散乱の総和を差し、非干渉性散乱は単一の原子核からの影響をあらわす。
ほかにも非弾性散乱とかもまとめたいが追々。。。
また、吸収断面積は中性子の速度$v$を用いて $$ \sigma \propto 1/v $$ に従う．
In [17]:
%%render
rate_10B = 0.199;
rate_11B = 0.801;
M_10B = 10.0129370; # g/mol
M_11B = 11.0093054;# g/mol
M_natB = M_10B*rate_10B + M_11B*rate_11B # g/mol
rho_natB = 2.08 #g/cm^3 by wiki
N_natB = rho_natB/M_natB*N_A*1e6 #/m^3
sigma_10B = 3835 #barn
sigma_11B = 0.0055 #barn
sigma_natB = sigma_10B*rate_10B + sigma_11B*rate_11B #barn
t_natB = 0.12 #mm
factor = -N_natB*t_natB*1e-3*sigma_natB*1e-28
T = exp(factor) #ratio of transmit
\[
\begin{aligned}
\mathrm{rate}_{10B} &= 0.199 \;
\\[10pt]
\mathrm{rate}_{11B} &= 0.801 \;
\\[10pt]
M_{10B} &= 10.013 \; \;\textrm{(g/mol)}
\\[10pt]
M_{11B} &= 11.009 \; \;\textrm{(g/mol)}
\\[10pt]
M_{natB} &= M_{10B} \cdot \mathrm{rate}_{10B} + M_{11B} \cdot \mathrm{rate}_{11B} = 10.013 \cdot 0.199 + 11.009 \cdot 0.801 &= 10.811 \; \;\textrm{(g/mol)}
\\[10pt]
\rho_{natB} &= 2.080 \; \;\textrm{(g/cm^3 by wiki)}
\\[10pt]
N_{natB} &= \frac{ \rho_{natB} }{ M_{natB} } \cdot N_{A} \cdot 1 \times 10 ^ {6} = \frac{ 2.080 }{ 10.811 } \cdot 602214075999999987023872.000 \cdot 1 \times 10 ^ {6} &= 115863659574061999693863321600.000 \; \;\textrm{(/m^3)}
\\[10pt]
\sigma_{10B} &= 3835 \; \;\textrm{(barn)}
\\[10pt]
\sigma_{11B} &= 0.005 \; \;\textrm{(barn)}
\\[10pt]
\sigma_{natB} &= \sigma_{10B} \cdot \mathrm{rate}_{10B} + \sigma_{11B} \cdot \mathrm{rate}_{11B} = 3835 \cdot 0.199 + 0.005 \cdot 0.801 &= 763.169 \; \;\textrm{(barn)}
\\[10pt]
t_{natB} &= 0.120 \; \;\textrm{(mm)}
\\[10pt]
\mathrm{factor} &= \left( - N_{natB} \right) \cdot t_{natB} \cdot 1 \times 10 ^ {-3} \cdot \sigma_{natB} \cdot 1 \times 10 ^ {-28} \\&= \left( - 115863659574061999693863321600.000 \right) \cdot 0.120 \cdot 1 \times 10 ^ {-3} \cdot 763.169 \cdot 1 \times 10 ^ {-28} \\&= -1.061 \\[10pt]
\\[10pt]
T &= \exp \left( \mathrm{factor} \right) = \exp \left( -1.061 \right) &= 0.346 \; \;\textrm{(ratio of transmit)}
\end{aligned}
\]
In [18]:
print(str((1.-T)*100) + "% が止まる")
65.39192682396823% が止まる
ちょっと現実的な系¶
次に，少し現実的な計算として，遮蔽によく使われるB4C焼結体の透過率を計算する．
In [19]:
%%render
roi_t = 8 # mm
M_B4C = 55.255 # g/mol by wiki
rho_B4C = 2.52 # g/cm^3
N_B4C = rho_B4C/M_B4C*N_A*1e6 #/m^3
E = 1 #eV
roi_sigma_B = sigma_natB/sqrt(E/0.025) #barn
sigma_C = 0.0035#barn
roi_sigma_C = sigma_C/sqrt(E/0.025) #barn
roi_sigma = roi_sigma_B*4 + roi_sigma_C #barn
factor = -N_B4C*roi_t*1e-3*(roi_sigma)*1e-28
T = exp(factor)
\[
\begin{aligned}
\mathrm{roi}_{t} &= 8 \; \;\textrm{(mm)}
\\[10pt]
M_{B4C} &= 55.255 \; \;\textrm{(g/mol by wiki)}
\\[10pt]
\rho_{B4C} &= 2.520 \; \;\textrm{(g/cm^3)}
\\[10pt]
N_{B4C} &= \frac{ \rho_{B4C} }{ M_{B4C} } \cdot N_{A} \cdot 1 \times 10 ^ {6} = \frac{ 2.520 }{ 55.255 } \cdot 602214075999999987023872.000 \cdot 1 \times 10 ^ {6} &= 27465016225137994954405052416.000 \; \;\textrm{(/m^3)}
\\[10pt]
E &= 1 \; \;\textrm{(eV)}
\\[10pt]
roi_{\sigma_{B}} &= \frac{ \sigma_{natB} }{ \sqrt { \frac{ E }{ 0.025 } } } = \frac{ 763.169 }{ \sqrt { \frac{ 1 }{ 0.025 } } } &= 120.668 \; \;\textrm{(barn)}
\\[10pt]
\sigma_{C} &= 0.004 \; \;\textrm{(barn)}
\\[10pt]
roi_{\sigma_{C}} &= \frac{ \sigma_{C} }{ \sqrt { \frac{ E }{ 0.025 } } } = \frac{ 0.004 }{ \sqrt { \frac{ 1 }{ 0.025 } } } &= 0.001 \; \;\textrm{(barn)}
\\[10pt]
roi_{\sigma} &= roi_{\sigma_{B}} \cdot 4 + roi_{\sigma_{C}} = 120.668 \cdot 4 + 0.001 &= 482.671 \; \;\textrm{(barn)}
\\[10pt]
\mathrm{factor} &= \left( - N_{B4C} \right) \cdot \mathrm{roi}_{t} \cdot 1 \times 10 ^ {-3} \cdot roi_{\sigma} \cdot 1 \times 10 ^ {-28} \\&= \left( - 27465016225137994954405052416.000 \right) \cdot 8 \cdot 1 \times 10 ^ {-3} \cdot 482.671 \cdot 1 \times 10 ^ {-28} \\&= -10.605 \\[10pt]
\\[10pt]
T &= \exp \left( \mathrm{factor} \right) = \exp \left( -10.605 \right) &= 0.000
\end{aligned}
\]
In [20]:
print(str((1.0-T)*100) + "% が止まる")
99.99752146874712% が止まる
遮蔽の厚さに対する吸収率を図にする
In [21]:
import matplotlib.pyplot as plt
import numpy as np
x = np.logspace(-1, 1, 100, 'base=10') ## mm
y_val = -N_B4C*x*1e-3*(roi_sigma)*1e-28
y =np.exp(y_val)
plt.plot(x, y, label="trans")
plt.xlabel('thickness [mm]')
plt.ylabel('ratio of transmit')
plt.grid()
plt.xscale('log')
plt.yscale('log')
plt.show()
In [22]:
factor = -8e28*10*1e-3*8*1e-28
T = exp(factor) #ratio of transmit
print(T)
0.5272924240430487
運動量移行の計算¶
運動量移行$q$はその波長$\lambda$，移動角$\theta$を用いて $$ q = 4\pi\frac{\sin\theta}{\lambda} $$ と書ける。 ここで$\theta$とは、透過ビームと反射ビームの角度差の1/2であることに留意。
In [23]:
#final update
!echo final update
!date '+%b %d %Y (%a)'
final update 12 27 2023 (水)
In [ ]:
In [ ]:
In [ ]: