nano_exit

基礎的なことこそ、簡単な例が必要だと思うのです。

光学定理(複素ポテンシャル):実波数

前回、実ポテンシャルに対する光学定理を導いた。
koideforest.hatenadiary.com

今回は、複素ポテンシャルだが、波数 kは実数である時の光学定理を考える。
これは、無限遠方ではポテンシャルの虚部が完全に無くなっている場合に対応する。

導出は前回とほぼほぼ同じで、違うのは連続の式にポテンシャルの虚部が入ってくることである。
koideforest.hatenadiary.com

\displaystyle
\frac{ \partial \rho }{ \partial t } + \nabla \cdot {\bf j} = \rho \frac{2}{\hbar} \Im V
定常状態では、密度の時間微分は落ちるので、

\displaystyle
\nabla \cdot {\bf j} = \rho \frac{2}{\hbar} \Im V
\\
\displaystyle
\therefore
\int_{\Omega_r} d\vec{r} \, \nabla \cdot {\bf j} = \int_{\Omega_r} d\vec{r} \, \rho \frac{2}{\hbar} \Im V
\\
\displaystyle
\lim_{ r \rightarrow \infty} r^2 \int d\hat{r} \, {\bf j} \cdot \hat{r} = \frac{2}{\hbar} \int d\vec{r} \, \rho \Im V

確率流密度の積分は、前回の結果をそのまま流用出来るので、

\displaystyle
\pm \sigma^{\pm}_{ela} \mp \frac{4\pi}{k} \Im\left[ f^{\pm}( \pi/2 \mp \pi/2) \right] = \frac{2}{\hbar} \int d\vec{r} \, \rho \Im V

したがって、

\displaystyle
\pm \sigma^{\pm}_{ela} - \frac{2}{\hbar} \int d\vec{r} \, \rho \Im V = \pm \frac{4\pi}{k} \Im\left[ f^{\pm}( \pi/2 \mp \pi/2) \right] 
\\
\displaystyle
\sigma^{\pm}_{ela} + \sigma^{\pm}_{ine} = \frac{4\pi}{k} \Im\left[ f^{\pm}( \pi/2 \mp \pi/2) \right] = \sigma^{\pm}_{tot}
\\
\displaystyle
\sigma^{\pm}_{ine} \equiv \mp \frac{2}{\hbar} \int d\vec{r} \, \rho \Im V

 k = \sqrt{ E - V } [\rm Ryd ]であることを考えば、 \Im V <0のときに e^{ikr} = e^{ik'r}e^{-k''r}\, (k = k' + i k'')となって減衰し、 \sigma^+_{ine} > 0と対応していることがわかる。
このとき、 \sigma^-_{ine} < 0となって確率が増幅するが、これは外向波の時間反転が内向波に対応することを考えると納得出来ると思う。

  • 入射波 -> 非弾性散乱で減衰 -> 外向散乱波
  • 内向散乱波 -> 非弾性散乱で増幅 -> 入射波

光学定理(実ポテンシャル):外向波と内向波

実ポテンシャル、つまり非弾性散乱が無い時の光学定理を導く。

散乱波を次のように定義する。

\displaystyle
\psi^{\pm}( \vec{ r } )
  = A \left( e^{ i \vec{k} \cdot \vec{r} } + \frac{ f^{\pm}( \theta, \phi ) }{ r } e^{ \pm i k r} \right)
  \equiv \psi_{in}( \vec{ r } ) + \psi_{out}^{\pm}( \vec{ r } )
 \pmはそれぞれ外向波、内向波を表す。

球座標に対する \nabla

\displaystyle
\nabla = \hat{r} \frac{ \partial }{\partial r} + \frac{ 1 }{ r } u( \theta, \phi )
\\
\displaystyle
 u( \theta, \phi ) = \vec{e}_{\theta} \frac{ \partial }{ \partial \theta } + \frac{ \vec{e}_{\phi} }{ \sin\theta } \frac{ \partial }{ \partial \phi }

入射波に対する勾配は、

\displaystyle
\nabla \psi_{in} = i k \hat{k} \psi_{in}

