[양자역학] 7.수소꼴 원자 모형 (Hydrogenic Atoms Model)과 시각화

지난 글에서는 구면에서의 회전 운동을 기술하는 슈뢰딩거 방정식을 풀고, 그 해인 구면 조화 함수(Spherical Harmonic Functions)에 대해서 알아보았다. 이번에는 수소꼴 원자 모형에 대해서 슈뢰딩거 방정식을 풀어보자. (지난번 글과 내용이 이어지니 보고오는 것을 추천한다.)

 

구면좌표 표기법에는 physics convention과 mathmatical convention 2가지가 있는데, 아래와 같은  physics convention을 따르고 있으니 했깔리지 않도록 유의바란다.

Hydrogenic Atoms (수소꼴 원자)

수소꼴 원자는 전자가 하나밖에 없는 원자나 이온으로 \( \mathrm{ H, He^+, Li^{2+} } \)등이 해당한다. 여기서 왜 일반적인 원자 모형이 아니라 수소꼴 원자인지 궁금할 수 있는데, 헬륨 원자는 핵 1개, 전자 2개로 삼체(三體) 문제가 되어 근사없이 해석적인 일반해를 구하는 것이 불가능하기 때문이다. (참고로 앙리 푸엥카레가 삼체 문제는 해석적인 일반해를 못구한다고 증명하였다.) 

 

Effective mass (유효 질량)

$$\mu = \dfrac{m_1 m_2}{m_1 + m_2} $$

1번 입자와 2번 입자에 대한 문제를 푼다고 생각해보자. 한번에 해를 구하기 보다는 1번 입자에 대한 운동과 1번 입자가 보았을 때 2번 입자가 어떻게 움직이는 지 상대적인 운동을 기술하는 것이 쉬울 것이다. 이때 만약 두 입자사이에 인력이 작용한다면, 1번 입자가 2번 입자의 운동을 바라볼 때의 운동은 두 물체의 질량 중심에서 2번 입자를 바라본 운동에 비해 서로가 다가오기 때문에 더 큰 가속도를 가지는 것으로 보이는데, 이는 질량이 작아진 효과로 생각할 수 있다. 유효 질량을 유도해봄으로써 이를 확인해보자.

 

먼저 질량이 각각 \(m_1, m_2\)인 두 물체의 질량중심을 가르키는 위치벡터는 다음과 같다.

$$ \mathbf{r}_{\mathrm{cm}} = \dfrac{m_1 \mathbf{r}_1 + m_2 \mathbf{r}_2}{m_1 + m_2} $$

그리고 질량 중심에서 보았을때 1번 물체에 작용하는 힘\(F= ma\)에서 아래와 같이 쓸 수 있다.

$$\begin{align} \mathbf{F} &= m_1 \dfrac{d^2}{dt^2} (\mathbf{r}_1 - \mathbf{r}_{\mathrm{cm} })\\ &= \dfrac{m_1}{m_1 + m_2} \dfrac{d^2}{dt^2} \left[ (m_1 + m_2)\mathbf{r}_1 - (m_1 \mathbf{r}_1 + m_2 \mathbf{r}_2) \right]\\ &= \dfrac{m_1 m_2}{m_1 + m_2} \dfrac{d^2}{dt^2} (\mathbf{r}_1 - \mathbf{r}_2) \end{align}$$

또한 만약에 \(m_1 >> m_2\)일 경우 \(\mu \approx m_2\)이므로 가벼운 물체만 움직이는 것처럼 보일 것이다.

 

수소 원자에서 유효질량을 계산해보자. 양성자의 질량과 전자의 질량은 각각 다음과 같다.

$$ m_p = 1.672621868 \times 10^{-27}\ \mathrm{kg} $$

$$ m_e = 9.10938356 \times 10^{-31}\ \mathrm{kg} $$

유효질량을 계산하면 다음과 같다.

$$ \mu = 9.104425135 \times 10^{-31}\ \mathrm{kg} $$

양성자가 전자에 비해 훨씬 무겁기 때문에 거의 전자의 질량과 같은 것을 확인할 수 있다.

(casio fx-570es plus에 있는 값을 사용하였습니다.)

 

수소꼴 원자모형의 슈뢰딩거 방정식

원자 번호 \(Z\)인 임의의 수소꼴 원자의 원자핵에서 전자의 운동을 바라본다고 생각해보자. 무한대에서 포텐셜을 0이라고 할 때, 전자가 핵에 붙잡혀 있으므로 아래의 coulomb attraction에 의한 포텐셜이 생긴다. (\(e\)는 기본전하량, \( \varepsilon_0\)는 진공에서의 유전율, 쿨롱 상수가 \(k = \frac{1}{4\pi \varepsilon_0} \)의 관계가 있음을 상기하자)

$$ V(r) = - \dfrac{Ze^2}{4\pi \varepsilon_0 r} $$

그러므로 슈뢰딩거 방정식은

$$ -\dfrac{\hbar^2}{2\mu} \nabla^2 \psi + V(r) \psi = E \psi$$

이때 포텐셜 항이 반지름에만 의존하므로 구면좌표계에서 푸는 것이 유리할 것이다. 구면에서 라플라시안은 다음과 같다.

$$ \begin{align} \nabla^2 &= \dfrac{1}{r^2} \dfrac{\partial}{\partial r} \left(r^2 \dfrac{\partial}{\partial r} \right) + \dfrac{1}{r^2} \Lambda^2 \\ \Lambda^2 &= \dfrac{1}{\sin \theta} \left( \sin\theta \dfrac{\partial}{\partial \theta} \right) + \dfrac{1}{\sin^2 \theta} \dfrac{\partial^2}{\partial \phi^2} \end{align}$$

(이에 대한 유도는 5.델 연산자와 라플라시안의 좌표 변환에서 다루고 있다.) 이를 대입하면

$$ -\dfrac{\hbar^2}{2\mu r^2} \left[ \dfrac{\partial}{\partial r}\left(r^2 \dfrac{\partial \psi}{\partial r} \right) + \Lambda^2 \psi \right]  \psi + V(r) \psi = E \psi$$

\(\psi(r, \theta, \phi) = R(r)Y(\theta, \phi) \)로 변수 분리를 한뒤, \(RY\)로 나누면

$$ -\dfrac{\hbar^2}{2\mu r^2} \left[ \dfrac{1}{R} \dfrac{d}{d r}\left(r^2 \dfrac{d R}{d r} \right) + \dfrac{1}{Y} \Lambda^2 Y \right]  \psi + V(r) = E $$

지난 6.회전 운동(Rotational Motion)과 구면 조화 함수(Spherical Harmonic Function)에서 각에 대한 부분의 해가 구면 조화 함수이며, 아래와 같은 성질을 가지는 것을 증명하였다.

$$ \Lambda^2 Y_l^m = -l(l+1) Y_l^m$$

이를 대입하면 다음과 같이 반지름에 관한 미분 방정식을 얻을 수 있다.

$$ -\dfrac{\hbar^2}{2\mu r^2} \dfrac{d}{dr} \left(r^2 \dfrac{dR}{dr} \right) + V_{\mathrm{eff} }(r) R = ER $$

$$ V_{\mathrm{eff} } (r) = - \dfrac{Ze^2}{4\pi \varepsilon_0 r } + \dfrac{l(l+1) \hbar^2}{2\mu r^2} $$

