Phương trình dừng Schrodinger

Vấn đề cơ bản

Trong bài “Sự hình thành trạng thái dừng“, ta đã đi đến kết luận rằng, chỉ khi năng lượng có giá trị cụ thể ở một vài mức nhất định, rời rạc, sóng trong hố thế mới ổn định và đạt đến trạng thái dừng. Khi ấy, tại mỗi điểm trong không gian, sóng chỉ dao động tại chỗ, không di chuyển. Hàm sóng đặc trưng cho trạng thái phải có dạng:

\[\psi_E(x,t)=\Psi(x)e^{-i(E/\hbar)t}.\label{eq:1}\tag{1}\]

Chỉ số \(E\) kí hiệu ở đây ý nói rằng \(\psi_E(x,t)\) là hàm tương ứng với trạng thái dừng, có mức năng lượng \(E\) xác định. Hàm \(\Psi(x)\) chỉ phụ thuộc vào toạ độ, không phụ thuộc vào thời gian. Nó chỉ ra biên độ dao động của hàm sóng tại mỗi điểm trong không gian. Tại mỗi vị trí \(x\), sóng dao động tại chỗ với biên độ \(\Psi(x)\) và tần số \(\omega=E/\hbar\). Bản thân hàm \(\Psi(x)\) được gọi là hàm biên độ. Trong nhiều tài liệu khác, \(\Psi(x)\) cũng được gọi một cách chưa chính xác là hàm sóng, bởi vì \(\psi_E(x,t)=\Psi(x)e^{-i(E/\hbar)t}\) với sự vận động theo thời gian mới thực sự là sóng.

Một trong những bài toán cơ bản của cơ học lượng tử là đi tìm dạng của hàm biên độ \(\Psi(x)\).

Phương trình dừng Schrodinger

Sóng dừng \eqref{eq:1} phải là nghiệm của phương trình Schrodinger:

\[i\hbar\frac{\partial}{\partial t}\psi_E(x,t)=\left[-\frac{\hbar^2}{2m}\frac{\partial^2}{\partial x^2}+U(x)\right]\psi_E(x,t).\label{eq:2}\tag{2}\]

Thế \eqref{eq:1} vào \eqref{eq:2} ta có:

\[i\hbar\frac{\partial}{\partial t}\left(\Psi(x)e^{-i(E/\hbar)t}\right)=\left[-\frac{\hbar^2}{2m}\frac{\partial^2}{\partial x^2}+U(x)\right]\left(\Psi(x)e^{-i(E/\hbar)t}\right).\]

Dấu \(\partial/\partial t\) ở đây biểu thị cho đạo hàm theo thời gian, còn \(\partial^2/\partial x^2\) biểu thị đạo hàm hai lần theo toạ độ. Sau vài phép biến đổi đơn giản, thu được phương trình:

\[-\frac{\hbar^2}{2m}\frac{\partial^2}{\partial x^2}\Psi(x)+U(x)\Psi(x)=E\Psi(x).\label{eq:3}\tag{3}\]

Ta có thể viết lại phương trình \eqref{eq:3} dưới dạng thực hành:

\[\Psi”(x)=-k\left[E-U(x)\right]\Psi(x),\label{eq:4}\tag{4}\]

với \(k=\dfrac{2m}{\hbar^2}\). Phương trình \eqref{eq:3} hoặc \eqref{eq:4} được gọi là phương trình dừng Schrodinger. Đây là một phương trình vi phân bậc hai, theo lý thuyết hoàn toàn có thể giải được một khi đã biết dạng hố thế \(U(x)\) và năng lượng \(E\) của hạt. Tuy vậy không phải mức năng lượng \(E\) nào cũng cho ra kết quả phù hợp với tự nhiên. Một hàm biên độ \(\Psi(x)\) phù hợp với tự nhiên cần thoả điều kiện biên:

\[\begin{aligned}\Psi(x=-\infty)&=0\\ \Psi(x=+\infty)&=0\end{aligned}\]

Mật độ của đám mây hạt cho trạng thái dừng có thể tính theo hàm biên độ:

\[\psi_E(x,t)^*\psi_E(x,t)=\Psi(x)e^{i(E/\hbar)t}\cdot\Psi(x)e^{-i(E/\hbar)t}=[\Psi(x)]^2.\]

Rõ ràng mật độ vi hạt tại trạng thái dừng chỉ phụ thuộc vào toạ độ, không phụ thuộc vào thời gian.

Giải phương trình dừng Schrodinger trên máy tính

