[찐따 에세이] 파이썬을 활용한 양자 현상의 수치적 연구 - 슈뢰딩거 방정식을 중심으로

 

 

 

파이썬을 활용한 양자 현상의 수치적 연구 - 슈뢰딩거 방정식을 중심으로

 

들어가기 전에

 

<< 본 문서는 고졸 찐따가 작성하였습니다. >>

 

양자 현상의 수치적 연구를 하기 위한 목적으로, The Crank Nicholson 기법을 통해 슈뢰딩거 방정식 증명하는 과정을 바탕으로 Thomas 알고리즘을 활용하여 파이썬으로 시뮬레이터 프로토타입 개발을 목적으로 작성된 에세이입니다.

 

개요

$$
i\hbar\frac{\partial}{\partial t} \Psi(x, t) = \hat{H}\Psi(x, t)
$$

슈뢰딩거 방정식은 양자역학에서 파동 함수의 시간 변화를 좌우한다. 다양한 전위 분포에 대한 슈뢰딩거 방정식의 해를 풀고 연구하는 것은 입문 양자역학 학습과 연구에 중요한 부분이다. 그러나 슈뢰딩거 방정식의 해는 해석학적 폐쇄형 해(analytic closed-form solutions)를 명확하게 제시할 수 있는 경우가 몇 가지 밖에 되지 않는다.

이렇듯 손에 꼽을 정도의 소수의 경우에서도 해석학적 해를 찾기 위해서는 슈뢰딩거 방정식의 연구가 불투명한 노력이 되어버리는 일종의 '수학적 트릭'을 사용해야만 한다. 이 '수학적 트릭'은 슈뢰딩거 방정식에 대한 연구를 무의미한 노력으로 만들어버리는 문제가 발생한다. 따라서 임의의 초기 조건 및 전위 분포에 대해서 슈뢰딩거 방정식을 풀기 위한 수치해석은 다양한 양자역학 시나리오를 시뮬레이션하고 연구하는데 도움을 줄 수 있기 때문에 가치가 있다고 판단된다.

이 문서에서는 1D에서 슈뢰딩거 방정식을 해결하기 위해 수치해석 기법으로 파이썬 코드를 사용하여 양자 현상에 대한 응용에 대해 논의하고자 한다.

수치해석 (Numerical Method)

1D에서의  TDSE(Time Dependent Schrodinger Equation)는 편미분 방정식이다. 따라서, TDSE를 수치적분을 하기 위해서는 분명 유한차분법이 필요하다. Robertson은 TDSE를 수치해석 하기 위해서 The Crank-Nicholson 기법을 사용한 이유에 대해 설명한다.

이 문서에서는 The Crank-Nicholson 기법에 대한 설명과 해당 기법을 통해 파이썬 코드를 작성하여 알고리즘을 구현하고자 한다.

The Crank Nicholson Technique

Crank Nicholson 기법은 ImEx(Implicit-Explicit) 기법 중 하나이다. Crank Nicholson 기법은 완전행렬(혹은 완성행렬; complete matrix) 방정식을 풀지 않고도 절대 안정적인(unconditionally stable) 장점이 있다. 이는 적어도 어떤 방면으로는 명시적 방법과 암묵적 방법 모두 최고의 특징을 제공한다. (Thomas W. Baumgarte and Stuart L. Shapiro.Numerical Relativity: Solving Einstein’s Equa-tions on the Computer. 2010.)

아래의 예시에서 보게 되겠지만, 각각의 타임 인스턴스에서 프로그램은 띠행렬(tri-diagonal matrix)만 풀면 되며, 이것은 임의의 행렬 방정식을 푸는 것보다 훨씬 더 효율적이다. 이것이 어떤 방식으로 설계되었는지 간략히 증명하자면 다음과 같다.

이 증명을 위해 먼저, 슈뢰딩거 방정식부터 시작해야 한다.

$$
i\hbar\frac{\partial}{\partial t} \Psi(x, t) = \hat{H}\Psi(x, t)

\quad\quad\quad\quad (1)
$$

$\Delta{x}$ 에 대한 테일러 급수 전개(Taylor series expansion)를 사용하여 다음의 방정식을 도출해낼 수 있다.

$$
\Psi(x + \Delta{x}) = \Psi(x) + \Delta{x}\frac{\partial \Psi(x)}{\partial x} + \frac{1}{2} \Delta{x^2} \frac{\partial^2 \Psi(x)}{\partial x^2} + O(\Delta{x}^3)

\quad\quad\quad\quad (2)
$$

이와 비슷하게, $x$ 에 대한 테일러 급수를 $-\Delta{x}$ 를 사용하여 확장하면 다음의 방정식을 얻을 수 있다.

$$
\Psi(x - \Delta{x}) = \Psi(x) - \Delta{x} \frac{\partial x}{\partial \Psi(x)} + \frac{1}{2}\Delta{x^2} \frac{\partial x^2}{\partial^2 \Psi(x)} - O(\Delta{x^3})