*\( V_{\mathrm{eff} } (r) \)에서 두번째 항은 고전적인 개념에 빗대어 설명한다면 원심력으로 인해 감소하는 포텐셜이라고 해석할 수 있다.

\(V_{\text{eff} }\)의 개형

 

이때 반지름이 0일 때는 양성자와 부딪혀 전자가 쌍소멸(pair annihilation)을 일으킬 것이니 전자가 존재해선 안되며, 원자핵에 전자가 속박되어 있으므로 무한대에서도 전자가 존재해선 안된다. 따라서 아래와 같은 경계조건을 가진다.

$$ R(0) = \lim_{r \rightarrow \infty} R(r) = 0 $$

 

변수 치환하기

이 미분방정식에서 \( \dfrac{\partial}{\partial r}\left(r^2 \dfrac{\partial R}{\partial r} \right) \)안에 있는 \(r^2\)이 풀이를 어렵게하고 있다. 그래서 \(R = \frac{y}{r}\)이라고 치환을 한다면, 몫의 미분에서 \(r^2\)의 역수가 생겨 없어질 것이다.

$$\dfrac{d}{dr} \left(r^2 \dfrac{dR}{dr} \right) = \dfrac{\partial}{\partial r}\left(\dfrac{\partial y}{\partial r}y -y \right) = \dfrac{dy^2}{dr^2} $$

이를 방정식에 적용하고 \(V_{\mathrm{eff}}\)도 대입한 뒤 정리하면 다음과 같은 식을 얻을 수 있다.

$$ \dfrac{d^2y}{dr^2} + \left[ \dfrac{2\mu E}{\hbar^2} + \dfrac{\mu Z e^2}{2\pi \varepsilon_0 \hbar^2 r} - \dfrac{l(l+1)}{r^2} \right]y = 0 $$

보어 반지름 \(a_0 = \dfrac{4\pi \varepsilon_0 \hbar^2}{m_e e^2} \), \(a = \dfrac{m_e}{\mu}a_0\)를 이용하여 물리상수들을 정리하면 아래와 같다.

$$ \dfrac{d^2y}{dr^2} + \left[ \dfrac{2\mu E}{\hbar^2} + \dfrac{2Z }{a r} - \dfrac{l(l+1)}{r^2} \right]y = 0 $$

\(\kappa = \dfrac{ \sqrt{-2\mu E }}{\hbar}\), \( s = \kappa r \), \(\lambda = \dfrac{2Z}{a \kappa} \)로 치환하면 아래와 같이 깔끔한 식이 나온다.

$$ \dfrac{d^2y}{ds^2} + \left[ -1 + \dfrac{\lambda }{s} - \dfrac{l(l+1)}{s^2} \right]y = 0 $$

 

점근해 구하기

\(\kappa\)가 복소수이니(에너지는 양수이니) \(s\)도 복소수지만, 목적은 미분 방정식을 간단하게 만드는 것이니 실수처럼 생각하고 점근해를 구해보자

\(s \rightarrow \infty\)때

$$ \dfrac{d^2y}{ds^2}  -y = 0 $$

\(y = e^{\pm s}\)가 나오지만 발산하는 해를 버리고 \(y = e^{- s}\)를 택한다.

\(s \rightarrow 0\)때

$$ \dfrac{d^2y}{ds^2}  - \dfrac{l(l+1)}{s^2} y = 0 $$

$$ s^2 \dfrac{d^2y}{ds^2}  - l(l+1) y = 0 $$

이는 Euler-Cauchy Equation이므로 \(y = s^t\)이라고 둔다면 \(s^{-l}\), \(s^{l+1}\)이라는 해를 구할 수 있다. 이중 발산하는 해를 버리고 \(s^{l+1}\)택한다.

 

그럼 이제 \(y = q(s) L(s) \) , \(q(s) = s^{l+1} e^{-s} \)라고 두고 방정식에 집어넣으면 아래의 결과를 얻을 수 있다.

$$ q(s)\left[ \dfrac{d^2 L}{ds^2} + 2 \left( \dfrac{l+1}{s} -1 \right) \dfrac{dL}{ds} + \left( \dfrac{\lambda}{s} - \dfrac{2(l+1)}{s} \right)L \right] = 0 $$

\(q(s) \neq 0 \)이므로 (반지름이 0일수는 없으니까)

$$ \dfrac{d^2 L}{ds^2} + (2l + 2 -2s) \dfrac{dL}{ds} + [\lambda - 2(l+1)]L = 0 $$

\(2s\)가 불편하니 \(2s = \rho\)라고 치환하자, \(n = \frac{\lambda}{2} \)라 두면 아래와 같이 쓸 수 있다.

$$ s\dfrac{d^2 L}{ds^2} + [(2l + 1) + 1 -s] \dfrac{dL}{ds} + (n - l -1)L = 0 $$

이때 Associated Laguerre Equation은 아래와 같고(이 방정식의 풀이는 밑에서 하고 지금은 결론을 먼저 보자)

$$ x \dfrac{d^2 y}{dx^2} + (k+1-x) \dfrac{dy}{dx} + ny\ (n \geq 0, n \in \mathbb{Z}) $$

이 방정식의 해는 Associated Laguerre Function으로 다음과 같다. (해가 발산하지 않으려면 \(n\)이 음이 아닌 정수여야 한다.)

$$ L_n^{(k)} (x) = \sum_{m=0}^n (-1)^{m} \binom{n+k}{m+k} \dfrac{x^m}{m!} = \sum_{m=0}^n (-1)^{m} \dfrac{ (n+k)! }{ (n-m)! (m+k)! m! }x^m $$

지난  6.회전 운동(Rotational Motion)과 구면 조화 함수(Spherical Harmonic Function)에서 경계조건을 만족하기 위해서 \(l, m\)은 정수이며, \(m = 0 \pm 1, \pm 2, \cdots \), \(0 \leq |m| \leq l \) 이어야 했다.  그러므로 \(n - l -1 \geq 0\)에서 \(n \geq l+1 \)이라는 관계를 만족하여야 된다는 것을 알 수 있다.

 

이로부터 \(L = L_{n-l-1}^{(2l+1)} (\rho) \)임을 알 수 있으며, 이를 적용하여 \(R(r)\)을 쓰면 다음과 같다.

$$ R(r) = \dfrac{1}{r} s^{l+1} e^{-s} L_{n-l-1}^{(2l+1)} (\rho) $$

치환했던 변수들을 정리하자 \( \lambda = 2n\) , \(2s = \rho\), \(\kappa = \dfrac{ \sqrt{-2\mu E }}{\hbar}\), \( s = \kappa r \), \(\lambda = \dfrac{2Z}{a \kappa} \), \(R = \frac{y}{r} \), \(a_0 = \dfrac{4\pi \varepsilon_0 \hbar^2}{m_e e^2} \), \(a = \dfrac{m_e}{\mu} a_0 \)로부터 아래와 같이 식을 정리할 수 있다.

$$\rho = \dfrac{2Zr}{na} $$

$$ R(r) = \dfrac{2Z}{2^l na} \cdot \rho^l e^{-\frac{\rho}{2}} L_{n-l-1}^{(2l+1)} (\rho) $$

이때 정규화를 거처야되고, 정규화를 거치는 과정에서 상수배가 될 것이기 때문에 지금은 상수항을 날리고, 정규화 상수 \(N\)을 사용하여 표기하자

$$ R(r) = N \rho^l e^{-\frac{\rho}{2}} L_{n-l-1}^{(2l+1)} (\rho) $$

