Dao động tử điều hoà - Phần 2

Chúng ta đã tính được hàm sóng của các trạng thái dừng trong hố thế parabol, tương ứng với phổ năng lượng rời rạc:

\begin{equation}
E_n=\hbar\omega_{classic}\left(n-\frac{1}{2}\right),\qquad n=1,2,3,\ldots
\label{eq:phonangluong_hoparabol}
\end{equation}

Trên mỗi mức năng lượng xác định đó, sóng chỉ dao động tại chỗ. Bây giờ, hãy khảo sát một bó sóng hình chuông, đặc trưng cho một hạt đang chuyển động với năng lượng E_0 và xung lượng p_0=\sqrt{2mE_0}. Tại thời điểm ban đầu t=0 sóng có dạng hàm:

\begin{equation}
\psi(x,0)=Ae^{-x^2/4\sigma_x^2}e^{i(\frac{p_0}{\hbar}x-\frac{E_0}{\hbar}0)}.
\label{eq:bosong_t_0}
\end{equation}

Bó sóng \eqref{eq:bosong_t_0} với năng lượng E_0=E_5 được diễn tả như hình 1, với độ bất định vị trí \sigma_x=0.5\,\mathrm{A}. Bó sóng này tượng trưng cho một hạt có năng lượng ngang với mức năng lượng của trạng thái dừng bậc 5. Mật độ của hạt lúc t=0

\begin{equation}
\psi(x,0)^*\psi(x,0)=Ae^{-x^2/2\sigma_x^2}
\label{eq:matdo_t_0}
\end{equation}

có dạng của phân bố Gauss với độ lệch chuẩn bằng \sigma_x, diễn tả qua đường đứt nét trên hình 1. Như vậy, hàm sóng \eqref{eq:bosong_t_0} diễn tả một "đám mây" hạt mà có đến 70\,\% khối lượng của nó hội tụ quanh vị trí x=x_0 trong vòng bán kính \sigma_x.

Hình 1: Bó sóng tại thời điểm t=0

Để biết được sự vận động của bó sóng theo thời gian, ta cần phân tích bó sóng thành sự chồng chập của các trạng thái dừng:

\begin{equation}
\psi(x,0)=\sum_n{C_n\Psi_n(x)},
\label{eq:bosong_parabol_t_0_khaitrien}
\end{equation}

với \Psi_n(x) là nghiệm bậc n của phương trình Schrodinger:

-\frac{\hbar^2}{2m}\frac{\partial^2}{\partial x^2}\Psi(x)+U(x)\Psi(x)=E\Psi(x).

Hệ số C_n trong khai triển \eqref{eq:bosong_parabol_t_0_khaitrien} tính toán qua phép biến đổi Fourier:

\begin{equation}
C_n=\int\limits_{-\infty}^{\infty}{\!\Psi_n(x)^*\psi(x,0)\,dx}.
\label{eq:biendoi_Fourier}
\end{equation}

Kết quả tính toán cho C_n thể hiện qua đồ thị hình 2. Đồ thị biên độ - amplitude thể hiện tỉ lệ góp mặt của mỗi trạng thái dừng. Ta thấy trạng thái dừng bậc 5 đóng vai trò chủ đạo nhất. Điều đó hoàn toàn hợp lý khi chúng ta đang mô phỏng bó sóng đặc trưng cho một hạt có năng lượng E_0=E_5.

Hình 2: Phổ phân tích hàm sóng ra các trạng thái dừng

Sau khi có được bộ hệ số C_n, ta làm phép tổng hợp và cho các thành phần sóng vận động theo thời gian:

\begin{equation}
\psi(x,t)=\sum_n{C_n\Psi_n(x)}e^{i(E_n/\hbar)t}.
\label{eq:bosong_parabol_t_tonghop}
\end{equation}

Kết quả sẽ có một sóng vận động như biểu diễn trong video dưới. Đó không phải là một sóng dừng, mà mang hình ảnh của một hạt dao động điều hoà qua lại giữa hai bờ hố thế, giới hạn bởi mức năng lượng toàn phần, có tần số dao động đúng bằng \omega_{classic}. Đúng như trên thực tế, orbitan di chuyển chứ không rung động tại chỗ như các trạng thái dừng. Điều đó xảy ra bởi vì vi hạt đồng thời nằm trên nhiều mức năng lượng, tương ứng với nhiều trạng thái song song cùng một lúc.

 

Code chương trình Matlab

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
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
function wave_packet_parabol_potential_well
% Created by Tran Hai Cat
% 2019.04.15
clc;
clear variables
close all

global h m k k_parabol E

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

%% INPUT PARAMETRS FOR USER
omega0 = sqrt(k_parabol/m);
E0 = h*omega0*(7-1/2);

N = 30;
k_parabol = 400;
x0 = 0;
delta_x = 0.7*Angstrom;

Emin = -20*q_element;
Emax = 200*q_element;
xmin = -6*Angstrom;
xmax = 6*Angstrom;
x1 = -10*Angstrom;
x2 = 10*Angstrom;
x_inter = 0.5769*Angstrom;

