有限方势阱

现代电子结构理论作业

  • 作业:有限方势阱
  • 日期:2021年4月14日
1
2
3
4
5
6
import numpy as np
import matplotlib as mpl
from matplotlib import pyplot as plt

mpl.rcParams['font.size'] = 24
mpl.rcParams['figure.figsize'] = (16, 10)

波函数的形式

一维空间中运动的粒子,势垒在有限区域($0 < x < a$)为常量 $U_0 > 0$,区域外为 $0$。即:
\begin{align}
U(x) &= 0,& x < 0\ {\rm or}\ x > a\\
U(x) &= U_0,& 0 < x < a\\
\end{align}
粒子的波函数 $\Psi$ 满足薛定谔方程:
$$- \frac{\hbar^2}{2m} \frac{d^2 \psi}{dx^2} + U\psi = E\psi$$
带入上式,可得:
\begin{align}
\frac{d^2 \psi}{dx^2} + \frac{2m}{\hbar^2} E\psi &= 0,& x < 0\ {\rm or}\ x > a\\
\frac{d^2 \psi}{dx^2} + \frac{2m}{\hbar^2} (E - U_0)\psi &= 0,& 0 < x < a\\
\end{align}

$E > U_0$ 时

令:
\begin{align}
k_1 &= \sqrt{\frac{2m}{\hbar^2} E}\\
k_2 &= \sqrt{\frac{2m}{\hbar^2} (E - U_0)}\\
\end{align}
则有:
\begin{align}
\frac{d^2 \psi}{dx^2} + k_1^2\psi &= 0,& x < 0\ {\rm or}\ x > a\\
\frac{d^2 \psi}{dx^2} + k_2^2\psi &= 0,& 0 < x < a\\
\end{align}
其中 $k_1$ 和 $k_2$ 均为实数且大于零。
方程的解为:
\begin{align}
\psi_1 &= A e^{ik_1x} + A’ e^{-ik_1x},& x < 0\\
\psi_2 &= B e^{ik_2x} + B’ e^{-ik_2x},& 0 < x < a\\
\psi_3 &= C e^{ik_1x} + C’ e^{-ik_1x},& x > a\\
\end{align}
其中第一项是正向传播(入射)的波函数,第二项是反向传播(反射)的波函数。由于在 $x > a$ 区域内不存在反射,所以 $C’ = 0$。

由条件 $\psi_1 \big \lvert_{x = 0} = \psi_2 \big \lvert_{x = 0}$,得:
$$A + A’ = B + B’$$
由条件 $\frac{d\psi_1}{dx} \big \lvert_{x = 0} = \frac{d\psi_2}{dx} \big \lvert_{x = 0}$,得:
$$k_1 A - k_1 A’ = k_2 B - k_2 B’$$
由条件 $\psi_2 \big \lvert_{x = a} = \psi_3 \big \lvert_{x = a}$,得:
$$B e^{ik_2a} + B’ e^{-ik_2a} = C e^{ik_1a}$$
由条件 $\frac{d\psi_2}{dx} \big \lvert_{x = a} = \frac{d\psi_3}{dx} \big \lvert_{x = a}$,得:
$$k_2 B e^{ik_2a} - k_2 B’ e^{-ik_2a} = k_1 C e^{ik_1a}$$
解得:
$$A’ = \frac{2i (k_1^2 - k_2^2) sin(k_2a)}{(k_1 - k_2)^2 e^{ik_2a} - (k_1 + k_2)^2 e^{-ik_2a}} A$$
$$B = \frac{1}{2} \big(\frac{k_1(A - A’)}{k_2} + A + A’\big)$$
$$B’ = A + A’ - B$$
$$C = \frac{4k_1k_2 e^{-ik_1a}}{(k_1 + k_2)^2 e^{-ik_2a} - (k_1 - k_2)^2 e^{ik_2a}} A$$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
def get_coeff_by_E_gt_U_0(init_a, init_A, init_E, init_U_0, para_m, para_hbar):
assert init_E > init_U_0

k_1 = np.sqrt((2 * para_m / para_hbar ** 2) * init_E)
k_2 = np.sqrt((2 * para_m / para_hbar ** 2) * (init_E - init_U_0))