또한 아래와 같이 에너지도 얻을 수 있다.

$$ E = - \dfrac{Z^2\hbar^2}{2\mu a^2} \cdot \dfrac{1}{n^2} = - \dfrac{\mu Z^2 e^4}{32 \pi^2 \varepsilon^2 \hbar^2} \cdot \dfrac{1}{n^2}  $$

수소의 경우 에너지를 계산하면 다음과 같다. (고등학교 물리 교과서에도 나오는 유명한 공식이다.)

$$ E_n - \dfrac{13.6}{n^2}\ \mathrm{eV}$$

 

\(R(r) \)는 위와 같은 개형을 가진다.

 

그럼 이제 미분 방정식을 풀어보자. 저번 Legendre와 마찬가지로 Laguerre Equation을 먼저 살펴본 뒤, 이를 이용하여 Associated Lagurre Equation을 풀 것이다.

Mathmatical Background: Laguerre Equation

Laguerre Equation은 다음의 미분방정식이다.

$$ x \dfrac{d^2 y}{dx^2} + (1-x) \dfrac{dy}{dx} + ny = 0\ (n \geq 0, n \in \mathbb{Z}) $$

$$ \dfrac{d^2 y}{dx^2} + \dfrac{1-x}{x} \dfrac{dy}{dx} + \dfrac{nx}{x^2}y $$

에서 \(1-x\)와 \(nx\)가 \(x=0\)에서 해석적(analytic) 하기 때문에 프로베니우스의 해법(Frobenius's Method)을 적용할 수 있다.

$$ \begin{align} \text{let}\ y &= \sum_{m=0}^\infty a_m x^{m+r}\\ y' &= \sum_{m=0}^\infty (m+r)a_m x^{m+r-1}\\ y'' &= \sum_{m=0}^\infty (m+r)(m+r-1)a_{m} x^{m+r-2} \end{align} $$

을 방정식에 대입한 뒤 정리하면 다음과 같은식을 얻을 수 있다.

$$ a_0 r^2 x^{r-1} + \sum_{m=0}^\infty \left[ (m+r+1)^2 a_{m+1} + (n-m-r) a_m \right] x^{m+r} = 0 $$

\(a_0 = 0\)이면 trivial solution이 되므로, \(a_0 \neq 0 \)이라면 위 식이 성립하기 위해선 먼저 \(r^2 = 0\), 즉 \(r=0\)이라는 중근을 갖는다. 이를 대입하면 아래와 같은 점화식을 얻을 수 있다.

$$ a_{m+1} = - \dfrac{n-m-r}{(m+r+1)^2} a_m = - \dfrac{n-m}{(m+1)^2} a_m $$

이 경우에는 아래와 같이 점화식으로부터 쉽게 일반항을 얻을 수 있다.

$$ \begin{align} a_m &= - \dfrac{n-m+1}{m^2}a_{m-1}\\ &= (-1)^m \times \dfrac{n-m+1}{m^2} \times \dfrac{n-m+2}{(m-1)^2} \times \dfrac{n-m+3}{(m-2)^2} \times \cdots \times \dfrac{n-1}{2^2} \times \dfrac{n}{1^2} a_0 \\ &= (-1)^m \dfrac{n!}{ (n-m)! \cdot (m!)^2} a_0 \\ &= (-1)^m \binom{n}{m} \dfrac{a_0}{m!} \end{align} $$

참고로 \(\binom{n}{m}\)은 조합(conbination)의 표기법 중 하나이다.

\(a_0\)는 적분 상수에 해당하므로, 이부분을 때면 미분 방정식의 해는 다음과 같다.

$$ y = \sum_{m=0}^\infty (-1)^m \binom{n}{m} \dfrac{x^m}{m!} $$

그런데 \(n\)이 음이 아닌 정수이므로 점화식을 살펴보면 다항식으로 끝난다는 것을 알 수 있다. 그러므로 아래와 같이 Lagurre Function이 정의된다.

$$ L_n (x) = \sum_{m=0}^n (-1)^m \binom{n}{m} \dfrac{x^m}{m!} $$

참고로 \(n\)이 음이 아닌 정수일 때만 non-singular solution을 가진다는 것이 알려져 있다. 그리고 2계 미분 방정식이니 2개의 기저가 존재하지만, 중근이 나왔기 때문에 구한 것은 하나이다. 다만 물리적으로 의미가 있지는 않기 때문에 제 2의 해에 대해서는 생략한다.

https://en.wikipedia.org/wiki/Laguerre_polynomials#

Laguerre polynomials - Wikipedia

 

●Rodrigues's Fomular of Laguerre Function

Lagurre Function은 다음과 같이 쓸 수 있다.

$$ L_n (x) = \dfrac{e^x}{n!} \dfrac{d^n}{dx^n} (x^n e^{-x}) $$

 

눌러서 증명 펼치기

증명)

\(x^n e^x \)의 미분하였을 때, 계수가 이항 계수(binomial coefficient)임은 쉽게 떠올릴 수 있을 것이다. (잘 모르겠으면 고등학교 확률과 통계에 나오는 이항정리를 공부하고 오자) 여기서 \( (-1)^m \)이 있으니 \(x^n e^{-x} \)의 미분을 생각해보자. 이항 정리를 사용하여 곱미분을 전개하면 다음과 같다.

$$ \begin{align} \dfrac{d^n}{dx^n}(x^n e^{-x}) &= \sum_{m=0}^{n} \binom{n}{m} \left( \dfrac{d^m}{dx^m} e^{-x} \right) \left( \dfrac{d^{n-m} }{dx^{n-m}} x^n \right)\\ &= \sum_{m=0}^n \binom{n}{m} \cdot (-1)^m e^{-x} \cdot \dfrac{n!}{((n- (n-m))!}x^{n-(n-m)}\\ &= \sum_{m=0}^n \binom{n}{m} \cdot (-1)^m e^{-x} \cdot \dfrac{n!}{m!}x^{m}\\ &= n! e^{-x} L_n (x) \end{align} $$

$$\therefore L_n (x) = \dfrac{e^x}{n!} \dfrac{d^n}{dx^n} (x^n e^{-x}) $$

 

●Generating Function of Laguerre Function

Laguerre Function의 Generating Function은 다음과 같다.

$$ g(x, t) = \dfrac{1}{1-t} \exp \left(- \dfrac{xt}{1-t} \right) = \sum_{n=0}^{\infty} L_n (x) t^n $$

 

눌러서 증명 펼치기

증명)

Generating Function을 직접 찾아내는 것은 쉬운 일이 아니지만, 주어진 함수가 Generating Function인지 확인하는 것은 그보다 쉬운 일이다. \(t\)에 대해서 Taylor Expension을 하면 다음과 같다.

$$ g(x, t) = \sum_{n=0}^\infty \dfrac{t^n}{n!} \dfrac{\partial^n}{\partial t^n} \left[ \dfrac{1}{1-t} \exp \left( -\dfrac{xt}{1-t} \right) \right] \Biggr|_{t=0} $$

exp를 \(x\)에 대해 Taylor Expension을 하면

