Nguyên tử hidro trường hợp tổng quát

Bài toán hidro tổng quát

Trong mục “Nguyên tử hidro trường hợp đối xứng cầu” ta đã phân tích các trạng thái dừng của nguyên tử hidro mà không xét đến sự quay của đám mây electron. Nói cách khác, ta đã khảo sát nghiêm túc nguyên tử hidro, nhưng chỉ với trường hợp moment quay bằng không. Giờ đây vấn đề nguyên tử hidro cần nhìn nhận lại một cách tổng quát hơn, khi tìm các trạng thái dừng có mức năng lượng xác định của electron trong nguyên tử hidro có tính đến cả sự quay. Các trạng thái dừng này tương ứng với những sóng dừng \(\Psi(x,y,z)\), thoả mãn phương trình Schrodinger:

\[-\frac{\hbar^2}{2m}\left(\frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial y^2}+\frac{\partial^2}{\partial z^2}\right)\psi(x,y,z,t)+U(r)\psi(x,y,z,t)=E\psi(x,y,z,t),\label{eq:1}\tag{1}\]

và tất nhiên, sẽ thuận tiện hơn nếu biểu diễn đạo hàm trong hệ toạ độ cầu:

\[\frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial y^2}+\frac{\partial^2}{\partial z^2}=\frac{1}{r^2}\frac{\partial}{\partial r}\left(r^2\frac{\partial}{\partial r}\right)+\frac{1}{r^2\sin\vartheta}\frac{\partial}{\partial\vartheta}\left(\sin\vartheta\frac{\partial}{\partial\vartheta}\right)+\frac{1}{r^2\sin^2\vartheta}\frac{\partial^2}{\partial\varphi^2}.\]

Từ đây phương trình Schrodinger \eqref{eq:1} viết thành:

\[-\frac{\hbar^2}{2m}\frac{1}{r^2}\frac{\partial}{\partial r}\left(r^2\frac{\partial}{\partial r}\right)\psi(r,\vartheta,\varphi,t)\\ +\frac{\hbar^2}{2m_e}\frac{1}{r^2}\left[\frac{1}{\sin\vartheta}\frac{\partial}{\partial\vartheta}\left(\sin\vartheta\frac{\partial}{\partial\vartheta}\right)+\frac{1}{\sin^2\vartheta}\frac{\partial^2}{\partial\varphi^2}\right]\psi(r,\vartheta,\varphi,t)\\ +U(r)\psi(r,\vartheta,\varphi,t)\\ =E\psi(r,\vartheta,\varphi,t).\label{eq:2}\tag{2}\]

Ở đây \(r\) – khoảng cách đến tâm nguyên tử, \(\vartheta\) – góc cực, \(\varphi\) – góc phương vị. Phương trình \eqref{eq:2} nói trên chứa trọn biểu thức của toán tử bình phương moment động lượng:

\[\hat{M^2}=-\hbar^2\left[\dfrac{1}{\sin\vartheta}\dfrac{\partial}{\partial\vartheta}\left(\sin\vartheta\dfrac{\partial}{\partial\vartheta}\right)+\dfrac{1}{\sin^2\vartheta}\dfrac{\partial^2}{\partial\varphi^2}\right].\]

Nên có thể viết gọn:

\[-\frac{\hbar^2}{2m}\frac{1}{r^2}\frac{\partial}{\partial r}\left(r^2\frac{\partial}{\partial r}\right)\psi(r,\vartheta,\varphi,t)+\frac{\hat{M^2}}{2m_er^2}\psi(r,\vartheta,\varphi,t)+U(r)\Psi(r,\vartheta,\varphi)=E\psi(r,\vartheta,\varphi,t).\label{eq:3}\tag{3}\]

Rõ ràng hàm sóng \(\psi(r,\vartheta,\varphi,t)\) phải là hàm riêng của phương trình:

\[\hat{M^2}\psi(r,\vartheta,\varphi,t)=M^2\psi(r,\vartheta,\varphi,t).\]

Trong bài “Moment động lượng”, ta đã đi đến khẳng định rằng, hạt có năng lượng \(E\), moment động lượng toàn phần \(M\) cũng như moment động lượng \(M_z\) hoàn toàn xác định sẽ có hàm sóng dạng:

\[\psi(r,\vartheta,\varphi,t)=R(r)\Theta_{lm}(\vartheta)e^{im\varphi}\cdot e^{-i\frac{E}{\hbar}t}.\label{eq:4}\tag{4}\]

Dạng của hàm \(\Theta_{lm}(\vartheta)\) đã khảo sát chi tiết trong bài “Moment động lượng”. Moment quay có giá trị bằng \(M=\hbar\sqrt{l(l+1)}\) với \(l=0,1,2,\ldots\) Còn \(m\) là một số nguyên có độ lớn không vượt quá \(l\). \(R(r)\) là hàm bán kính, phụ thuộc vào hình dạng của hố thế. Thế hàm sóng \eqref{eq:4} và giá trị của moment quay vào \eqref{eq:3}:

\[-\frac{\hbar^2}{2m}\frac{1}{r^2}\frac{\partial}{\partial r}\left(r^2\frac{\partial}{\partial r}\right)\left[R(r)\Theta_{lm}(\vartheta)e^{im\varphi}\cdot e^{-i\frac{E}{\hbar}t}\right]\\ +\frac{\hbar^2l(l+1)}{2m_er^2}\left[R(r)\Theta_{lm}(\vartheta)e^{im\varphi}\cdot e^{-i\frac{E}{\hbar}t}\right]\\ +U(r)\left[R(r)\Theta_{lm}(\vartheta)e^{im\varphi}\cdot e^{-i\frac{E}{\hbar}t}\right]\\ =E\left[R(r)\Theta_{lm}(\vartheta)e^{im\varphi}\cdot e^{-i\frac{E}{\hbar}t}\right].\]

Đạo hàm trong phương trình trên chỉ tác động lên biến \(r\) nên dễ dàng đơn giản hoá, thu được phương trình chỉ chứa hàm bán kính \(R(r)\):

\[-\frac{\hbar^2}{2m}\frac{1}{r^2}\frac{\partial}{\partial r}\left(r^2\frac{\partial}{\partial r}\right)R_{nl}(r)+\frac{\hbar^2l(l+1)}{2m_er^2}R_{nl}(r)+U(r)R_{nl}(r)=ER_{nl}(r).\label{eq:5}\tag{5}\]

Từ đây hàm bán kính \(R_{nl}(r)\) kí hiệu thêm chỉ số \(nl\), bởi nó chỉ phụ thuộc vào hai tham số này, với \(n\) – số thứ tự của hàm \(R(r)\) tìm thấy tính từ mức năng lượng thấp nhất.

Tìm hàm bán kính \(R(r)\) bằng phương pháp số

Thế năng của trường Coulomb do hạt nhân proton tạo ra tiến đến sâu vô cùng tại \(r=0\). Biểu hiện của hàm sóng \(\psi(r)\) tại vị trí đặc biệt này còn chưa rõ. Vì vậy ta sẽ đi tìm \(R_{nl}(r)\) dưới dạng:

\[R_{nl}(r)=\frac{u(r)}{r}.\]

Qua hệ thức này, ta biết chắc điều kiện biên của hàm \(u(r)\) tại \(r=0\) phải bằng \(0\). Bởi nếu \(u(0)>0\), \(\psi(r=0)\) sẽ tiến đến \(\infty\). Nhờ biến đổi toán học:

\[\frac{1}{r^2}\frac{\partial}{\partial r}\left(r^2\frac{\partial}{\partial r}\right)\frac{u(r)}{r}=\frac{1}{r}\frac{\partial^2u(r)}{\partial r^2},\]