\quad\quad\quad\quad (3)
$$

이렇게 확장된 두 방정식을 결합하면 다음과 같은 식을 도출해낼 수 있다.

$$
\Psi(x + \Delta{x}) + \Psi(x - \Delta{x}) = + 2\Psi(x) + \Delta{x^2} \frac{\partial x^2}{\partial^2 \Psi(x)} + O(\Delta{x^4})

\quad\quad\quad\quad (4)
$$

고차항을 소거한 후 재배열함으로써 2차 도함수를 유한으로 분리할 수 있다.

$$
\frac{\partial x^2}{\partial^2 \Psi(x)} \approx \frac{\Psi(x + \Delta{x}) + \Psi(x - \Delta{x}) - 2\Psi(x)}{\Delta{x^2}}

\quad\quad\quad\quad (5)
$$

이제 $n \Delta{t}$ 와 $i \Delta{x}$ 점의 이산 메쉬(discrete mesh)에 대한 그리드 표기법(grid notation)으로 전환한다.

$\Psi_{i}^{n}$ 은 공간적으로 인덱스 $i$ 로 인덱싱되고, 시간적으로는 인덱스 $n$ 으로 인덱싱된다. 이는 $\Psi(t, x \pm \Delta{x}) = \Psi_{i \pm 1}^{n}$ 혹은 $\Psi(t \pm \Delta{t}, x) = \Psi_{i}^{n \pm 1}$ 을 의미하며, 참고로 이산 방정식을 공식화하는 방법을 선택하는 데 있어 자유로울 수 있다는 점에 유의해야 한다. 유일한 제약 조건은 극한값 $\Delta{x} = \Delta{t} \rightarrow 0$ 이므로 공식을 슈뢰딩거 방정식으로 줄여야 한다는 것이다. 따라서 두 시간 단계에서 파동 함수의 평균을 사용하여 슈뢰딩거 방정식을 만족하도록 할 수 있다.

$$
i\hbar(\frac{\Psi^{n+1} - \Psi^n)}{\Delta{t}}) = \frac{1}{2}(\hat{H}\Psi^{n+1} + \hat{H}\Psi^n)

\quad\quad\quad\quad (6)
$$

다음을 얻기 위해 $\Psi^{n+1}$ 항을 분리하도록 재배열할 수 있다.

$$
\Psi^{n+1} = (\frac{1}{1 + \frac{i\Delta{t}}{2\hbar}\hat{H}})(1 - \frac{i\Delta{t}}{2\hbar}\hat{H})\Psi^n

\quad\quad\quad\quad (7)
$$

해밀턴 연산자 $\hat{H}$ 는 다음과 같이 정의된다.

$$
\hat{H} = \frac{-h^2}{2m} \frac{\partial^2}{\partial x^2} + V(x)

\quad\quad\quad\quad (8)
$$

해밀턴 연산자에서 $(5)$ 의 방정식의 유한차분 근사치를 대입한 다음, $(7)$ 의 방정식에 다시 대입하여 시공간 모두 유한 근사치를 얻을 수 있다. 치환과 약간의 대수 후에 다음과 같은 방정식을 도출해낼 수 있다.

$$
\Psi_{i+1}^{n+1} + \Psi_{i-1}^{n+1} + A_{i} \Psi_{i}^{n+1} = B_{i}

\quad\quad\quad\quad (9)
$$

여기서 $A_{i}$ 는 다음과 같다.

$$
A_{i} = -2 + \frac{4 i m \Delta{x^2}}{\hbar \Delta{t}} - \frac{2 m \Delta{x^2}}{h^2} V_{i}

\quad\quad\quad\quad (10)
$$

그리고 $B_{i}$ 는 다음과 같다.

$$
B_{i} = -\Psi_{i+1}^{n} - \Psi_{i-1}^{n} + \Psi_{i}^{n}(2 + \frac{4 i m \Delta{x^2}}{\hbar \Delta{t}} + \frac{2 m \Delta{x^2}}{\hbar^2} V_{i})

\quad\quad\quad\quad (11)
$$

$(9)$ 의 방정식은 LHS에서 오로지 $\Psi^{n+1}$ 만 가지며, RHS에서는 오로지 $\Psi{n}$ 만 갖는다. 하지만 Crank Nicholson 기법은 명시적인 방법이 아니며, $(9)$ 의 방정식을 풀려면 family equations을 풀어야 한다. 만약 초기조건을 $\Psi_{i}^{0}$ , 경계조건을 $\Psi_{0}^{n}$ 와 $\Psi_{L}^{n}$ 로 알고 있다면, $(9)$ 방정식을 통해 다음과 같은 띠행렬을 도출해낼 수 있다.

$$
\begin{pmatrix}
 &A_{1}  &1  &0  &0 &... &  \\
 &1  &A_{2}  &1  &0 &... &  \\
 &0  &1  &A_{2}  &1 &... &  \\
 &..  &..  &.. &.. &... &  \\
 &0  &0  &..  &1 &A_{L-1} &  \\