A = init_A
A_p = (2j * (k_1 ** 2 - k_2 **2) * np.sin(k_2 * init_a) / ((k_1 - k_2) ** 2 * np.exp(1j * k_2 * init_a) - (k_1 + k_2) ** 2 * np.exp(-1j * k_2 * init_a))) * A

B = 1 / 2 * (k_1 * (A - A_p) / k_2 + A + A_p)
B_p = A + A_p - B

C = (4 * k_1 * k_2 * np.exp(-1j * k_1 * init_a) / ((k_1 + k_2) ** 2 * np.exp(-1j * k_2 * init_a) - (k_1 - k_2) ** 2 * np.exp(1j * k_2 * init_a))) * A
C_p = 0

return k_1, k_2, A, A_p, B, B_p, C, C_p

$E < U_0$ 时

$k_2$ 是虚数,令:
$$k_2 = ik_3$$
其中 $k_3$ 是实数。满足:
$$k_3 = \sqrt{\frac{2m}{\hbar^2} (U_0 - E)}$$
且有:
$$C = \frac{2ik_1k_3 e^{-ik_1a}}{(k_1^2 - k_3^2) sh(k_3a) + 2ik_1k_3 ch(k_3a)} A$$
等式中 $sh$ 和 $ch$ 分别为双曲正弦和双曲余弦函数:
\begin{align}
sh(x) &= \frac{e^x - e^{-x}}{2}\\
ch(x) &= \frac{e^x + e^{-x}}{2}\\
\end{align}
其他参数 $A’$、$B$ 和 $B’$ 同上。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
def get_coeff_by_E_lt_U_0(init_a, init_A, init_E, init_U_0, para_m, para_hbar):
assert init_E < init_U_0

k_1 = np.sqrt((2 * para_m / para_hbar ** 2) * init_E)
k_3 = np.sqrt((2 * para_m / para_hbar ** 2) * (init_U_0 - init_E))
k_2 = 1j * k_3

A = init_A
A_p = (2j * (k_1 ** 2 - k_2 **2) * np.sin(k_2 * init_a) / ((k_1 - k_2) ** 2 * np.exp(1j * k_2 * init_a) - (k_1 + k_2) ** 2 * np.exp(-1j * k_2 * init_a))) * A

B = 1 / 2 * (k_1 * (A - A_p) / k_2 + A + A_p)
B_p = A + A_p - B

C = ((2j * k_1 * k_3 * np.exp(-1j * k_1 * init_a)) / ((k_1 ** 2 - k_3 ** 2) * (np.exp(k_3 * init_a) - np.exp(k_3 * init_a * -1)) / 2 + 2j * k_1 * k_3 * (np.exp(k_3 * init_a) + np.exp(k_3 * init_a * -1)) / 2)) * A
C_p = 0

return k_1, k_2, A, A_p, B, B_p, C, C_p

波函数隧穿示意图

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
def plot_wave(init_a, init_A, init_E, init_U_0, para_m, para_hbar):
if init_E > init_U_0:
k_1, k_2, A, A_p, B, B_p, C, C_p = get_coeff_by_E_gt_U_0(init_a, init_A, init_E, init_U_0, para_m, para_hbar)

else:
k_1, k_2, A, A_p, B, B_p, C, C_p = get_coeff_by_E_lt_U_0(init_a, init_A, init_E, init_U_0, para_m, para_hbar)

print("Coefficients: ")
print("-" * 20)
print("k1 = {}\nk2 = {}".format(k_1, k_2))
print()
print("-" * 5)
print("A = {}\nA' = {}\n\nB = {}\nB' = {}\n\nC = {}\nC' = {}".format(A, A_p, B, B_p, C, C_p))

plt.fill_between((0, init_a), (0, 0), (init_U_0, init_U_0), color="grey", alpha=0.3)

plt.plot((0, 0), (0, init_U_0), ls="--", color="grey")
plt.plot((init_a, init_a), (0, init_U_0), ls="--", color="grey")
plt.plot((0, init_a), (init_U_0, init_U_0), ls="--", color="grey")

plt.axhline(0, ls="--", color="grey")
plt.axhline(init_E, ls="-.", color="grey")

x_1 = np.linspace(-5, 0, 100)
x_2 = np.linspace(0, init_a, 100)
x_3 = np.linspace(init_a, init_a + 5, 100)