phương trình \eqref{eq:5} có thể viết lại theo hàm \(u(r)\):

\[-\frac{\hbar^2}{2m}\frac{\partial^2u(r)}{\partial r^2}+\frac{\hbar^2l(l+1)}{2m_er^2}u(r)+U(r)u(r)=Eu(r).\label{eq:6}\tag{6}\]

Để giải phương trình \eqref{eq:6}, cần viết lại theo thể thực hành:

\[u”(r)=-\frac{2m}{\hbar^2}\left[E-\frac{\hbar^2l(l+1)}{2m_er^2}-U(r)\right]u(r).\label{eq:7}\tag{7}\]

Phương trình \eqref{eq:6} trình bày trong Matlab theo code:

1
2
3
4
5
6
function dy = Schrodinger(r,y)
%% Schrodinger Equation
global h E k m_e l
dy = zeros(2,1);
dy(1) = y(2);
dy(2) = -k<em>(E-h</em>h<em>l</em>(l+1)/2/m_e./r./r-U_hydro(r))*y(1);

Với hàm thế năng hyperbol:

1
2
3
4
5
6
7
8
function y = U_hydro(r)
%% Potential Energy
global epsilon0 q_element
lenr = length(r);
y = zeros(1,lenr);
for i = 1:lenr
y(i) = -1/4/pi/epsilon0<em>q_element</em>q_element/r(i);
end

Phạm vi giải phương trình của chúng ta trải dài từ tâm nguyên tử \(r=0\) đến “rìa” nguyên tử, một khoảng cách lựa chọn khéo léo để có thể xem là đủ xa vô cực. Nghiệm của phương trình \eqref{eq:6} hoặc \eqref{eq:7} phải tuân theo điều kiện biên:

\[\begin{cases}u(0)=0,\\ u(+\infty)=0.\end{cases}\]

Áp dụng phương pháp “bắn tên”, ta sẽ chọn một mức năng lượng \(E\) nào đó rồi bắn thử từ biên bên này và nghiệm xem mũi tên liệu có trúng đích, tức điều kiện biên của phía biên còn lại hay không. Để thuận tiện, ta sẽ bắn mũi tên không phải từ \(r=0\), mà từ ngoài rìa bắn vào với điều kiện ban đầu:

\[u(+\infty)\approx 0,\qquad u'(+\infty)\approx 0.\]

Code Matlab tương ứng có thể viết:

1
2
3
u0 = 1e-10;
u10 = 0;
[r,U] = ode45(@Schrodinger,linspace(rmax_cal,rmin_cal,1000),[u0 u10]);

Ở đây \(\mathrm{rmax\_cal}\) tương ứng với một điểm nằm cách đủ xa hạt nhân nguyên tử. \(\mathrm{rmin\_cal}\) ta chọn sát hạt nhân, nhưng không được trùng với vị trí \(r=0\) của hạt nhân, bởi lực hút sẽ tiến về vô cùng, tâm hố thế sâu vô cực. Trên thực tế, hạt nhân cũng có kích thước, tuy rất nhỏ có thể xem gần như chất điểm, nhưng cũng vừa đủ để tránh lực điện từ vô hạn. Bên trong phạm vi rất nhỏ của hạt nhân có hố thế đặc biệt khác và diễn biến khác, ta không xét ở đây. Có thể chọn \(\mathrm{rmin\_cal}\) tầm một phần nghìn Angstrom. Nhớ rằng bán kính Bohr của hidro bằng \(0.529\,\mathrm{A}\).

Mục đích của chúng ta là phải bắn sao cho hàm \(u(r)\) rơi trúng gốc toạ độ: \(u(0)\approx 0\). Kết quả tính toán cho thấy, không phải mọi giá trị của năng lượng \(E\) đều thoả điều kiện biên này, mà chỉ tại những mức \(E_n\) rời rạc nhất định. Hoá ra những mức rời rạc này đúng theo công thức theo thuyết nguyên tử Bohr đã đề cập trước đó:

\[E_n=-\frac{m_ee^4}{32\pi^2\epsilon_0^2\hbar^2}\cdot\frac{1}{n^2},\qquad n=1,2,3,\ldots\label{eq:8}\tag{8}\]

\(n\) là một số nguyên lớn hơn \(0\) và đặc biệt: phải lớn hơn \(l\). Công thức \eqref{eq:8} không hề chứa các tham số \(l\) và \(m\), có nghĩa năng lượng của electron không phụ thuộc vào sự quay! Đây là điều hết sức đặc biệt, chỉ xảy ra ở nguyên tử hidro. Chúng ta sẽ thấy tình hình khác đi ở những loại nguyên tử khác. Công thức \eqref{eq:8} cũng khẳng định tính đúng đắn của thực nghiệm quang phổ, cũng như trùng khớp với các mức năng lượng theo mẫu nguyên tử Bohr, rằng mức năng lượng tỉ lệ nghịch với bình phương của \(n\).

Những hình dưới đây thể hiện kết quả tính toán hàm bán kính \(R_{nl}(r)\) cho \(6\) trường hợp đầu tiên. Đường màu vàng biểu diễn hàm \(u(r)\), luôn đi qua gốc toạ độ và hội tụ về \(0\) ở rìa nguyên tử, đảm bảo điều kiện biên. Vùng tô màu cam thể hiện mật độ của đám mây electron theo khoảng cách tính từ tâm nguyên tử.

Trường hợp \(n=1,l=0\), electron quanh hạt nhân nguyên tử tạo ra \(1\) bụng sóng.

Hình 1: Trường hợp \(n=1, l=0\)

Trường hợp \(n=2,l=0\), electron tạo ra \(2\) bụng sóng.

Hình 2: Trường hợp \(n=2,l=0\)

Trường hợp \(n=2,l=1\), electron tạo ra \(1\) bụng sóng.

Hình 3: Trường hợp \(n=2,l=1\)

Trường hợp \(n=3,l=0\), electron tạo ra \(3\) bụng sóng.

Hình 4: Trường hợp \(n=3,l=0\)

Trường hợp \(n=3,l=1\), electron tạo ra \(2\) bụng sóng.

Hình 5: Trường hợp \(n=3,l=1\)

Trường hợp \(n=3,l=2\), electron tạo ra \(1\) bụng sóng.

Hình 6: Trường hợp \(n=3,l=2\)

Một cách tổng quát, số bụng sóng do electron tạo ra quanh hạt nhân nguyên tử bằng đúng hiệu số \(n-l\).

“Bán kính” nguyên tử

Từ những kết quả tính toán trên cho hàm bán kính \(R_{nl}(r)\), ta nhận thấy đám mây electron loang dần ra khỏi tâm khi mức năng lượng dâng cao. Như đã một lần nói trong bài “Nguyên tử hidro trường hợp đối xứng cầu”, một cách hình tượng, có thể xem electron như nước rót vào trong cốc, một chiếc cốc hyperbol tròn xoay. Nước càng dâng cao, năng lượng càng đầy, mặt nước càng rộng, bán kính nguyên tử càng lớn. Lý thuyết nguyên tử của Niels Bohr đi đến kết luận thô sơ rằng: bán kính quỹ đạo electron trong nguyên tử hidro tăng dần theo quy luật:

\[r_n=n^2r_{Bohr},\qquad r_{Bohr}=0.529\,\mathrm{A},\]

có nghĩa electron giống như một chất điểm di chuyển trên các quỹ đạo tròn có bán kính \(1\,r_{Bohr}\), \(4\,r_{Bohr}\), \(9\,r_{Bohr}\ldots\). Nhưng cơ học lượng tử nói khác: electron tản ra thành “đám sương mù”, khái niệm quỹ đạo không còn phù hợp nữa. Tuy nhiên vẫn có những nét tương đồng hết sức đặc biệt. Hãy điểm lại các trạng thái, khi đám mây electron chỉ hội tụ tại một đỉnh, khi \(n-l=1\).

Với trường hợp \(n=1,l=0\) trên hình 1, đám mây electron có xu hướng hội tụ đậm đặc nhất quanh vị trí cách tâm nguyên tử đúng bằng \(1\) bán kính \(r_{Bohr}\). Trường hợp \(n=2,l=1\) trên hình 3, đám mây electron có xu hướng hội tụ đậm đặc nhất quanh vị trí cách tâm nguyên tử đúng bằng \(4\) lần bán kính \(r_{Bohr}\). Trường hợp \(n=3,l=2\) trên hình 6, đám mây electron có xu hướng hội tụ đậm đặc nhất quanh vị trí cách tâm nguyên tử đúng bằng \(9\) lần bán kính \(r_{Bohr}\).

Những kết quả trên rất đương đồng với lý thuyết Bohr! Tuy nhiên đó cũng chỉ là sự trùng hợp ngẫu nhiên. Khi electron thay đổi quỹ đạo từ mức năng lượng cao xuống mức năng lượng thấp, nó sẽ phát ra photon để làm nên quang phổ vạch trong máy phân tích quang phổ. Những nghiên cứu về sau cho thấy, bản thân hạt photon cũng mang một moment quay riêng, còn gọi là spin, có giá trị bằng đúng \(1\hbar\). Theo định luật bảo toàn moment động lượng, electron chỉ có thể di chuyển xuống quỹ đạo nào đó, sao cho moment quay giảm đi lượng \(1\hbar\). Điều đó dẫn đến nguyên tắc lọc lựa:

Sự bức xạ và hấp thụ photon chỉ xảy ra giữa hai trạng thái có chỉ số \(l\) chênh nhau đúng \(1\) đơn vị.

Cần nhớ rằng số \(l\) đặc trưng cho moment quay của electron quanh nguyên tử.

Kết luận

Như vậy, vận động theo không gian và thời gian của electron trong nguyên tử hidro được miêu tả qua hàm sóng:

\[\psi_{nlm}(r,\vartheta,\varphi,t)=R_{nl}(r)\Theta_{lm}(\vartheta)e^{im\varphi}\cdot e^{-i\frac{E}{\hbar}t}.\]

Với hàm \(\Theta_{lm}(\vartheta)\) là nghiệm của phương trình vi phân đã khảo sát kĩ ở bài “Moment động lượng”:

\[\frac{1}{\sin\vartheta}\frac{\partial}{\partial\vartheta}\left(\sin\vartheta\frac{\partial}{\partial\vartheta}\right)\Theta(\vartheta)-\dfrac{m^2}{\sin^2\vartheta}\Theta(\vartheta)=-\frac{M^2}{\hbar^2}\Theta(\vartheta).\]

Còn hàm bán kính \(R_{nl}(r)\) là nghiệm của phương trình vi phân:

\[-\frac{\hbar^2}{2m}\frac{1}{r^2}\frac{\partial}{\partial r}\left(r^2\frac{\partial}{\partial r}\right)R_{nl}(r)+\frac{\hbar^2l(l+1)}{2m_er^2}R_{nl}(r)+U(r)R_{nl}(r)=ER_{nl}(r).\]

\(n\) được gọi là số lượng tử chính. \(l=0,1,2,\ldots n-1\) là một số nguyên không âm, ảnh hưởng đến moment động lượng \(M=\hbar\sqrt{l(l+1)}\), nên còn gọi là số lượng tử quỹ đạo.

\(m=-l,-l+1,\ldots,-1,0,1,\ldots,l-1,l\) ảnh hưởng đến giá trị của moment động lượng theo một trục cố định: \(M_z=m\hbar\). Moment theo trục \(z\) này lại ảnh hưởng đến moment từ khi hệ nguyên tử đặt trong từ trường, nên \(m\) còn có tên gọi là số lượng tử từ.

Bảng dưới liệt kê một cách hệ thống một vài trạng thái có mức năng lượng thấp nhất.

Bảng 1: Liệt kê trạng thái nguyên tử hidro theo các số lượng tử

Code 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
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
function R_hidro_general
% Created by Tran Hai Cat
% Lecturer in physics, quantum mechanics, computing techniques
% 2019.04.27
clc;
clear variables
close all

global h E k epsilon0 q_element m_e l

%% CONSTANTS
h = 1.054e-34;
m_e = 9.1095e-31;
q_element = 1.6e-19;
Angstrom = 1e-10;
k = 2*m_e/h/h;
epsilon0 = 8.854e-12;
r_bohr = 0.529e-10;

%% INPUT PARAMETRS
n = 2;
l = 1;
E = -m_e*q_element^4/32/pi^2/epsilon0^2/h^2/n^2;

% E = -1.5*q_element;

%% DATA PROCESSING
Emin = -20*q_element;
Emax = 20*q_element;
rmin = -5*r_bohr*n^2;
rmax = 5*r_bohr*n^2;
rmin_cal = 0.1*r_bohr;
rmax_cal = 2*rmax;

%% CALCULATION
u0 = 1e-10;
u10 = 0;

[r,U] = ode45(@Schrodinger,linspace(rmax_cal,rmin_cal,1000),[u0 u10]);
u = U(:,1);

r = [-r;flipud(r)];
u = [-u;flipud(u)];
R = u./r;
P = 4*pi*r.^2.*R.^2;
norm = abs(trapz(r,P));
R = R/sqrt(norm);
u = R.*r;
P = P/norm;

u_draw = u*3e-4;
R_draw = R*2e-14;
P_draw = P*3e-9;

%% FIGURE
r_plot = [linspace(-rmax,-0.5*Angstrom,1000),...
linspace(0.5*Angstrom,rmax,1000)];
U_plot = U_hydro(abs(r_plot));

figure('name','R-function hydrogen atom',...
'color','black','numbertitle','off')
hold on
line([0,0],[Emin/q_element,Emax/q_element],'color','w');
patch([r/r_bohr,fliplr(r/r_bohr)],...
[E/q_element+P_draw,fliplr(E/q_element+zeros(size(P_draw)))],...
[0.6 0.3 0.0]);

plot(r/r_bohr,E/q_element+u_draw,'linewidth',1.2,'color','y');
plot(r/r_bohr,E/q_element+R_draw,'linewidth',1.2,'color','m');
plot(r_plot/r_bohr,U_plot/q_element,'linewidth',2,'color','w');
line([rmin/r_bohr,rmax/r_bohr],[E/q_element,E/q_element],...
'linewidth',1.5,'color','b');

title(sprintf('E = %0.3f [eV]',E/q_element),'color','w');
axis([rmin/r_bohr rmax/r_bohr Emin/q_element Emax/q_element]);
set(gca,'color','k','xcolor','w','ycolor','w')
xlabel('r/r_{bohr}');
ylabel('E, U [eV]');
legend(...
'\color{white} r=0',...
'\color{white} 4\pi r^2R^2(r)',...
'\color{white} u(r)',...
'\color{white} R(r)',...
'\color{white} U(r)',...
'\color{white} E'...
);

function dy = Schrodinger(r,y)
%% Schrodinger Equation
global h E k m_e l
dy = zeros(2,1);
dy(1) = y(2);
dy(2) = -k*(E-h*h*l*(l+1)/2/m_e./r./r-U_hydro(r))*y(1);

function y = U_hydro(r)
%% Potential Energy
global epsilon0 q_element
lenr = length(r);
y = zeros(1,lenr);
for i = 1:lenr
y(i) = -1/4/pi/epsilon0*q_element*q_element/r(i);
end