$$ \begin{align} g(x, t) &=  \sum_{n=0}^\infty \dfrac{t^n}{n!} \dfrac{\partial^n}{\partial t^n} \left[ \sum_{m=0}^\infty (-1)^m \dfrac{x^m}{m!} \cdot t^m (1-t)^{-m-1} \right] \Biggr|_{t=0}\\ &= \sum_{n=0}^\infty \dfrac{t^n}{n!} \sum_{m=0}^{\infty} (-1)^m \dfrac{x^m}{m!} \sum_{k=0}^n \binom{n}{k} \dfrac{d^k t^m}{dt^k} \Bigr|_{t=0} \cdot \dfrac{d^{n-k} }{dt^{n-k} } (1-t)^{-m-1} \Bigr|_{t=0} \end{align} $$

여기서 \(t\)에 \(0\)을 대입할 것이니 \(m=k\)때 \(t^m\)이 \(m\)번 미분된 경우만 살아남는다. 그러므로

$$ \begin{align} g(x, t) &= \sum_{n=0}^\infty \dfrac{t^n}{n!} \sum_{m=0}^{\infty} (-1)^m \dfrac{x^m}{m!} \binom{n}{m} m! \cdot \dfrac{n!}{m!}\\ &= \sum_{n=0}^\infty t^n \sum_{m=0}^\infty (-1)^m \binom{n}{m} \dfrac{x^m}{m!}\\ &= \sum_{n=0}^\infty L_n (x) t^n \end{align} $$

임을 알 수 있다.

 

Mathmatical:Background Associated Laguerre Function

Associated Laguerre Equation은 다음의 미분방정식이다.

$$ x \dfrac{d^2 y}{dx^2} + (k+1-x) \dfrac{dy}{dx} + ny=0\ (n \geq 0, n \in \mathbb{Z}) $$

Lagurre Equation \( x \dfrac{d^2 y}{dx^2} + (1-x) \dfrac{dy}{dx} + ny \)을 한 번 미분하면 다음과 같다.

$$ x \dfrac{d^3 y}{dx^3} + (2-x) \dfrac{dy^2}{dx^2} + (n-1)\dfrac{dy}{dx} $$

여기에 착안하여 \( x \dfrac{d^2 y}{dx^2} + (1-x) \dfrac{dy}{dx} + (n+k)y \)를 \(k\)번 미분하면 다음과 같다. (곱미분에 이항정리를 적용한 Leibnitz Rule을 적용하여 미분하면 된다. 비슷한 과정을 앞에서 여러번 반복하였으니 과정은 생략한다.)

$$ x y^{(k+2)} + (k+1-x) y^{(k+1)} + ny^{(k)} $$

이는 Laguerre Equation이므로 해는 다음과 같이 Laguerre Function을 이용하여 쓸 수 있다.

$$ y = L_n^{(k)} = (-1)^k \dfrac{d^k}{dx^k} L_{n+k} $$

*여기서 왜 \( (-1)^k\)이 붙는지 의문을 가질 수 있는데, 이는 부호를 Lagurre Function과 맞추어 식 전개를 간편하게 하기 위함이다.

$$\begin{align} \dfrac{d^k}{dx^k} L_{n+k} (x) &= \dfrac{d^k}{dx^k} \sum_{m=0}^n (-1)^m \binom{n+k}{m} \dfrac{1}{m!} x^m\\ &= \sum_{m=k}^n (-1)^m \binom{n+k}{m} \dfrac{x^{m-k}}{ (m-1)!}\\ &= (-1)^k \sum_{m=0}^n (-1)^{m} \binom{n+k}{m+k} \dfrac{x^m}{m!} \end{align} $$

이므로 Associated Laguerre Function은 다음과 같다.

$$ L_n^{(k)} (x) = \sum_{m=0}^n (-1)^{m} \binom{n+k}{m+k} \dfrac{x^m}{m!} = \sum_{m=0}^n (-1)^{m} \dfrac{ (n+k)! }{ (n-m)! (m+k)! m! }x^m $$

●Rodrigues's Formular of Associated Laguerre Function

Laguerre Function에 대한 Rodrigues's Formular를 유도할때와 같은 방법으로 다음의 공식을 얻을 수 있다.

$$ L_n^{(k)} = \dfrac{x^{-k} e^x}{n!} \dfrac{d^k}{d x^k} (x^{n+k} e^{-x}) $$

아니면 Reibnitz Rule로 위 식을 전개해봄으로써 확인해볼 수 도 있다.

●Generating Function of Associated Laguerre Function

Laguerre Function의 generating function에서

$$ \dfrac{1}{1-t} \exp \left(- \dfrac{xt}{1-t} \right) = \sum_{n=0}^{\infty} L_{n+k} (x) t^{n+k} $$

이를 \(x\)에 대해 \(k\)번 미분하면 아래 식을 얻는다.

$$ \dfrac{(-t)^k}{(1-t)^{k+1}} \exp \left(- \dfrac{xt}{1-t} \right) = \sum_{n=0}^{\infty} (-1)^k L_{n}^{(k)}(x) t^{n+k} $$

양변을 \((-t)^k \)로 나눠주면 아래와 같은 Generating Function을 얻는다.

$$ \dfrac{1}{(1-t)^{k+1}} \exp \left(- \dfrac{xt}{1-t} \right) = \sum_{n=0}^{\infty} L_{n}^{(k)}(x) t^{n} $$

●Recurrence Relation of Associated Laguerre Function

$$(n+1)L_{n+1}^{(k)} = (2n + k + 1 - x) L_n^{(k)} - (n+k) L_{n-1}^{(k)} $$

밑에서 Nomailization할 때 필요한 관계식이다.

 

눌러서 증명 펼치기

증명)

Generating function을 \(t\)에 대해 미분하면 다음 식을 얻을 수 있다.

$$ \left( \dfrac{k+1}{1-t} - \dfrac{x}{(1-t)^2} \right) \sum_{n=0}^\infty L_n^{(k)} t^n = \sum_{n=0}^\infty nL_n^{(k)} t^{n-1} $$

이는 \(t\)에 대한 항등식이니 양변에 \( (1-t)^2 \)을 곱한뒤 \(t^m\)항만 적으면 다음과 같다.

$$(k+1) (L_m^{(k)} - L_{m-1}^{(k)} ) - x L_m^{(k)} = (m+1) L_{m+1}^{(k)} - 2m L_m^{(k)} + (m-1)L_{m-1}^{(k)} $$

이를 정리한 뒤, \(m\)을 \(n\)으로 바꾸면 위의 관계식을 얻을 수 있다.

●Orthogonality of Associated Laguerre Function

Associated Laguerre Function은 가중함수 \(x^k e^{-x} \)에 대하여 다음과 같은 직교성을 가진다.

$$ \int_0^{\infty} x^k e^{-x} L_n^{(k)} L_m^{(k)} dx = \dfrac{(n+m)!}{n!} \delta_{n, m} $$

여기서 \(\delta_{m, n}\)은 Kronecker's delta이다.

 

눌러서 증명 펼치기

증명)

$$(x^{k+1} e^{-x} )' = (k+1-x)x^k e^{-x} $$

에서 Associated Laguerre Equation에  가중함수 \(x^k e^{-x} \)를 곱해주면 다음과 같이 쓸 수 있다.

$$\dfrac{d}{dx}\left[ x^{k+1} e^{-x} \dfrac{dy}{dx} \right] + nx^k e^{-x} y = 0 $$

그리고 경계조건 \( R(0) = \lim_{r \rightarrow \infty} R(r) = 0 \)을 가지므로 Strum-Liouvile problem이며, Strum-Liouvile problem의 해는 직교함으로 Associated Laguerre Function은 가중함수 \(x^k e^{-x} \)에 대하여 직교성을 가지고 있음을 알 수 있다. (Strum-Liouvile problem에 대해서는 2.슈뢰딩거 방정식(Schrodinger Equation)과 해의 직교성에서 다루고 있다.) 그러므로 \(n=m\)때를 살펴보자