散乱波に対する勾配は、 r \rightarrow \inftyにおいて O( r^{-1} ) だけ残すと、

\displaystyle
\frac{1}{A} \nabla \psi^{\pm}_{out}
  = \frac{ f^{\pm} }{ r } \hat{r} \frac{ \partial e^{\pm i kr } }{\partial r}
    + f^{\pm} e^{\pm i kr } \hat{r} \frac{ \partial \frac{1}{r} }{\partial r} 
    + \frac{e^{\pm i kr }}{r} \frac{1}{r} ( u f^{\pm} )
\\
\displaystyle
\quad
  = \pm i k \hat{r} \frac{ f^{\pm} }{ r } e^{\pm i kr } + O\left( \frac{1}{r^2} \right)
  \rightarrow \pm i k \hat{r} \psi^{\pm}_{out}

したがって、確率流密度は

\displaystyle
{\bf J}^{\pm}
  = \frac{\hbar}{m} \Im
    \left( ( \psi_{in} + \psi^{\pm}_{out} )^* i k ( \hat{k} \psi_{in} \pm \hat{r} \psi^{\pm}_{out} )\right)
\\
\displaystyle
\quad
  = \frac{\hbar k }{m} |\psi_{in}|^2 \hat{k} \pm \frac{\hbar k }{m} |\psi^{\pm}_{out}|^2 \hat{r}
    + \frac{\hbar }{m} \Im \left( \pm i k \psi_{in}^* \psi^{\pm}_{out} \hat{r} + i k (\psi^{\pm}_{out})^* \psi_{in} \hat{k} \right)
\\
\displaystyle
\quad
  = \frac{\hbar k }{m} |\psi_{in}|^2 \hat{k} \pm \frac{\hbar k }{m} |\psi^{\pm}_{out}|^2 \hat{r}
    + \frac{\hbar }{m} \left( \pm \Re( k \psi_{in}^* \psi^{\pm}_{out} ) \hat{r} + \Re ( k (\psi^{\pm}_{out})^* \psi_{in} ) \hat{k} \right)
\\
\displaystyle
\quad
  = \frac{\hbar k }{m} |\psi_{in}|^2 \hat{k} \pm \frac{\hbar k }{m} |\psi^{\pm}_{out}|^2 \hat{r}
    + \frac{\hbar k }{m} \left( \pm \Re( \psi_{in}^* \psi^{\pm}_{out} ) \hat{r} + \Re ( (\psi^{\pm}_{out})^* \psi_{in} ) \hat{k} \right)
\\
\displaystyle
\quad
  = \frac{\hbar k }{m} |\psi_{in}|^2 \hat{k} \pm \frac{\hbar k }{m} |\psi^{\pm}_{out}|^2 \hat{r}
    + \frac{\hbar k }{m} \Re( \psi_{in}^* \psi^{\pm}_{out} ) \left(\hat{k} \pm \hat{r} \right)
\\
\displaystyle
\quad
  \equiv {\bf J}_{in} + {\bf J}^{\pm}_{out} + {\bf J}^{\pm}_{inter}

干渉項を計算するために、平面波の球面波による漸近形を利用する。
koideforest.hatenadiary.com

\displaystyle
\frac{1}{A} \psi_{in}^*
  = e^{ - i \vec{k} \cdot \vec{r} }
  \rightarrow \frac{ 4\pi }{ 2 i k r } \left( e^{ i k r } \delta( \hat{r} - (- \hat{k} ) ) - e^{ - i k r } \delta( \hat{r} + ( - \hat{k} ) ) \right)

これにより、

\displaystyle
\psi_{in}^* \psi^+_{out} ( \hat{ r } + \hat{k} )
  \rightarrow |A|^2 \frac{ 4\pi }{ 2 i k r } \left( f^+( \pi ) e^{ 2 i k r } \delta( \hat{r} + \hat{k} ) - f^+( 0 ) \delta( \hat{r} - \hat{k} ) \right) ( \hat{ r } + \hat{k} )