Việc giải phương trình dừng Schrodinger \eqref{eq:4} thường rất khó khăn, cần đến rất nhiều sự hỗ trợ toán học từ các hàm siêu việt. Đôi khi việc giải thậm chí không thể tiến hành bằng biến đổi giải tích nếu hố thế có dạng phức tạp. Bởi vậy trong giáo trình này chúng ta tiếp cận lời giải phương trình Schrodinger theo phương pháp số, tức giải trên máy tính.

Ta sẽ giải phương trình Schrodinger được viết dưới dạng \eqref{eq:4}:

\[\Psi”(x)=-k\left[E-U(x)\right]\Psi(x),\]

với \(k=\dfrac{2m}{\hbar^2}\). \(U(x)\) là hàm thế năng có phương trình phụ thuộc vào hình dạng của hố thế. Hình 1 minh hoạ cho trường hợp hố thế vuông, rộng 4 A, sâu 80 eV.

Hình 1: Mức năng lượng \(E\) thấp hơn chiều cao hố thế vuông

Nếu đặt:

\[\begin{cases}y_1=\Psi(x),\\ y_2=\Psi'(x),\end{cases}\]

có thể viết lại phương trình vi phân bậc hai trên thành hệ phương trình vi phân bậc nhất:

\[\begin{cases}\dfrac{dy_1}{dx}=\Psi'(x)=y_2,\\ \dfrac{dy_2}{dx}=\Psi”(x)=-k\left[E-U(x)\right]y_1.\end{cases}\label{eq:5}\tag{5}\]

Hệ phương trình \eqref{eq:5} trình bày theo code Matlab như sau:

1
2
3
4
5
6
function dy = syst(x,y)
%% Schrodinger Equation
global E k U0 x1 x2
dy = zeros(2,1);
dy(1) = y(2);
dy(2) = -k<em>(E-U_square(U0,x1,x2,x))</em>y(1);

Kí hiệu \(\mathrm{dy(1)}\) biểu thị cho đạo hàm bậc nhất \(\Psi'(x)\), còn \(\mathrm{dy(2)}\) biểu thị cho đạo hàm bậc hai \(\Psi”(x)\). Hàm \(\mathrm{U\_square(U0,x1,x2,x)}\) sử dụng ở đây chỉ là một ví dụ dành cho hố thế vuông với định nghĩa:

\[\begin{cases}U=U_0&\mathrm{khi}\ x\in(x_1,x_2),\\ U=0&\mathrm{khi}\ x\notin(x_1,x_2).\end{cases}\]

Trên Matlab hàm thế năng \(U(x)\) trên thực hiện qua function:

1
2
3
4
5
6
7
8
9
10
11
function y = U_square(U0,x1,x2,x)
%% Potential Energy
lenx = length(x);
y = zeros(1,lenx);
for i = 1:lenx
if(x(i)<=x1)||(x(i)>=x2)
y(i) = 0;
else
y(i) = U0;
end
end

Với những trường hợp hố thế khác, như hố parabol, hố tuần hoàn… người viết chương trình có thể định nghĩa hàm thế năng \(U(x)\) theo dạng khác.

Phương trình vi phân giải trên Matlab bằng phương pháp Runge-Kutta bậc 4-5:

1
[x,psi] = ode45(@syst,linspace(xmin_cal,xmax_cal,1000),[psi psi1]);

Trong đó \(\mathrm{xmin\_cal}\) và \(\mathrm{xmax\_cal}\) là hai vị trí nằm hai bên rìa hố thế, đủ xa về hai phía sao cho có thể xem như:

\[\begin{cases}\mathrm{xmin\_cal}=-\infty,\\ \mathrm{xmax\_cal}=+\infty.\end{cases}\]

Còn \(\mathrm{psi}\) và \(\mathrm{psi1}\) giúp khai báo điều kiện ban đầu của \(\Psi\) và \(\Psi’\), với giá trị đủ nhỏ:

\[\begin{cases}\Psi(-\infty)\approx\Psi(\mathrm{xmin\_cal})=10^{-10},\\ \Psi'(-\infty)\approx\Psi'(\mathrm{xmin\_cal})=10^{-10}.\end{cases}\]

Kết quả tính toán thể hiện trên hình 2. Khi năng lượng \(E=-75\,\mathrm{eV}\), việc giải phương trình vi phân \eqref{eq:4} cho ra nghiệm có đồ thị màu đỏ. Đường màu xanh tương ứng với nghiệm khi \(E=-72.368\,\mathrm{eV}\). Cả hai đều thoả điều kiện biên \(\Psi(-\infty)\approx 0\), bởi vì ta cố tình đặt như thế. Nhưng điều kiện biên \(\Psi(+\infty)\approx 0\) chỉ đường màu xanh thoả mãn. Như vậy trạng thái ứng với năng lượng \(E=-72.368\,\mathrm{eV}\) là một trạng thái dừng, tương ứng với một nghiệm đích thực của phương trình Schrodinger.

giải phương trình Schrodinger
Hình 2: So sánh nghiệm hội tụ và nghiệm phân kì

Có thể tổng kết lại quy trình giải phương trình dừng Schrodinger như sau:

  • Bước 1: Gán giá trị năng lượng \(E\) sát gần đáy hố thế, cao hơn đáy chỉ một chút.
  • Bước 2: Chọn vị trí \(\mathrm{xmin\_cal}\) đủ xa về phía bên trái của hố thế, sao cho có thể xem như \(\mathrm{xmin\_cal}\approx -\infty\). Cũng như chọn vị trí \(\mathrm{xmax\_cal}\) đủ xa về phía bên phải của hố thế, sao cho có thể xem như \(\mathrm{xmax\_cal}\approx +\infty\). Gán giá trị ban đầu cho \(\Psi(\mathrm{xmin\_cal})=10^{-10}\approx 0\) và \(\Psi'(\mathrm{xmin\_cal})=10^{-10}\approx 0\).
  • Bước 3: Giải phương trình vi phân (4) với giá trị năng lượng \(E\) đã cho.
  • Bước 4: Đánh giá nghiệm \(\Psi(x)\) vừa thu được. Nếu \(\Psi(x)\) bị phân kì mạnh ở rìa bên phải hố thế, có thể suy rằng nó không phải là nghiệm của phương trình Schrodinger, do không thoả mãn điều kiện biên \(\Psi(+\infty)\approx 0\). Ngược lại, nếu \(\Psi(x)\) hội tụ, áp sát vào trục hoành, chứng tỏ nó là nghiệm.

Nếu \(\Psi(x)\) phân kì, chứng tỏ năng lượng \(E\) không phải giá trị phù hợp, ta có thể nâng \(E\) lên một chút và lặp lại quy trình giải, cho đến khi nào tìm thấy nghiệm hội tụ.

Mẫu chương trình hoàn chỉnh

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
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
function Schrodinger_solve_model

clc;
clear variables
close all

global E k U0 x1 x2

%% CONSTANTS
h = 1.054e-34;
m = 9.1095e-31;
q_element = 1.6e-19; % eV -> J
Angstrom = 1e-10; % m
k = 2*m/h/h;

%% INPUT PARAMETRS
E = -72.37; % eV

U0 = -80*q_element;
x1 = -2*Angstrom;
x2 = 2*Angstrom;

Emin = -100*q_element;
Emax = 100*q_element;

xmin = -6*Angstrom;
xmax = 6*Angstrom;
xmin_cal = -6*Angstrom;
xmax_cal = 3*Angstrom;

%% DATA PROCESSING
E = E*q_element;
U = U_square(U0,x1,x2,x1);

psi = 1e-10;
psi1 = 1e-10;

%% CALCULATION
[x,psi] = ode45(@syst,linspace(xmin_cal,xmax_cal,1000),[psi psi1]);
psi_max = max(psi(:,1));
psi_draw = psi(:,1)/psi_max*Emax/q_element;
psi_psi = (psi(:,1)/psi_max*Emax/q_element).^2*2e-2;

%% FIGURE
x_plot = linspace(xmin,xmax,1000);
U_plot = U_square(U0,x1,x2,x_plot);

figure('name','Psi function','color','black','numbertitle','off');
hold on
plot(x_plot/Angstrom,U_plot/q_element,'linewidth',2,'color','w');
line([xmin/Angstrom,xmax/Angstrom],[E/q_element,E/q_element],'linewidth',2,'color','b');
plot(x/Angstrom,psi_draw,'linewidth',1,'color','y');

axis([xmin/Angstrom xmax/Angstrom Emin/q_element Emax/q_element]);
set(gca,'color','k','xcolor','w','ycolor','w')
xlabel('x [Angstrom]');
ylabel('E, U [eV]');

function dy = syst(x,y)
%% Schrodinger Equation
global E k U0 x1 x2
dy = zeros(2,1);
dy(1) = y(2);
dy(2) = -k*(E-U_square(U0,x1,x2,x))*y(1);

function y = U_square(U0,x1,x2,x)
%% Potential Energy
lenx = length(x);
y = zeros(1,lenx);
for i = 1:lenx
if(x(i)<=x1)||(x(i)>=x2)
y(i) = 0;
else
y(i) = U0;
end
end