psi_1 = A * np.exp(1j * k_1 * x_1) + A_p * np.exp(-1j * k_1 * x_1)
psi_2 = B * np.exp(1j * k_2 * x_2) + B_p * np.exp(-1j * k_2 * x_2)
psi_3 = C * np.exp(1j * k_1 * x_3) + C_p * np.exp(-1j * k_1 * x_3)

psi_r_1 = (psi_1.real) + init_E
psi_r_2 = (psi_2.real) + init_E
psi_r_3 = (psi_3.real) + init_E

plt.plot(x_1, psi_r_1, color="red", label="real part")
plt.plot(x_2, psi_r_2, color="yellow")
plt.plot(x_3, psi_r_3, color="orange")

psi_i_1 = (psi_1.imag) + init_E
psi_i_2 = (psi_2.imag) + init_E
psi_i_3 = (psi_3.imag) + init_E

plt.plot(x_1, psi_i_1, color="red", ls="--", alpha=0.5, label="image part")
plt.plot(x_2, psi_i_2, color="yellow", ls="--", alpha=0.5)
plt.plot(x_3, psi_i_3, color="orange", ls="--", alpha=0.5)

plt.legend()

plt.xlim(-5, init_a + 5)

plt.xlabel("Position")
plt.ylabel("Potential")

示例 1

1
2
3
4
5
6
7
8
9
init_a = 0.8
init_A = 1.0
init_E = 18.0
init_U_0 = 20.0

para_m = 1.0
para_hbar = 1.0

plot_wave(init_a, init_A, init_E, init_U_0, para_m, para_hbar)
Coefficients: 
--------------------
k1 = 6.0
k2 = 2j

-----
A = 1.0
A' = (0.7520265700863319-0.6119552683751504j)

B = (1.7939461876058913-0.6779377790580774j)
B' = (-0.04191961751955953+0.06598251068292704j)

C = (-0.17568770176257886+0.17058903493914326j)
C' = 0

示例 1

示例 2

1
2
3
4
5
6
7
8
9
init_a = 2.4
init_A = 1.0
init_E = 18.0
init_U_0 = 20.0

para_m = 1.0
para_hbar = 1.0

plot_wave(init_a, init_A, init_E, init_U_0, para_m, para_hbar)
Coefficients: 
--------------------
k1 = 6.0
k2 = 2j

-----
A = 1.0
A' = (0.7999219735365254-0.6000227522136864j)

B = (1.7999951150887923-0.6001284158020551j)
B' = (-7.314155226678665e-05+0.0001056635883687207j)

C = (0.006089318000316135-0.007775171099988885j)
C' = 0

示例 2

示例 3

1
2
3
4
5
6
7
8
9
init_a = 2.4
init_A = 1.0
init_E = 18.0
init_U_0 = 10.0

para_m = 1.0
para_hbar = 1.0

plot_wave(init_a, init_A, init_E, init_U_0, para_m, para_hbar)
Coefficients: 
--------------------
k1 = 6.0
k2 = 4.0

-----
A = 1.0
A' = (0.013645635554689332-0.0711485628631865j)

B = (1.2465885911113275+0.01778714071579663j)
B' = (-0.23294295555663824-0.08893570357898314j)

C = (0.07308505153279372+0.9946910343796513j)
C' = 0

示例 3

示例 4

1
2
3
4
5
6
7
8
9
init_a = 2.4
init_A = 1.0
init_E = 18.0
init_U_0 = 40.0

para_m = 1.0
para_hbar = 1.0

plot_wave(init_a, init_A, init_E, init_U_0, para_m, para_hbar)
Coefficients: 
--------------------
k1 = 6.0
k2 = 6.6332495807108j

-----
A = 1.0
A' = (-0.09999999999999416-0.9949874371065912j)

B = (0.89999999999999-0.9949874371066029j)
B' = (1.5765166949677223e-14+1.1657341758564144e-14j)

C = (-8.615766211993087e-08-2.268307815735802e-07j)
C' = 0

示例 4

1
2


文章作者: 喵函数
文章链接: https://eigenmiao.site/2021/04/14/homework-01/
版权声明: 本博客所有文章除特别声明外,均采用 CC BY-NC-SA 4.0 许可协议。转载请注明来自 本征喵的小站