\\
\displaystyle
\quad
  = |A|^2 \frac{ 4\pi }{ 2 i k r } \left( - f^+( 0 ) \delta( \hat{r} - \hat{k} ) \cdot 2 \hat{k} \right)
\\
\displaystyle
\psi_{in}^* \psi^-_{out} ( \hat{ r } - \hat{k} )
  \rightarrow |A|^2 \frac{ 4\pi }{ 2 i k r^2 } \left( f^-( \pi ) \delta( \hat{r} + \hat{k} ) - f^-( 0 ) e^{ - 2 i k r }  \delta( \hat{r} - \hat{k} ) \right) ( \hat{ r } - \hat{k} )
\\
\displaystyle
\quad
  = |A|^2 \frac{ 4\pi }{ 2 i k r^2 } \left( f^-( \pi ) \delta( \hat{r} + \hat{k} ) \cdot ( - 2 \hat{k} ) \right)
\\
\displaystyle
\therefore
{\bf J}_{inter}^{\pm}
  = - |A|^2 \frac{\hbar k}{m} \frac{ 4\pi }{ 2 k r^2 } 2 \Im \left[ f^{\pm}( \pi/2 \mp \pi/2 )\right] \delta( \hat{r} \mp \hat{k} ) \hat{k}
\\
\displaystyle
\quad
  = - \frac{ 4\pi }{ k } |A|^2 \frac{ \hbar k  }{ m } \frac{ 1 }{ r^2 } \Im \left[ f^{\pm}( \pi/2 \mp \pi/2 )\right] \delta( \hat{r} \mp \hat{k} ) \hat{k}

確率流密度が求まったので、連続の式で確率の保存を考える。
定常状態かつ実ポテンシャルの下で、

\displaystyle
\frac{ \partial \rho }{ \partial t } = 0,
\quad
\Im V = 0,
\\
\displaystyle
\therefore
\frac{ \partial \rho }{ \partial t } + \nabla \cdot {\bf J} = 0
\rightarrow \nabla \cdot {\bf J} = 0

したがって、半径 rの体積積分を取り、ガウスの定理を使えば、

\displaystyle
\int_{\Omega_r} d\vec{r} \, \nabla \cdot {\bf J}
  = r^2 \int d\hat{r} \,{\bf J} \cdot \hat{r}

したがって、各確率流密度が与える寄与は、

\displaystyle
r^2 \int d\hat{r} \,{\bf J}_{in} \cdot \hat{r}
  = r^2 \frac{ \hbar k }{ m } |A|^2 \int d\hat{r} \, \hat{k} \cdot \hat{r}
  = 0
\\
\displaystyle
r^2 \int d\hat{r} \,{\bf J}^{\pm}_{out} \cdot \hat{r}
  = \pm \frac{ \hbar k }{ m } |A|^2 \int d\hat{r} \, | f^{\pm} |^2
  = \pm \frac{ \hbar k }{ m } |A|^2 \sigma^{\pm}_{ela}
  = \pm \frac{ \hbar k }{ m } |A|^2 \sigma^{\pm}_{tot}
\\
\displaystyle
r^2 \int d\hat{r} \, {\bf J}^{\pm}_{inter} \cdot \hat{r}
\\
\displaystyle
\quad
  = - r^2 \frac{ 4\pi }{ k } |A|^2 \frac{ \hbar k }{ m } \frac{ 1 }{ r^2 } \Im \left[ f^{\pm}( \pi/2 \mp \pi/2 )\right] \int d\hat{r} \, \delta( \hat{r} \mp \hat{k} ) \hat{k} \cdot \hat{r}
\\
\displaystyle
\quad
  = \mp \frac{ 4\pi }{ k } |A|^2 \frac{ \hbar k }{ m } \Im \left[ f^{\pm}( \pi/2 \mp \pi/2 )\right]
実ポテンシャルの場合には、非弾性散乱が無いので、 \sigma_{ela} = \sigma_{tot}である。

したがって、上記の確率保存の観点から、光学定理と呼ばれる次の関係式が得られる。