Rodrigues's Formular \( L_n^{(k)} = \dfrac{x^{-k} e^x}{n!} \dfrac{d^k}{d x^k} (x^{n+k} e^{-x}) \)를 이용하면 다음과 같이 식을 전개할 수 있다.

 

$$ \int_0^{\infty} x^k e^{-x} \left[L_n^{(k)} \right]^2 dx = \dfrac{1}{n!} \int_0^\infty L_n^{(k)} \dfrac{d^n}{d x^n}(x^{n+k} e^{-x}) dx $$

\(L_n^{(k)}\)가 최고차항의 계수가 \( (-1)^m \)인 \(n\)차 다항식이라는 것을 상기하여 부분적분을 해보자. 대괄호 안에 있는 식들은 \(x\)와 \(e^{-x}\)가 남아있어 0과 \(\infty\)에서 0이 된다. 그리고 같은 방법으로 부분적분을 한 번 더 해보면 아래와 같이 계산된다.

$$\begin{align} \int_0^{\infty} x^k e^{-x} \left[L_n^{(k)} \right]^2 dx &= \dfrac{(-1)^n }{n!} \int_0^\infty (-1)^n x^{n+k} e^{-x} dx \\ &= \dfrac{1}{n!} \int_0^\infty x^{n+k} e^{-x} dx \\ &= \dfrac{(-1)^{n+k}}{n!} \int_0^\infty (n+k)! \cdot (-1)^{n+k} e^{-x} dx\\ &= \dfrac{(n+k)!}{n!} \int_0^{\infty} e^{-x} dx\\ &= \dfrac{(n+k)!}{n!} \end{align} $$

 

Nomalization

지금가지 Laguerre Function의 Rodrigues's formular와 generating function, 직교성을 살펴본 것은 모두 정규화를 위해서 였다. 구한 수소꼴 원자의 파동함수에 대해서 정규화를 진행해보자. 부피에 대해서 정규화를 진행한다면

$$ \int_V |\psi(r, \theta, \phi)|^2 dV = 1 $$

\(dV = r^2 \sin\theta dr d\theta d\phi = r^2 d\Omega\)에서 ( \(\Omega\)는 입체각이다.)

$$ \int_0^\infty r^2|R(r)|^2 dr \int_\Omega |Y_l^m (\theta, \phi)|^2 d\Omega = 1 $$

구면 조화 함수 \( Y_m^l \)에 대해서는 지난 글에서 정규화를 진행하였으므로 아래 식만 계산하면 된다.

$$ \int_0^\infty r^2|R(r)|^2 dr  = 1 $$

$$ \begin{align} \int_0^\infty r^2|R(r)|^2 dr  &= N^2 \int_0^\infty r^2 \rho^{2l} e^{-\rho} \left[ L_{n-l-1}^{(2l+1)}(\rho)\right]^2 dr\\ &= \left( \dfrac{na}{2Z} \right)^3 N^2 \int_0^\infty \rho^{(2l+1) + 1} e^{-\rho} \left[ L_{n-l-1}^{(2l+1)} \right]^2\\ &= \left( \dfrac{na}{2Z} \right)^3 N^2 \cdot \dfrac{(n+l)!}{(n-l-1)!} \cdot 2n  \end{align} $$

$$ \therefore N = \sqrt{ \left( \dfrac{2Z}{na} \right)^3 \dfrac{(n-l-1)!}{2n(n+l)!} } $$

 

여기서

$$ \int_0^{\infty} x^{k+1} e^{-x} \left[L_n^{(k)} \right]^2 dx $$

형태의 적분은 Recurrence Relation \( (n+1)L_{n+1}^{(k)} = (2n + k + 1 - x) L_n^{(k)} - (n+k) L_{n-1}^{(k)} \)을 이용하면

$$ x L_n^{(k)} = (2n+k+1) L_n^{(k)} - (n+1) L_{n+1}^{(k)} - (n+k) L_{n-1}^{(k)} $$

그럼 직교성 때문에 아래와 같이 식이 전개된다.

$$ \begin{align} \int_0^{\infty} x^{k+1} e^{-x} \left[L_n^{(k)} \right]^2 dx &= (2n+k+1) \int_0^\infty x^k e^{-x} \left[L_n^{(k)} \right]^2 dx - (n+1) \int_0^\infty x^k e^{-x} L_n^{(k)} L_{n+1}^{(k)} dx\\ &-(n+k) \int_0^\infty x^k e^{-x} L_n^{(k)} L_{n-1}^{(k)} dx\\ &= (2n+k+1) \int_0^\infty x^k e^{-x} \left[ L_n^{(k)} \right]^2 dx\\ &= (2n+k+1) \dfrac{(n+k)!}{n!} \end{align} $$

 

그러므로 최종적인 파동함수는 아래와 같다.

$$ \psi_{n, l, m}(r, \theta, \phi) = \sqrt{ \left( \dfrac{2Z}{na} \right)^3 \dfrac{(n-l-1)!}{2n(n+l)!} } \rho^l e^{-\frac{\rho}{2}} L_{n-l-1}^{(2l+1)} (\rho) Y_m^{l} (\theta, \phi) $$

\(\rho = \dfrac{2Zr}{na} \), 보어 반지름 \(a_0 =  \dfrac{4\pi \varepsilon_0 \hbar^2}{m_e e^2} \), \(a = \dfrac{m_e}{\mu}a_0\)

 

Orbital Visualization and Real Orbital

 

https://energywavetheory.com/atoms/orbital-shapes/

오비탈이 위의 그림처럼 생겼다라고 고등학교 화학시간이나 일반화학 등에서 배웠을 것이다. 이를 한 번 시각화 해보자. 그런데 파동함수를 제곱하여 나타나는 확률밀도함수가 4차원 함수라서 그대로 시각화할 수 없다는 문제가 있다. 이를 해결하는 한가지 방법은 아래와 같이 점들로 공간을 가득 채운 뒤, 확률밀도값에따라 점의 투명도인 알파값(1이면 불투명, 0이면 투명)을 조절하는 방법으로 오비탈을 모양을 보이는 것이다.

 

 

import matplotlib.pyplot as plt

fig = plt.figure()
ax = fig.add_subplot(projection="3d")
ax.set_xlabel("x")
ax.set_ylabel("y")

X, Y, Z = np.meshgrid(np.linspace(-1, 1, 30), np.linspace(-1, 1, 30), np.linspace(-1, 1, 30))
scatter = ax.scatter(X, Y, Z)

점의 edgecolor와 그림자때문에 제대로 시각화가 어려우니 아래와 같은 옵션을 추가하여야 한다.

점의 edgecolor와 그림자 제거

fig = plt.figure()
ax = fig.add_subplot(projection="3d")
ax.set_xlabel("x")
ax.set_ylabel("y")

X, Y, Z = np.meshgrid(np.linspace(-1, 1, 30), np.linspace(-1, 1, 30), np.linspace(-1, 1, 30))
scatter = ax.scatter(X, Y, Z, edgecolor="none", depthshade=False)