\end{pmatrix}

\begin{pmatrix}
\Psi_{1}^{n+1}    \\
\Psi_{2}^{n+1}    \\
\Psi_{3}^{n+1}    \\
...               \\
\Psi_{L-1}^{n+1}  \\
\end{pmatrix}

=

\begin{pmatrix}
B_{1}    \\
B_{2}    \\
B_{3}    \\
...      \\
B_{L-1}  \\
\end{pmatrix}

\quad\quad\quad\quad (12)
$$

위와 같은 종류의 띠행렬은 전진 스위프(Forward sweep)와 역 스위프(Backward sweep)를 포함하는 Thomas 알고리즘을 통해 매우 효율적으로 풀어낼 수 있다. (Koonin) 띠행렬을 풀 때 편리한 점은 애초에 행렬을 다룰 필요가 없다는 사실에 있다. 이는 벡터를 통해 구할 수 있다. 이와 같은 경우, 비대각선은 단순히 1s이다. 따라서, 주 대각선 벡터와 RHS 상의 벡터만으로 미지수 $\Psi^{n+1}$ 를 풀면 된다. $\Psi$ 에 대한 방정식 집합을 풀기 위해서 알고리즘의 간략한 스냅샷을 제공하고자 한다.

먼저, 두 보조 배열 $R$ 과 $U$ 를 생성해야 한다.

$$
U_{1} = \frac{1}{A_{1}} \, \text{and} \, R_{1} = B_{1}U_{1}

\quad\quad\quad\quad (13)
$$

여기서 $i \in \left [ 2, 3, ..., L-1 \right ]$ 에 대한 전진 스위프를 시작한다.

$$
U_{i} = \frac{1}{A_{i} - U_{i-1}} \, \text{and} \, R_{i} = (B_{i} - R_{i-1}) U_{i}

\quad\quad\quad\quad (14)
$$

이제 $\Psi^{n+1}$ 을 풀기 위해 $R$ 과 $U$ 를 이용해서 역 스위프를 한다.

$$
\Psi_{i}^{n+1} = 
\begin{cases}
R_{L-1}                        &i = L - 1 \\
R_{i} - U_{i} \Psi_{i+1}^{n+1} &i = L < 1 \\
\end{cases}

\quad\quad\quad\quad (15)
$$

여기서 $i$ 는 $L-1, L-2, L-3...1$ 부터 시작한다. 이런 방식으로 $\Psi^{n+1}$에 대해 풀었을 것이다. 이제 간단히 이 과정을 반복하여 $\Psi^{n+1}$을 사용하여 $\Psi^{n+2}$를 얻을 수 있다. 참고로 아래의 시뮬레이션에서는 복소수를 사용하고 있다.

 

# 시간변수에 대한 인덱스 n을 반복한다.
for n in range(0, len(t)-2):
    # 고정된 경계조건을 사용한다.
    # 각각의 변은 상자 속 입자, 즉 무한 포텐셜 우물(infinite potential walls)로 해석될 여지가 있다.
    # 다시 말해, 입자가 무한히 깊은 포텐셜 우물에 갇혀 있어서 나가지 못하는 시스템을 의미한다.
    Psi[n, 0], Psi[n,  -1] = 0, 0

    # 공간변수에 대한 인덱스 i를 반복한다.
    for i in range(1, len(x)-2):
        A[n, i] = -2 + (4j*m*dx**2)/(hbar*dt)\
                     - (2*m*dx**2)/(hbar **2)*V(x[i])
        B[n, i] = -Psi[n, i+1] - Psi[n, i-1]\
                    + Psi[n, i] * (2 + ((4j*m*dx**2)/(hbar*dt))\
                    + ((2*m*dx**2)/(hbar **2))*V(x[i]))

    # 보조 행렬 R과 U 정의한다.
    U[n, 1] = 1 / A[n, 1]
    R[n, 1] = B[n, 1] * U[n, 1]

    # 전진 스위프 (Foward Sweep)
    for i in range(1, len(x) -2):
        U[n, i] = 1/(A[n, i] - U[n, i-1])
        R[n, i] = (B[n, i] - R[n, i-1])*U[n, i]

    N = len(x) - 1
    i = N - 1

    Psi[n+1, N] = R[n, N]

    # 역 스위프 (Backward Sweep)
    while i >= 1:
        Psi[n+1, i] = R[n, i] - U[n, i]*Psi[n+1, i+1]
        i -= 1

 

참고

• Miguel Alcubierre.Introduction to 3+1 numerical relativity. International series of monographson physics. Oxford: Oxford Univ. Press, 2008.
• David G. Robertson. “Solving the Time-Dependent Schrodinger Equation.” In: (Oct. 2011).
• Thomas W. Baumgarte and Stuart L. Shapiro.Numerical Relativity: Solving Einstein’s Equa-tions on the Computer. 2010.
• Tiwari, Kartik. Quantum Mechanics on Python. (2021).