\displaystyle
\pm \sigma^{\pm}_{tot} \mp \frac{ 4\pi }{ k } \Im \left[ f^{\pm}( \pi/2 \mp \pi/2 )\right] = 0
\\
\displaystyle
\therefore
\sigma^{\pm}_{tot} = \frac{ 4\pi }{ k } \Im \left[ f^{\pm}( \pi/2 \mp \pi/2 )\right]

多くの場合、外向波の条件でしか光学定理が導かれていないが、内向波だと後方散乱の虚部が全断面積を与えることがわかる。
これは、入射波と同じ向きの波が、外向波が前方だが、内向波では後方に対応することに因る。

連続の式

連続の式を導出する。

発想として、密度の時間依存性と時間依存のシュレーディンガーを方程式を結び付けることを考える。

\displaystyle
\frac{\partial \rho}{ \partial t }
  = \frac{\partial \psi^*}{ \partial t } \psi + \psi^* \frac{\partial \psi}{ \partial t }
  = \frac{1}{ i \hbar }
      \left(
          - \left( i \hbar \frac{\partial \psi } { \partial t } \right)^* \psi
          + \psi^* \left( i \hbar \frac{\partial \psi}{ \partial t }\right)
      \right)
\\
\displaystyle
i \hbar \frac{\partial \psi } { \partial t }
  = \left( \frac{ \hbar^2 }{ 2m } \nabla^2 + V \right) \psi
\\
\displaystyle
\therefore
\frac{\partial \rho}{ \partial t }
  = \frac{1}{ i \hbar }
      \left(
          - \psi \left( \frac{ \hbar^2 }{ 2m } \nabla^2 + V^* \right) \psi^*
          + \psi^* \left( \frac{ \hbar^2 }{ 2m } \nabla^2 + V \right) \psi
      \right)
\\
\displaystyle
\quad
  = - \frac{1}{ i \hbar } \frac{ \hbar^2 }{ 2m }( \psi^* \nabla^2 \psi - \psi \nabla^2 \psi^* )
    + \frac{1}{ i \hbar } ( V - V^* ) \rho
\\
\displaystyle
\quad
  = - \frac{1}{ i \hbar } \frac{ \hbar^2 }{ 2m } 2 i \Im ( \psi^* \nabla^2 \psi )
    + \frac{\rho}{ i \hbar } 2 i \Im V

ここで、ベクトル解析の関係式から、

\displaystyle
\psi^* \nabla^2 \psi = \nabla \cdot ( \psi^* \nabla \psi ) - ( \nabla \psi )^* \cdot ( \nabla \psi )  
\\
\displaystyle
\therefore
\Im ( \psi^* \nabla^2 \psi ) = \nabla \cdot \Im ( \psi^* \nabla \psi )

したがって、

\displaystyle
\frac{\partial \rho}{ \partial t }
  = - \nabla \cdot \left( \frac{1}{ i \hbar } \frac{ \hbar^2 }{ 2m } 2 i \Im ( \psi^* \nabla \psi ) \right)
    + \rho \frac{2}{ \hbar } \Im V
  \equiv - \nabla \cdot {\bf j} + \rho \frac{2}{ \hbar } \Im V
\\
\displaystyle
{\bf j}
  = \frac{1}{ i \hbar } \frac{ \hbar^2 }{ 2m } 2 i \Im ( \psi^* \nabla \psi )
  = \frac{ \hbar }{ m } \Im ( \psi^* \nabla \psi )

よって、連続の式は、

\displaystyle
\frac{\partial \rho}{ \partial t } + \nabla \cdot {\bf j}
  = \rho \frac{2}{ \hbar } \Im V

ポテンシャルが虚部を含む時( \Im V \neq 0)、確率が保存しないことがわかる。

二階微分演算子の変数変換 for Numerov method