그리고 파동함수를 파이썬 코드로 옮기자. 이때 보어 반지름 \(a_0 = 53\ \mathrm{pm} \)이 너무 작은 값이기 때문에 값을 그대로 사용하면 시각화하기 어렵다. 그래서 \(r' = \frac{Z}{a} \)로 무차원화하여 파동함수를 나타내면

$$ \psi_{n, l, m}(r', \theta, \phi) = \sqrt{ \left( \dfrac{2}{n} \right)^3 \dfrac{(n-l-1)!}{2n(n+l)!} } \rho^l e^{-\frac{\rho}{2}} L_{n-l-1}^{(2l+1)} (\rho) Y_m^{l} (\theta, \phi),\ \rho = \dfrac{2r'}{n} $$

아래를 이를 코드로 옮긴 것이다.

import math
import numpy as np
import scipy.special as spe
import matplotlib.pyplot as plt

#radial function
def R(n, l, r):
    coeff = np.sqrt( (2.0/n)**3 * math.factorial(n-l-1)/(2.0*n*math.factorial(n+l)) )
    
    #associated laguerre
    #https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.assoc_laguerre.html
    laguerre = spe.assoc_laguerre(2.0*r/n, n-l-1, 2.0*l+1)
    return coeff * (2.0*r/n)**l * np.exp(-r/n) * laguerre

#wavefunction
def psi(n:int, l:int, m:int, r:float, elev:float, azim:float)->complex:
    """Hydrogenic Wave Function

    Args:
        n (int): principal quantum number
        l (int): azimutal quantum number
        m (int): magnetic quantum number
        r (float): radius
        elev (float): elevation angle
        azim (float): azimutal angle

    Returns:
        complex: value
    """
    #spherical harmonic
    #https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.sph_harm.html
    return R(n, l, r) * spe.sph_harm(m, l, azim, elev)

def probd(f):
    return np.abs(f)**2

#spherical coordinate to cartersian coordinate 
def sph2cart_probd(wavefunc, n:int, l:int, m:int, x:float, y:float, z:float)->float:
    #arctan가 아닌 arctan2를 사용해야 회전 방향을 알 수 있음
    r = np.sqrt(x**2 + y**2 + z**2)
    elev = np.arctan2(np.sqrt(x**2 + y**2), z)
    azim = np.arctan2(y, x)
    return probd(wavefunc(n, l, m, r, elev, azim))

 

또한 랜더링할 반경은 오비탈의 반경에 대한 확률밀도함수 \(|rR(r)|^2 \)의 적분값이 0.999인 반지름으로 정하였다. 수치 적분을 통해 확률밀도함수의 적분을 구하고, 뉴튼-랩슨 방법등의 해를 찾는 수치해석 알고리즘을 파이썬에서 직접 구현해볼 수도 있지만(그리 어렵지 않다.), 아래와 같이 FORTRAN 라이브러리를 warping한 scipy에서 제공해주는 함수를 이용하는 것이 훨씬 성능이 좋다. 

from scipy.optimize import fsolve
import scipy.integrate 
def get_render_radius(n, l):
    #확률 밀도가 0.999인 반지름 구하기 
    f = lambda x : scipy.integrate.quad(lambda r : probd(r*R(n, l, r)), 0, x)[0] - 0.999
    return fsolve(f, 3*n)[0]

그다음에 아래와 같이 시각화를 해보자

render_r = get_render_radius(n, l)
resol = 70
X, Y, Z = np.meshgrid(np.linspace(-render_r, render_r, resol), np.linspace(-render_r, render_r, resol), np.linspace(-render_r, render_r, resol))
P = sph2cart_probd(psi, n, l, m, X, Y, Z)
    
fig = plt.figure(dpi=500)
ax = fig.add_subplot(projection="3d")
ax.set_title(f"Electron Cloud of ({n}, {l}, {m}) Hydrogenic Orbital")
ax.set_xlabel(r"$\frac{Z}{a}x$")
ax.set_ylabel(r"$\frac{Z}{a}y$")
ax.set_zlabel(r"$\frac{Z}{a}z$")

#making box
scatter = ax.scatter3D(X, Y, Z, edgecolors="none", depthshade=False)

#coloring
wa = P.flatten()/P.max()
scatter.set_facecolor([[0.7, 0.3, 0.6, a] for a in wa])
plt.subplots_adjust(right=0.5)

*set_facecolor method에 0~1범위의 RGBA값을 담은 1차원 리스트를 전달하면 전체 점의 색이 바뀌고, 2차원으로 전달하면 점별로 설정을 할 수 있다.

명암 없이 단일 색상으로만 칠하니 상당히 보기 밋밋하다. \(x, y, z\)값에 따라 색상이 약간씩 바뀌는 명암을 주자. 조명이 카메라 윗쪽에 있다고 생각한다면, \(z\)값이 크면 조금더 색깔이 밝아질 것이고, \(x\)값이 커도 밝을 것이며, \(y\)는 작으면 밝아질 것이다. 이를 이용하여 적당히 색상에 가중치를 주어 명암을 표현하면 다음과 같다.

 

 

def nomaliztion(x):
    temp = x + abs(x.min())
    return temp/temp.max()

def render_3d_cloud(n, l, m, mode, resol=70):
    if mode == "complex":
        f = psi
    elif mode == "real":
        f = psi_real
    else:
        raise ValueError("Unsurported mode")
    
    render_r = get_render_radius(n, l)
    X, Y, Z = np.meshgrid(np.linspace(-render_r, render_r, resol), np.linspace(-render_r, render_r, resol), np.linspace(-render_r, render_r, resol))
    P = sph2cart_probd(f, n, l, m, X, Y, Z)
    
    fig = plt.figure(dpi=500)
    ax = fig.add_subplot(projection="3d")
    ax.set_title(f"Electron Cloud of ({n}, {l}, {m}) Hydrogenic Orbital ({mode.capitalize()})")
    ax.set_xlabel(r"$\frac{Z}{a}x$")
    ax.set_ylabel(r"$\frac{Z}{a}y$")
    ax.set_zlabel(r"$\frac{Z}{a}z$")

    #making box
    scatter = ax.scatter3D(X, Y, Z, edgecolors="none", depthshade=False)
    
    #flatten
    xf = X.flatten()
    yf = Y.flatten()
    zf = Z.flatten()

    #RGBA weight
    wr = 0.7 + 0.1*zf/abs(zf).max() + 0.1*xf/abs(xf).max() - 0.1*yf/abs(yf).max()
    wg = 0.3 + 0.1*zf/abs(zf).max() + 0.1*xf/abs(xf).max() - 0.1*yf/abs(yf).max()
    wb = 0.6 + 0.1*zf/abs(zf).max() + 0.1*xf/abs(xf).max() - 0.1*yf/abs(yf).max()
    wa = P.flatten()/P.max()

    wr = nomaliztion(wr)
    wg = nomaliztion(wg)
    wb = nomaliztion(wb)
    
    #coloring
    color_map = np.array([wr, wg, wb, wa]).T
    scatter.set_facecolor(color_map)
    plt.subplots_adjust(right=0.5)

 

2s와 2p 오비탈들도 시각화해보면 아래와 같다.

2s complex orbital

2p orbital

2p orbital

2p orbital

 

알고 있었던 오비탈의 모양이랑 다를 것이다. 고등학교 화학이나 일반화학등에서 보여주는 오비탈은 Real Orbital로 Complex Orbital을 선형 결합하여 파동함수에서 허수부를 소거한 것이다. 오비탈에서 복소수가 있는 부분은 \(e^{im\phi}\)이므로, \(m\)의 크기가 같고, 부호가 반대인 파동함수는 서로 켤래(conjugate) 관계임을 알 수 있다. 따라서 실수인 p오비탈들은 다음과 같이 구한다.

$$ p_0 = \psi_{n, 1, 0} $$

$$ p_{1} = \dfrac{1}{\sqrt{2}} (\psi_{n, 1, 1} + \psi_{n, 1, -1}) = \sqrt{2} \Re (\psi_{n, 1, 1})  $$

$$ p_{-1} = \dfrac{1}{i\sqrt{2} } (\psi_{n, 1, 1} - \psi_{n, 1, -1}) = \sqrt{2} \Im (\psi_{n, 1, 1})  $$

*사람에 따라 \(p_{1}, p_{-1} \)을 반대로 정의할 수도 있다.

앞에 \(1/\sqrt{2}\)가 붙는 이유는 nomailiztion 때문으로

$$ \int |\psi + \psi^*|^2 d\tau = \int \psi \psi d\tau + \int \psi^* \psi^* d\tau + 2 $$

인데 복소함수의 성질에 의해서

$$ \int_0^{2 \pi} e^{\pm 2 i m \phi} d \phi = 0 $$

이기 때문이다. 그러므로 아래와 같이 코드를 추가해서 2s와 2p 오비탈들을 시각화 해보면 다음과 같다.

def psi_real(n:int, l:int, m:int, r:int, elev:float, azim:float)->complex:
    """Hydrogenic Wave Function (Real)

    Args:
        n (int): principal quantum number
        l (int): azimutal quantum number
        m (int): magnetic quantum number
        r (float): radius
        elev (float): elevation angle
        azim (float): azimutal angle

    Returns:
        complex: value
    """
    
    if m == 0:
        return psi(n, l, m, r, elev, azim)
    elif m > 0:
        return np.sqrt(2)*np.real(psi(n, l, abs(m), r, elev, azim))
    elif m < 0:
        return np.sqrt(2)*np.imag(psi(n, l, abs(m), r, elev, azim))

 

위에 있는 코드를 모으면 아래와 같다.

파동함수 관련

import math
import numpy as np
import scipy.special as spe

#radial function
def R(n, l, r):
    coeff = np.sqrt( (2.0/n)**3 * math.factorial(n-l-1)/(2.0*n*math.factorial(n+l)) )
    
    #associated laguerre
    #https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.assoc_laguerre.html
    laguerre = spe.assoc_laguerre(2.0*r/n, n-l-1, 2.0*l+1)
    return coeff * (2.0*r/n)**l * np.exp(-r/n) * laguerre

#wavefunction
def psi(n:int, l:int, m:int, r:float, elev:float, azim:float)->complex:
    """Hydrogenic Wave Function

    Args:
        n (int): principal quantum number
        l (int): azimutal quantum number
        m (int): magnetic quantum number
        r (float): radius
        elev (float): elevation angle
        azim (float): azimutal angle

    Returns:
        complex: value
    """
    #spherical harmonic
    #https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.sph_harm.html
    return R(n, l, r) * spe.sph_harm(m, l, azim, elev)

def psi_real(n:int, l:int, m:int, r:int, elev:float, azim:float)->complex:
    """Hydrogenic Wave Function (Real)

    Args:
        n (int): principal quantum number
        l (int): azimutal quantum number
        m (int): magnetic quantum number
        r (float): radius
        elev (float): elevation angle
        azim (float): azimutal angle

    Returns:
        complex: value
    """
    
    if m == 0:
        return psi(n, l, m, r, elev, azim)
    elif m > 0:
        return np.sqrt(2)*np.real(psi(n, l, abs(m), r, elev, azim))
    elif m < 0:
        return np.sqrt(2)*np.imag(psi(n, l, abs(m), r, elev, azim))

def probd(f):
    return np.abs(f)**2

#spherical coordinate to cartersian coordinate 
def sph2cart_probd(wavefunc, n:int, l:int, m:int, x:float, y:float, z:float)->float:
    #arctan가 아닌 arctan2를 사용해야 회전 방향을 알 수 있음
    r = np.sqrt(x**2 + y**2 + z**2)
    elev = np.arctan2(np.sqrt(x**2 + y**2), z)
    azim = np.arctan2(y, x)
    return probd(wavefunc(n, l, m, r, elev, azim))

시각화 관련

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve
import scipy.integrate 
def get_render_radius(n, l):
    #확률 밀도가 0.999인 반지름 구하기 
    f = lambda x : scipy.integrate.quad(lambda r : probd(r*R(n, l, r)), 0, x)[0] - 0.999
    return fsolve(f, 3*n)[0]

def nomaliztion(x):
    temp = x + abs(x.min())
    return temp/temp.max()

def render_3d_cloud(n, l, m, mode, resol=70):
    if mode == "complex":
        f = psi
    elif mode == "real":
        f = psi_real
    else:
        raise ValueError("Unsurported mode")
    
    render_r = get_render_radius(n, l)
    X, Y, Z = np.meshgrid(np.linspace(-render_r, render_r, resol), np.linspace(-render_r, render_r, resol), np.linspace(-render_r, render_r, resol))
    P = sph2cart_probd(f, n, l, m, X, Y, Z)
    
    fig = plt.figure(dpi=500)
    ax = fig.add_subplot(projection="3d")
    ax.set_title(f"Electron Cloud of ({n}, {l}, {m}) Hydrogenic Orbital ({mode.capitalize()})")
    ax.set_xlabel(r"$\frac{Z}{a}x$")
    ax.set_ylabel(r"$\frac{Z}{a}y$")
    ax.set_zlabel(r"$\frac{Z}{a}z$")

    #making box
    scatter = ax.scatter3D(X, Y, Z, edgecolors="none", depthshade=False)
    
    #flatten
    xf = X.flatten()
    yf = Y.flatten()
    zf = Z.flatten()

    #RGBA weight
    wr = 0.7 + 0.1*zf/abs(zf).max() + 0.1*xf/abs(xf).max() - 0.1*yf/abs(yf).max()
    wg = 0.3 + 0.1*zf/abs(zf).max() + 0.1*xf/abs(xf).max() - 0.1*yf/abs(yf).max()
    wb = 0.6 + 0.1*zf/abs(zf).max() + 0.1*xf/abs(xf).max() - 0.1*yf/abs(yf).max()
    wa = P.flatten()/P.max()

    wr = nomaliztion(wr)
    wg = nomaliztion(wg)
    wb = nomaliztion(wb)
    
    #coloring
    color_map = np.array([wr, wg, wb, wa]).T
    scatter.set_facecolor(color_map)
    plt.subplots_adjust(right=0.5)

 

참고: 다전자 원자의 오비탈

다전자 원자의 경우 삼체 이상이기 때문에 그냥은 해석적인 풀이가 불가능하다. 그러므로 아래와 같이 각 전자간의 상호작용을 무시하여, 각 궤도별 파동함수의 곱으로 표현하는 Orbital Approximation을 적용한다.

$$\psi(\mathbf{r}_1, \mathbf{r}_2, \cdots) = \psi(\mathbf{r}_1) \psi(\mathbf{r}_2) \cdots $$

여기서 각 궤도별 파동함수는 우리가 앞에서 구한 수소꼴 원자의 오비탈이다. 그리고 유효 핵전하(effective charge)를 사용하여 전자간의 상호작용(반발력)으로 인한 구심력의 감소를 핵의 전하가 줄어든 효과로 치환하여 표현한다.

여기서 유효 핵전하는 주양자수 뿐만이 아니라 오비탈의 모양을 결정하는 부양자수의 영향을 받기 때문에 다전자 원자의 에너지는 수소원자와 달리 주양자수와 부양자수의 영향을 받게 된다.

 

참고2: 상대론적 효과

우리가 구한 오비탈은 원자 번호가 작을때는 현상과 잘 맞지만, 원자 번호가 커지면 잘 맞지 않게된다. 원자 번호가 커질 수록, 핵의 인력으로 인한 구심력이 강해지기 때문에 전자의 속력이 커진다. 그런데 특수 상대성 이론에 따르면, 운동하는 속도가 광속에 가까울수록 질량이 커지는 효과를 내기 때문에 전자가 우리가 계산한 것보다 조금 안쪽 궤도를 돌게된다. 이를 보여주는 대표적인 원소가 6주기의 금과 수은이다. 주기율표에서 77번 이리듐, 78번 백금, 79번 금은 순서대로 d 오비탈에 전자가 하나씩 늘어난다. 하지만 이리듐과 백금은 여타 금속과 같이 자외선 영역의 전자기파를 흡수하고, 가시 광선 영역의 전자기파를 반사하여 은백색의 광택을 가지는 반면, 금의 경우 에너지 준위가 내려와서 푸른색의 전자기파까지 흡수하기 때문에 우리가 알고 있는 금색을 띈다. 그리고 금 옆에있는 수은은 인접한 원소를 보았을 때 상온에서 고체로 예측할 수 있지만, 실제로는 상온에서 액체이다.

Summary

지난 구면 조화 함수의 내용까지 포함되었다.

 

슈뢰딩거 방정식

$$ -\dfrac{\hbar^2}{2\mu r^2} \nabla^2 \psi + V(r) \psi = E \psi$$

$$ V(r) = - \dfrac{Ze^2}{4\pi \varepsilon_0 r} $$

$$ \begin{align} \nabla^2 &= \dfrac{1}{r^2} \dfrac{\partial}{\partial r} \left(r^2 \dfrac{\partial}{\partial r} \right) + \dfrac{1}{r^2} \Lambda^2 \\ \Lambda^2 &= \dfrac{1}{\sin \theta} \left( \sin\theta \dfrac{\partial}{\partial \theta} \right) + \dfrac{1}{\sin^2 \theta} \dfrac{\partial^2}{\partial \phi^2} \end{align}$$

 

경계조건

$$ \psi(0 , \theta, \phi) = \lim_{r \rightarrow \infty} \psi(r, \theta, \phi) = 0 $$

$$ \psi(r, 0, \phi) = \psi(r, 2\pi, \phi ) $$

$$ \psi(r, \theta, 0) = \psi(r, \theta, 2 \pi) $$

 

각에 대한 미분 방정식의 해: 구면 조화 함수

정규화된 구면 조화 함수 ( \( 0 \leq |m| \leq l \), \(m = 0, \pm 1, \pm 2, \cdots \) )

$$ Y_l^m (\theta, \phi) =  \sqrt{ \dfrac{(2l+1) (l-|m|)!}{4\pi (l+ |m|)!} } P_l^m ( \cos\theta ) e^{im\phi} $$

제 1종 버금 르장드르 함수

$$ P_l^m (x) = (1-x^2)^{\frac{|m|}{2} } \dfrac{d^{|m|} P_l}{dx^{|m|} } = \dfrac{(1-x^2)^{\frac{|m|}{2}}}{2^l \cdot l!} \dfrac{d^{(l+|m|)} }{dx^{(l+|m|)} } (x^2-1)^l $$

 

반지름에 대한 미분방정식

$$ -\dfrac{\hbar^2}{2\mu r^2} \dfrac{d}{dr} \left(r^2 \dfrac{dR}{dr} \right) + V_{\mathrm{eff} }(r) R = ER $$

유효 포텐셜

$$ V_{\mathrm{eff} } (r) = - \dfrac{Ze^2}{4\pi \varepsilon_0 r } + \dfrac{l(l+1) \hbar^2}{2\mu r^2} $$

 

수소꼴 원자의 파동함수 \(n \geq l +1 \)

$$ \psi_{n, l, m}(r, \theta, \phi) = \sqrt{ \left( \dfrac{na}{2Z} \right)^3 \dfrac{(n-l-1)!}{2n(n+l)!} } \rho^l e^{-\frac{\rho}{2}} L_{n-l-1}^{(2l+1)} (\rho) Y_m^{l} (\theta, \phi) $$

\(\rho = \dfrac{2Zr}{na} \), 보어 반지름 \(a_0 =  \dfrac{4\pi \varepsilon_0 \hbar^2}{m_e e^2} \), \(a = \dfrac{m_e}{\mu}a_0\)

 

Associated Laguerre Function

$$ L_n^{(k)} (x) = \sum_{m=0}^n (-1)^{m} \binom{n+k}{m+k} \dfrac{x^m}{m!} = \sum_{m=0}^n (-1)^{m} \dfrac{ (n+k)! }{ (n-m)! (m+k)! m! }x^m = \dfrac{x^{-k} e^{x} }{n!} \dfrac{d^n}{dx^n} (x^{n+k} e^{-x} ) $$

 

에너지

$$ E = - \dfrac{Z^2\hbar^2}{2\mu a^2} \cdot \dfrac{1}{n^2} = - \dfrac{\mu Z^2 e^4}{32 \pi^2 \varepsilon^2 \hbar^2} \cdot \dfrac{1}{n^2}  $$

 

수소일 경우 에너지

$$ E_n - \dfrac{13.6}{n^2}\ \mathrm{eV}$$

 

각운동량의 크기

$$ |L|^2 = l(l+1) \hbar^2 $$

 

각운동량의 \(z\)축 성분

$$ \hat{L}_z \psi_{n, l, m} (r, \theta, \phi) =  m\hbar \psi_{n, l, m} (\theta, \phi) $$

 

그러므로 주양자수 \(n\)은 전자의 에너지를, 부 양자수 \(l\)은 전자의 각운동량을, 자기 양자수 \(m\)은 각운동량의 \(z\) 성분을 결정한다.

 

Real Orbital

$$ p_0 = \psi_{n, 1, 0} $$

$$ p_{1} = \dfrac{1}{\sqrt{2}} (\psi_{n, 1, 1} + \psi_{n, 1, -1}) = \sqrt{2} \Re (\psi_{n, 1, 1})  $$

$$ p_{-1} = \dfrac{1}{i\sqrt{2} } (\psi_{n, 1, 1} - \psi_{n, 1, -1}) = \sqrt{2} \Im (\psi_{n, 1, 1})  $$

 

참고문헌

Atkins' Physical Chemistry 11e

Erwin Keyszig Advanced Engineering 10th

https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.assoc_laguerre.html

assoc_laguerre — SciPy v1.14.0 Manual

https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.sph_harm.html

scipy.special.sph_harm — SciPy v1.14.0 Manual

 

https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fsolve.html

fsolve — SciPy v1.14.0 Manual

 

 

파이썬 시각화중 set_facecolor 부분은 아래를 참고했습니다.

https://github.com/liam-ilan/electron-orbitals/blob/master/generator/render_3d.py

electron-orbitals/generator/render_3d.py at master · liam-ilan/electron-orbitals