2018-05-05

球対称関数を別の位置で球平均する

$r$ を原点から測った距離とし、原点から見て球対称な関数を $F( r )$ とする。
原点から見て位置 ${\bf R}_i$ にあるサイト $i$ があるとする。
サイト $i$ を原点に取り直した任意の位置ベクトルを ${\bf r}_i$ と定義する。
やりたいことは、 $F( r )$ をサイト $i$ から見ると球対称ではないので、サイト $i$ から見たときの球対称成分 $f( r_i )$ を作りたいということである。
それはすなわちサイト $i$ 上で球平均することに対応する。
${\bf r}_i = {\bf r} - {\bf R}_i$ 、もしくは ${\bf r} = {\bf r}_i + {\bf R}_i$ （こっちの方が直観的かも）であるから、

$\displaystyle \begin{align} f( r_i ) = \left( \int d \hat{\bf r}_i \, F( | {\bf r}_i + {\bf R}_i | ) \right) / \left( \int d \hat{\bf r}_i \right) = \frac{ 1 }{ 4 \pi } \int d \hat{\bf r}_i \, F( | {\bf r}_i + {\bf R}_i | ) \end{align}$

球対称であるから、 ${\bf R}_i$ を $z$ 軸上に取っても結果は変わらない。
この時、 $F( | {\bf r}_i + {\bf R}_i | )$ は方位角に依存しないから、

$\displaystyle \begin{align} f( r_i ) = \frac{ 1 }{ 2 } \int^{\pi}_{0} {\rm sin} \theta \, d \theta \, F( | {\bf r}_i + {\bf R}_i | ) \end{align}$

ここで、 $\theta$ は ${\bf r}_i$ と $- {\bf R}_i$ の成す角である。
マイナスが気になるかも知れないが、最終的にはその動径成分の $R_i$ のみを使うので、符号はあまり重要でない。
$| {\bf r}_i + {\bf R}_i | = \sqrt{ r^2_i + R^2_i - 2 r_i R_i {\rm cos} \theta }$ であることに注意すると、

$\displaystyle \begin{align} &\frac{ d | {\bf r}_i + {\bf R}_i | }{ d \theta } = \frac{ r_i R_i }{ \sqrt{ r^2_i + R^2_i - 2 r_i R_i {\rm cos} \theta } } {\rm sin} \theta = \frac{ r_i R_i }{ | {\bf r}_i + {\bf R}_i | } {\rm sin} \theta \\ &\therefore \; {\rm sin} \theta \, d \theta = \frac{ | {\bf r}_i + {\bf R}_i | }{ r_i R_i } \, d | {\bf r}_i + {\bf R}_i | \end{align}$

よって、

$\displaystyle \begin{align} f( r_i ) &= \frac{ 1 }{ 2 r_i R_i } \int^{ R_i + r_i }_{ R_i - r_i } \, d | {\bf r}_i + {\bf R}_i | \, | {\bf r}_i + {\bf R}_i | \, F( | {\bf r}_i + {\bf R}_i | ) \\ &= \frac{ 1 }{ 2 r_i R_i } \int^{ R_i + r_i }_{ R_i - r_i } \, dr \, r \, F( r ) \end{align}$

$r_i$ の関数として $f( r_i )$ を求めたいので、ここでの扱いは言わばパラメータのように外から与えられるものであり、積分変数ではないため積分の外に出せる。

例として、クーロンポテンシャル $F( r ) = Z / r \equiv V( r )$ を考えると、

$\displaystyle \begin{align} V( r_i ) &= \frac{ Z }{ 2 r_i R_i } \left( R_i + r_i - ( R_i - r_i ) \right) = \frac{ Z }{ R_i } \equiv \bar{ V } \end{align}$

となり、違うサイトから見ると定数ポテンシャルとして与えられることになる。
したがって、点電荷が一次元に正負（絶対値は同じ大きさ $Z$ ）で交互に等間隔 $R$ で無限に並んでいて、ある点電荷に働くトータルの静電ポテンシャルは、球平均すれば

$\displaystyle \begin{align} V = - \frac{ 2 Z }{ R } \left( 1 - \frac{ 1 }{ 2 } + \frac{ 1 }{ 3 } - \cdots \right) = - \frac{ 2 Z }{ R } {\rm ln}2 \end{align}$

と求めることが出来る。
ただ、これは $ln( 1 + x )$ のマクローリン展開、つまり $x \approx 0$ に対して速く収束する展開に、 $x = 1$ を代入して得られたものを使用しているわけだから、本当に点電荷がほぼ無限に並んでいないと（ほぼ無限に展開しないと）値の一致は悪いと思われる。クーロンポテンシャルの収束が遅いのがよく分かる。