Numerov法は一階微分方程式を含まない二階微分方程式を解くのに便利な方法である。
この方法を適用するには、等間隔の変数刻み(メッシュ or グリッド)が必要である。
しかし、変化が穏やかなところでは粗いメッシュで十分だが、一方で変化が急なところでは細かなメッシュを用いなければ誤差が大きくなってしまうため、等間隔メッシュの使用は無駄が多い。
そこで、「細かいところは細かく、そうでないところは粗い、理想的なメッシュ」と対応する「等間隔なメッシュ」に変数変換し、変換後の方程式をNumerov法で解くことを考える。

変換前と後のメッシュを r, \rhoとして、 \rho = m[ r ] という関係で結ばれているとする。
 \rhoは等間隔なので、最初の点を \rho_0として、
 
\displaystyle
\rho_n = \rho_0 + n h \quad ( n = 0, 1, \cdots, N-1 )
と与えられる。 h, Nはそれぞれ \rhoのメッシュと点数である。
ここから、

\displaystyle
r = m^{-1}[ \rho ] \quad ( m ( m^{-1}( \rho ) ) = \rho )
と逆変換したいところだが、関数系が複雑だと、解析的に書けない。
そこで、 y(r) = m[ r ] - \rho_n ととして、 y(r) = 0を与える rを数値的に求めて、それを r_nとして保存しておけば良い。
 y(r) = 0を解くための数値計算法として、Newton法が簡便であるが、ここでは立ち入らない。
ニュートン法とは何か??ニュートン法で解く方程式の近似解 - Qiita

次に、二階の微分演算子を変数変換する。
数値計算をするのは \rhoに対してだから、 \rhoでの二階微分から出発して、 rの二階微分とどう繋がるのかを見る。
 df/dx = f'とすると、

\displaystyle
\frac{d^2f}{d\rho^2}
  = \frac{d}{d\rho} \left( \frac{dr}{d\rho} f' \right)
  = \frac{d^2 r }{d \rho^2} f' + \frac{dr }{d \rho} \frac{df'}{d\rho}
  = \frac{d^2 r }{d \rho^2} f' + \left( \frac{dr }{d \rho} \right)^2 f''

重要なポイントとして、一次微分 f'が残っているとNumerov法が使えない。
そのため、 f(r) = a(r) g(r)として、上手いこと f'の項を消すような a(r)を考える
 fを欲しい関数とすると上手く行かないから、本当に欲しい関数を gとして、そこに aを掛けた仮の関数を fと置き直すという流れ)。
準備として、