%% CALCULATION
omega_array = zeros(1,N);
E_array = zeros(1,N);
psi = 1e-10;
psi1 = 0;
for i = 1:N
omega_classic = sqrt(k_parabol/m);
E_array(i) = h*omega_classic*(i-1/2);
omega_array(i) = E_array(i)/h;

E = E_array(i);
[x_left,psi_left] = ode45(@Schrodinger_equation,linspace(x1,x_inter,1500),[psi psi1]);
[x_right,psi_right] = ode45(@Schrodinger_equation,linspace(x2,x_inter,1500),[psi psi1]);
psi_right = sign(psi_left(end,2)*psi_right(end,2))*psi_right;
X = [x_left;flipud(x_right)];
Psi_E{i} = [psi_left(:,1);flipud(psi_right(:,1))];
S_norm = trapz(X,Psi_E{i}.*Psi_E{i});
Psi_E{i} = Psi_E{i}/sqrt(S_norm);
end

C = packet_analys(E0,x0,delta_x,N,X,Psi_E);

Psi_synthes = wavepacket_synthes(C,omega_array,Psi_E,0);
Psi_synthes_max = abs(max(Psi_synthes));
Psi_draw = 0.25*Psi_synthes/Psi_synthes_max*Emax/q_element;
psi_psi = 0.3*(conj(Psi_synthes).*Psi_synthes)...
/Psi_synthes_max/Psi_synthes_max*Emax/q_element;

Amp = sqrt(2*E0/k_parabol);
x_particle = 0;
y_particle = U_parabol(k_parabol,x_particle);

%% FIGURE
figure('name','Spectrum of Stationary States',...
'color','w','numbertitle','off');
subplot(2,1,1);
stem(abs(C),'fill','-','color','blue');
ylabel('Amplitude');

subplot(2,1,2);
stem(angle(C),'fill','-','color','blue');
xlabel('n (harmonic)');
ylabel('Phase');

x_plot = linspace(xmin,xmax,1000);
U_plot = U_parabol(k_parabol,x_plot);

figure('name','Harmonic oscillator',...
'color','black','numbertitle','off');
hold on
plot(x_plot/Angstrom,U_plot/q_element,'linewidth',2,'color','w');
line([xmin/Angstrom,xmax/Angstrom],[E0/q_element,E0/q_element],...
'linewidth',2,'color','b');
line_psi = plot(X/Angstrom,E0/q_element+real(Psi_draw),'linewidth',1,'color','y');
line_psi_psi = plot(X/Angstrom,E0/q_element+psi_psi,'linewidth',2);
draw_partile = plot(x_particle/Angstrom,y_particle/q_element,...
'ro','MarkerSize',10,'markerfacecolor','r');
title(sprintf('E0 = %0.3f [eV]',E0/q_element),'color','w');
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]');

t = 0;
dt = 1e-19;
while 1
t = t+dt;
Psi_synthes = wavepacket_synthes(C,omega_array,Psi_E,t);
Psi_draw = 0.25*Psi_synthes/Psi_synthes_max*Emax/q_element;
psi_psi = 0.3*(conj(Psi_synthes).*Psi_synthes)...
/Psi_synthes_max/Psi_synthes_max*Emax/q_element;
x_particle = Amp*sin(omega0*t);
y_particle = U_parabol(k_parabol,x_particle);

set(line_psi,'ydata',E0/q_element+real(Psi_draw));
set(line_psi_psi,'ydata',E0/q_element+real(psi_psi));
set(draw_partile,'xdata',x_particle/Angstrom);
set(draw_partile,'ydata',y_particle/q_element);
pause(0.001);
end

function dy = Schrodinger_equation(x,y)
%% Schrodinger Equation
global E k k_parabol
dy = zeros(2,1);
dy(1) = y(2);
dy(2) = -k*(E-U_parabol(k_parabol,x))*y(1);

function y = U_parabol(k_parabol,x)
%% Potential Energy
lenx = length(x);
y = zeros(1,lenx);
for i = 1:lenx
y(i) = 0.5*k_parabol*x(i)*x(i);
end

function psi = wavepacket_Gauss(E,x0,delta_x,x)
%% wave packet with Gauss's distribution (normal distribution)
global h m
k0 = sqrt(2*m*E)/h;
psi = 300*exp(-(x-x0).^2/4/delta_x.^2).*exp(1i*k0*(x-x0));

function C = packet_analys(E,x0,delta_x,N,X,Psi_E)
%% Calculate koefficients of Stationary State
psi0 = wavepacket_Gauss(E,x0,delta_x,X);
C = zeros(1,N);
for i = 1:N
C(i) = trapz(X,conj(Psi_E{i}).*psi0);
end

function psi = wavepacket_synthes(C,omega_array,Psi_E,t)
%% Synthesys initial wave packet from koefficients of Stationary State
N = length(C);
psi = 0;
for i = 1:N
psi = psi+C(i)*Psi_E{i}*exp(-1i*omega_array(i)*t);
end