2018-05-03

デルタ関数のエンタングルメントを完備関係式で断つ？

数学量子力学

例えば球面調和関数

$\displaystyle \begin{align} \delta( \hat{\bf r} - \hat{\bf r}' ) &= \delta( \hat{\bf r}' - \hat{\bf r} ) \\ &= \sum_L Y_{L} ( \hat{\bf r} ) Y^*_L ( \hat{\bf r}' ) = \sum_L Y^*_{L} ( \hat{\bf r} ) Y_L ( \hat{\bf r}' ) \\ &= \sum_L Y_{L} ( \hat{\bf r}' ) Y^*_L ( \hat{\bf r} ) = \sum_L Y^*_{L} ( \hat{\bf r}' ) Y_L ( \hat{\bf r} ) \end{align}$

これはある意味 $\hat{\bf r},\, \hat{\bf r}$ のエンタングルメントを断ったと言えるか？
添字が一個なので、ベクトルに拡張すれば、

$\displaystyle \begin{align} \delta( \hat{\bf r} - \hat{\bf r}' ) = \left( Y_{L_1} ( \hat{\bf r} ), Y_{L_2} ( \hat{\bf r} ), \cdots \right) \left( \begin{array}{c} Y^*_{L_1}( \hat{\bf r}' ) \\ Y^*_{L_2}( \hat{\bf r}' ) \\ \vdots \end{array} \right) = \left( Y^*_{L_1} ( \hat{\bf r} ), Y^*_{L_2} ( \hat{\bf r} ), \cdots \right) \left( \begin{array}{c} Y_{L_1}( \hat{\bf r}' ) \\ Y_{L_2}( \hat{\bf r}' ) \\ \vdots \end{array} \right) \\ = \left( Y_{L_1} ( \hat{\bf r}' ), Y_{L_2} ( \hat{\bf r}' ), \cdots \right) \left( \begin{array}{c} Y^*_{L_1}( \hat{\bf r} ) \\ Y^*_{L_2}( \hat{\bf r} ) \\ \vdots \end{array} \right) = \left( Y^*_{L_1} ( \hat{\bf r}' ), Y^*_{L_2} ( \hat{\bf r}' ), \cdots \right) \left( \begin{array}{c} Y_{L_1}( \hat{\bf r} ) \\ Y_{L_2}( \hat{\bf r} ) \\ \vdots \end{array} \right) \end{align}$

よって、ある種の行列化によって変数分離が可能である。ポイントは、もともと行列ではなくただの数だったので、「左に横ベクトル・右に縦ベクトル」の順番は変更出来ないという点である。
あまりエンタングルメントについて詳しくない（というか無知に等しい）が、エンタングルメントを断つ（解す）とは多分こういうことなのかと。
上では非局所相関を局所相関の積で表すというようなことをしたわけだが、この結果を用いると、いわゆる多重極展開と呼ばれる、動径成分と立体角成分に分ける展開が可能になる。
（これもある種のエンタングルメント解しか？）

$\displaystyle \begin{align} F( r, \hat{\bf r} ) &= \int F( r, \hat{\bf r}' ) \delta( \hat{\bf r} - \hat{\bf r}' ) d \hat{\bf r}' = \int F( r, \hat{\bf r}' ) \sum_L Y_{L} ( \hat{\bf r} ) Y^*_L ( \hat{\bf r}' ) d \hat{\bf r}' \\ &= \sum_L \left( \int F( r, \hat{\bf r}' ) Y^*_L ( \hat{\bf r}' ) d \hat{\bf r}' \right) Y_{L}( \hat{\bf r} ) \equiv \sum_L f_L( r ) Y_{L}( \hat{\bf r} ) \\ &= \left( f_{L_1}( r ), f_{L_2}( r ), \cdots \right) \left( \begin{array}{c} Y_{L_1}( \hat{\bf r} ) \\ Y_{L_2}( \hat{\bf r} ) \\ \vdots \end{array} \right) \end{align}$

これは「連続量の相関を離散量の話に置き換える」と言えるだろうか？
ちなみに逆もおそらく可能（離散量の相関を連続量の話に置き換える）。

$\displaystyle \begin{align} \delta_{ L, L' } &= \int d \hat{\bf r} Y^*_L( \hat{\bf r} ) Y_{L'}( \hat{\bf r} ) \\ g_{N L} &= \sum_{L'} g_{NL'} \delta_{L' L} = \sum_{L'} g_{N L'} \int d \hat{\bf r} Y^*_{L'}( \hat{\bf r} ) Y_{L}( \hat{\bf r} ) \\ &= \int d \hat{\bf r} \left( \sum_{L'} g_{N L'} Y^*_{L'}( \hat{\bf r} ) \right) Y_{L}( \hat{\bf r} ) \equiv \int d \hat{\bf r} G_{N}( \hat{\bf r} ) Y_{L}( \hat{\bf r} ) \end{align}$

連続と離散の絡み合いは楽しいので、もっと深めて行きたい。

2018-05-01

Gaunt積分のまとめ

数学量子力学

Gaunt積分の自分的メモ。
Clebsch–Gordan coefficients - Wikipedia
3-j symbol - Wikipedia

Gaunt積分（多分一般的な呼び名ではない。が、自分の分野ではよく出てくる。）
$\displaystyle \begin{align} & \int d\hat{\bf r} \, Y^*_{ l_3 m_3 }(\hat{\bf r}) Y_{ l_1 m_1 }(\hat{\bf r}) Y_{ l_2 m_2 }(\hat{\bf r}) \\ & \quad = \sqrt{ \frac{ (2 l_1 + 1 ) ( 2 l_2 + 1 ) }{ 4 \pi ( 2 l_3 + 1 ) } } < l_1,0;l_2,0|l_3,0 > < l_1,m_1;l_2,m_2|l_3,m_3 > \\ & \quad = (-1)^{ m_3 } \sqrt{ \frac{ (2 l_1 + 1 ) ( 2 l_2 + 1 ) ( 2 l_3 + 1 ) }{ 4 \pi } } \left( \begin{array}{ccc} l_1 & l_2 & l_3 \\ 0 & 0 & 0 \end{array} \right) \left( \begin{array}{ccc} l_1 & l_2 & l_3 \\ m_1 & m_2 & - m_3 \end{array} \right) \end{align}$

Gaunt係数（Wikipediaとかに載っているのはこっち。）
$\displaystyle \begin{align} & \int d\hat{\bf r} \, Y_{ l_1 m_1 }(\hat{\bf r}) Y_{ l_2 m_2 }(\hat{\bf r}) Y_{ l_3 m_3 }(\hat{\bf r}) \\ & \quad = \sqrt{ \frac{ (2 l_1 + 1 ) ( 2 l_2 + 1 ) ( 2 l_3 + 1 ) }{ 4 \pi } } \left( \begin{array}{ccc} l_1 & l_2 & l_3 \\ 0 & 0 & 0 \end{array} \right) \left( \begin{array}{ccc} l_1 & l_2 & l_3 \\ m_1 & m_2 & m_3 \end{array} \right) \end{align}$

Clebsh-Gordan係数と3j記号の関係
$\displaystyle \begin{align} < l_1,m_1;l_2,m_2|l_3,m_3 > = (-1)^{ l_1 - l_2 } (-1)^{ m_3 } \sqrt{ 2 l_3 + 1 } \left( \begin{array}{ccc} l_1 & l_2 & l_3 \\ m_1 & m_2 & - m_3 \end{array} \right) \end{align}$

3j記号の列の交換
$\displaystyle \begin{align} \left( \begin{array}{ccc} l_1 & l_2 & l_3 \\ m_1 & m_2 & m_3 \end{array} \right) &= (-1)^{ l_1 + l_2 + l_3 } \left( \begin{array}{ccc} l_2 & l_1 & l_3 \\ m_2 & m_1 & m_3 \end{array} \right) \\ &= (-1)^{ l_1 + l_2 + l_3 } \left( \begin{array}{ccc} l_1 & l_3 & l_2 \\ m_1 & m_3 & m_2 \end{array} \right) \end{align}$

2018-05-01

行列要素を視覚的にチェックする方法

プログラミング数学

行列要素をなんとか可視化したくて、三次元プロット散布図で表現してみた。
Pythonでデータを可視化するmatplotlib初級（3D散布図編） - しぷぜん

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np

N = 20
row = range( N )
column = range( N )

X, Y = np.meshgrid( row, column )
X = X.reshape( N * N )
Y = Y.reshape( N * N )

fig = plt.figure()
ax = fig.add_subplot( 111, projection='3d' )
for i in range(3):
    matrix_reshape_1d = np.random.rand( N * N )
    # default marker size = 20
    ax.scatter( X, Y, [i] * ( N * N ), s = matrix_reshape_1d * 20. )
plt.show()

f:id:koideforest:20180501190738p:plain

例えば行列全体が時間に依存していて、それの時間発展を見たい時には良さそう。
複素数の場合には、色で位相を表現することになると思うが、色は配列で渡せないので、一点一点をいちいちプロットしないと行けない気がする。行列サイズが大きいと、おそらく途中で限界が来るだろう。

2018-04-29

Legendre多項式の直交性を python (scipy) で確かめる

プログラミング数学

いつもnumpyの多項式の使い方を忘れるので、Legendre多項式の直交関係式で練習する。
ルジャンドル多項式 - Wikipedia

$\displaystyle \int^{1}_{-1} P_{ l' }( x ) \, P_{ l }( x ) \, dx = \frac{ 2 }{ 2 l + 1 } \delta_{ l' l }$

scipy.special.legendreでLegendre多項式を作ると、poly1dというタイプの関数が帰ってくるので、そのままxを代入できる（下記リンク or サンプルスクリプト参照）。
scipy.special.legendre — SciPy v1.0.0 Reference Guide
numpy.poly1d — NumPy v1.14 Manual

積分は一次元積分のquadを使用。
Python SciPy : SciPy の積分関数の基本的使い方 | org-技術

from scipy.special import legendre
from scipy.integrate import quad

def Plx( l, x ):
    # legendre( l ) >>> poly1d([ ... ])
    # poly1d([ ... ])( x ) >>> float
    return legendre( l )( x )

def orthogonal_relation( lp, l ):
    # quad >>> ( result, error )
    return quad( lambda x: Plx( lp, x ) * Plx( l, x ), -1., 1. ) [0]

if __name__ == '__main__'
    for lp in range( 5 ):
        for l in range( 5 ):
            print( 'lp, l = ', lp, l )
            print ( 'numerical', orthogonal_relation( lp, l ) )
            if lp == l:
                print( 'analytical', 2. / ( 2 * l + 1. ) )
            else:
                print( 'analytical', 0. )

2018-04-26

sympyの虚数単位を通常に戻してnumpyが使えるようにする

プログラミング数学量子力学

sympy moduleを一部使ってゴニョゴニョ計算させて、それを元にnumpy moduleで処理（例えば逆行列計算）させようとした時に虚数単位が引っかかって怒られることがあった。
sympyの代数表記を数値表記に戻す方法として"N()"や"~.evalf()"があるが、これではsympyの虚数単位"I"は通常python上で扱われる"j"に戻らない。
戻すためには、"complex()"を使う必要がある。
ちなみに、Clebsh-Gordan係数のような関数丸ごとがobjectとしてsympy上で扱われている場合、"~.doit()"で一度崩してあげる必要がある。
つまり、"complex( ~.doit() )"のような感じ。

Sympyのウェブマニュアルには一箇所だけ"complex()"が出てきて、後は"N()"とか別のを使っているので、流し読みしていると見つからなかった。
Numerical evaluation — SymPy 1.2 documentation
検索機能って本当大事ですね。

追記
例えば"N()"や"evalf()"で数値化したものが実数だったとして、それに"+1j"で虚数を足すと、sympyの虚数単位"I"に変換されてしまう。
そのため、実数だと分かっている場合には、"float()"で通常の実数にし、そのあとで"+1j"を足すことが必要。

2018-03-25

平方根の行列表現

数学

平方根を評価するのに、級数展開以外のもので何か無いか考えた時に、虚数の行列表現を思い出した。

とりあえず $\sqrt{2}$ に対応する行列を求めたい。
つまり、二乗したら単位行列に2をかけたものを返す行列を考える。

$\displaystyle \begin{align} \left( \begin{array}{cc} 2 & 0 \\ 0 & 2 \end{array} \right) = \left( \begin{array}{cc} a & b \\ c & d \end{array} \right)^2 = \left( \begin{array}{cc} a^2 + bc & b ( a + d ) \\ c ( a + d ) & d^2 + bc \end{array} \right) \end{align}$

ここで非対角成分を消すために $b = c = 0$ としてしまうと、 $a^2=d^2=\sqrt{2}$ となって全く面白くない。
そのため $a + d=0$ の時の行列を探す。
この時、 $bc = 2 - a^2$ という条件しか出て来ないため、一意には決まらないことがわかる。
とりあえず、 $a = 0, c = 1$ を課せば、

$\displaystyle \begin{align} \left( \begin{array}{cc} 0 & 2 \\ 1 & 0 \end{array} \right) \end{align}$

が $\sqrt{2}$ に対応した行列になる。
$b$ と $c$ を入れ替えた転置の表式でも二乗して同じ値になる。
もし仮に $a = 1$ とした時には、

$\displaystyle \begin{align} \left( \begin{array}{cc} 1 & 1 \\ 1 & -1 \end{array} \right) = \left( \begin{array}{cc} 1 & 0 \\ 0 & -1 \end{array} \right) + \left( \begin{array}{cc} 0 & 1 \\ 1 & 0 \end{array} \right) = \sigma_3 + \sigma_1 \end{align}$