\displaystyle
\frac{ dr }{ d \rho} = (\rho')^{-1}
\\
\displaystyle
\frac{ d^2r }{ d \rho^2}
  = \frac{ d }{ d \rho} (\rho')^{-1}
  = ( {\rho'} )^{-1} \left( (\rho')^{-1} \right)'
  = ( {\rho'} )^{-1} \left( - \frac{ \rho'' }{ (\rho')^2 } \right)
  = - \frac{ \rho'' }{ (\rho')^3 }
\\
f' = a' g + a g'
\\
f'' = a'' g + 2 a' g' + a g''

したがって、

\displaystyle
\frac{d^2 f}{ d \rho^2} 
  = \frac{d^2 r}{ d \rho^2 } ( a' g + a g' ) + \left( \frac{d r}{ d \rho } \right)^2 ( a'' g + 2 a' g' + a g'' )
\\
\displaystyle
\quad
  = \left( \frac{d r}{ d \rho } \right)^2 a g'' + \left( \frac{d^2 r}{ d \rho^2} a + 2 \left( \frac{d r}{ d \rho } \right)^2 a' \right) g' +  \left( \frac{d^2 r}{ d \rho^2} a' + \left( \frac{d r}{ d \rho } \right)^2 a'' \right) g

よって、 f'を残さないためには、

\displaystyle
\frac{d^2 r}{ d \rho^2} a + 2 \left( \frac{d r}{ d \rho } \right)^2 a' = 0
\\
\displaystyle
  - \frac{ \rho'' }{ (\rho')^3 } a + \frac{2}{ (\rho')^2 } a' = 0
\\
\displaystyle
  \frac{ a' }{ a } = \frac{1}{2} \frac{ \rho'' }{ \rho' }
両辺を積分すれば、

\displaystyle
\ln( a ) = \frac{1}{2} \ln( \rho' )
\\
\therefore
\displaystyle
a = \sqrt{ \rho' }
 \rho=m[r]という形で \rhoが分かっているから、その微分は容易である。

今、欲しい関数は gに置き直していて、 gの二階微分方程式 g''(r) = b(r) g(r) + G(r)で与えられているとすれば、

\displaystyle
\frac{ d^2 f }{d \rho^2}
  = \frac{ 1 }{ (\rho')^2 } a g'' + \left( - \frac{ \rho'' }{ (\rho')^3 } \frac{a'}{a} + \frac{1}{(\rho')^2} \frac{a''}{a} \right) a g
\\
\displaystyle
  = \frac{ 1 }{ (\rho')^2 } a ( b g + G ) + \left( - \frac{ \rho'' }{ (\rho')^3 } \frac{a'}{a} + \frac{1}{(\rho')^2} \frac{a''}{a} \right) f
\\
\displaystyle
  = ( c^{(1)}_g + c ) f + c^{(2)}_g G
\\
c
  = - \frac{ \rho'' }{ (\rho')^3 } \frac{a'}{a} + \frac{1}{(\rho')^2} \frac{a''}{a}
\\
c^{(1)}_g =  \frac{ b }{ (\rho')^2 }
\\
c^{(2)}_g =  \frac{ a }{ (\rho')^2 }
 c^{(1)}_g, c^{(2)}_gは元の微分方程式にあった係数が修正されたもので、一方 cは変数変換によって新たに付け足された項となる。
これをNumrov法で解けば良い。

具体例として、linear-log mesh  \rho = \alpha r + \beta \ln r を使ってみる。
係数 cを手計算するのは面倒くさいので、pythonのsympyにやらせることにする。

from sympy import *

r, alpha, beta = symbols( "r α β" )

rho = alpha * r + beta * log( r )
rhop = Derivative( rho, r ).doit()
rhopp = Derivative( rhop, r ).doit()

a = sqrt( rhop )
ap = Derivative( a, r ).doit()
app = Derivative( ap, r ).doit()

c = simplify( - rhopp / (rhop**3) * ap / a + 1/(rhop**2) * app / a )

結果

\displaystyle
c = \frac{ \beta ( \alpha r + \beta / 4 )}{ ( \alpha r + \beta )^4 }

sympyをもっと有効活用して行きたい。

存在をチェックするpythonコマンドまとめ

  • Path (File and/or Directory)
import os

path = "/test"
os.path.exists( path )  # True or False

path_file = "/test/test.dat"
os.path.isfile( path_file )  # True or False

path_dir = "/test/"
os.path.isdir( path_dir )  # True or False

地味に英語の文章になっていて、"path"が三人称単数であることを受けて、"is"だったり"exist"に三人称単数の"s"が付いたりしている。

  • List
l = [ "test" ]
target = "test"

target in l  # True

try:
    l.index( target )  # integer
except ValueError
    print( "{} does not exist".format( target ) )

l.count( target )  # integer
  • 文字
string = "test"
target = "t"

target in string  # True

断面積の単位:barn

たまに  Mb という単位を見かけたことがあったが、何だろうなと思ってスルーしていたので調べてみた。

バーン (単位) - Wikipedia

要は面積の単位でした。


\displaystyle
 1b = 10^{-28} m^2 = 10^{-24} cm^2

原子核物理の分野で、微分散乱断面積を表すのに使うそうです。
 10 fm \times 10 fmの面積だから、原子核には丁度良いらしい。
僕は普通に電子散乱なので、見てるスケールがだいぶ違うなぁと。

原子核物理で思い出したが、下の散乱理論の本は読み易くて良かった。
www.kyoritsu-pub.co.jp
物性だと、微分散乱断面積そのものが大事にされる場面は多くないので、あまり実態が感じられない。
その点、ラザフォードの流れを汲む原子核から散乱理論を扱うと、イメージと直接繋がっていて、理解し易く感じた。

他分野の手法を使うのは、そう言った意味でも、ある種のハードルがあることを改めて思った。