Bó sóng với rào thế bậc thang

Trong bài giảng về sự tương tác của sóng de-Broglie với rào thế bậc thang, có thể hiểu rằng đó là tương tác giữa một chùm hạt đồng nhất lên rào thế. Để hiểu rõ ý nghĩa của tương tác này, ta sẽ đi xây dựng mô hình bó sóng với rào thế bậc thang, đặc trưng cho một hạt lao về phía rào thế.

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\) và xung lượng \(p=\sqrt{2mE}\). Tại thời điểm ban đầu \(t=0\) sóng có dạng hàm:

\[\psi(x,0)=Ae^{-x^2/4\sigma_x^2}e^{i(\frac{p}{\hbar}x-\frac{E}{\hbar}0)}.\label{eq:1}\tag{1}\]

Hàm sóng \eqref{eq:1} được diễn tả như hình 1, với độ bất định vị trí \(\sigma_x=10\,\mathrm{A}\). Mật độ của hạt lúc \(t=0\)

\[\psi(x,0)^*\psi(x,0)=Ae^{-x^2/2\sigma_x^2}\]

có dạng của phân bố Gauss với độ lệch chuẩn bằng \(\sigma_x\), diễn tả qua đường màu cam trên hình 1. Như vậy, hàm sóng \eqref{eq:1} 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\).

bó sóng với rào thế bậc thang
Hình 1: Bó sóng trước khi va vào rào thế

Để biết được bó sóng sẽ di chuyển như thế nào, ta cần phân tích bó sóng thành sự chồng chập của nhiều sóng de-Broglie tại thời điểm \(t=0\), bởi vì sóng de-Broglie có quy luật vận động rõ ràng. Thực vậy, qua bài Sóng de-Broglie với rào thế bậc thang, sóng được hình thành từ sóng tới, sóng phản xạ và sóng truyền qua:

\[\begin{cases}\psi_{in}(x,t)=e^{i(\frac{p_1}{\hbar}x-\frac{E}{\hbar}t)}\\ \psi_{ref}(x,t)=\dfrac{p_1-p_2}{p_1+p_2}e^{i(-\frac{p_1}{\hbar}x-\frac{E}{\hbar}t)},\\ \psi_{tran}(x,t)=\dfrac{2p_1}{p_1+p_2}e^{i(\frac{p_2}{\hbar}x-\frac{E}{\hbar}t)},\end{cases}\label{eq:2}\tag{2}\]

với \(p_1=\sqrt{2mE}, p_2=\sqrt{2m(E-U_0)}\). Ta sẽ khai triển bó sóng \eqref{eq:1} tại thời điểm \(t=0\) thành sự chồng chập của nhiều sóng de-Broglie dạng \eqref{eq:2}:

\[\psi(x,0)=\int\limits_{-\infty}^{\infty}{C(p)\left[\psi_{in}(x,0)+\psi_{ref}(x,0)\right]\,dp},\]

Phép phân tích này chỉ bao gồm sóng tới \(\psi_{in}(x,0)\) và sóng phản xạ \(\psi_{ref}(x,0)\), bởi ban đầu khi \(t=0\) bó sóng hoàn toàn nằm trong phạm vi của sóng tới và sóng phản xạ, tức vùng 1. Các hệ số \(C(p)\) được tính toán theo phép biến đổi Fourier:

\[\begin{aligned}C(p)&=\frac{1}{2\pi\hbar}\int\limits_{-\infty}^{\infty}{\left[\psi_{in}^(x,0)+\psi_{ref}^(x,0)\right]\psi(x,0)\,dx},\\&=\frac{1}{2\pi\hbar}\int\limits_{-\infty}^{\infty}{\left[e^{-i\frac{p_1}{\hbar}x}+\frac{p_1-p_2}{p_1+p_2}e^{i\frac{p_1}{\hbar}x}\right]\psi(x,0)\,dx}.\end{aligned}\label{eq:3}\tag{3}\]

Kết quả tính toán cho \(C(p)\) thể hiện trên hình 2:

Hình 2: Phân bố biên độ \(C(p)\) của các thành tố de-Broglie tạo nên bó sóng

Trạng thái của bó sóng thời điểm \(t\) bất kì là sự tổng hợp các sóng de-Broglie trở lại:

\[\begin{cases}\psi(x,t)=\int\limits_{-\infty}^{\infty}{C(p)\left[\psi_{in}(x,t)+\psi_{ref}(x,t)\right]dp} & \mathrm{khi}\ x<0,\\ \psi(x,t)=\int\limits_{-\infty}^{\infty}{C(p)\psi_{tran}(x,t)\,dp} & \mathrm{khi}\ x>0.\end{cases}\]

với các hệ số \(C(p)\) đã tính được từ phép biến đổi \eqref{eq:3}. Sử dụng code chương trình Matlab bên dưới, ta có thể thực hiện chuỗi tính toán trên bằng máy tính.

Video sau đây biểu diễn kết quả tính toán cho trường hợp hạt có năng lượng cao hơn rào thế: \(E>U_0\). Theo cơ học cổ điển, hạt sẽ di chuyển dễ dàng qua rào thế và đi tiếp. Còn ở đây, trong thế giới lượng tử, hạt chỉ đi qua một phần, còn một phần bị phản xạ ngược!

Video tiếp theo lại minh hoạ kết quả tính toán cho trường hợp hạt có năng lượng \(E\) thấp hơn chiều cao rào thế. Theo cơ học cổ điển, hạt sẽ bị dội ngược lại do không đủ năng lượng để leo lên. Trong cơ học lượng tử cũng xảy ra tình huống tương tự: hạt bị phản xạ toàn phần!

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
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
function wavepacket_step_potential_barrier
% Created by Tran Hai Cat
% 2018.05.01
% 2019.03.14 edited by Cat
clc;
clear variables
close all

global j h m dp

%% CONSTANTS
h = 1.054e-34;
m = 9.1095e-31;
q_element = 1.6e-19;
Angstrom = 1e-10;
j = sqrt(-1);

%% INPUT PARAMETRS
E0 = 5*q_element;
x_step = 0*Angstrom;
U0 = 10*q_element;

xmin = -200*Angstrom;
xmax = 200*Angstrom;

x0 = -150*Angstrom; % Initial possition of wave packet
delta_x = 10*Angstrom; % standart deviation (radius of wave packet)

Np = 20;
dp = 2*h/delta_x/Np;

%% DATA PROCESSING
x = linspace(xmin,xmax,2000);
w = E0/h;
T = 2*pi/w;
dt = T/20;
tmax = 1*T;
t = 0;

%% ANALYS WAVE PACKET INTO DE BROGLIE
p0 = sqrt(2*m*E0);
[p,k1,k2,B1,S2,w,C] = packet_analys(E0,U0,x0,delta_x,x_step,Np); %packet_analys(E,x,x0,delta_x,Np);
psi = wavepacket_synthes(k1,k2,B1,S2,w,p,C,x_step,x,t);
E = p.^2/2/m;

%% FIGURE
figure('name','Analys koefficients of Fourier series',...
'color','white','numbertitle','off')
stem(p/p0,abs(C));
xlabel('p/p_0');
ylabel('C');

figure('name','Wavepaket Barrier','color','black','numbertitle','off');
hold on
plot(x/Angstrom,Energy(E0,x)/q_element,'b','linewidth',2);
for i = 1:length(E)
line([xmin/Angstrom,xmax/Angstrom],[E(i)/q_element,E(i)/q_element],...
'color','b','linewidth',0.5);
end
plot(x/Angstrom,U_step(U0,x_step,x)/q_element,'w','linewidth',2);

line_psi = plot(x/Angstrom,10+real(psi),'linewidth',1,'color','y');
line_psi_psi = plot(x/Angstrom,psi.*conj(psi),'linewidth',1,'color',[1,0.5,0]);
axis([xmin/Angstrom xmax/Angstrom -0 2*U0/q_element]);
set(gca,'color','k','xcolor','w','ycolor','w')
xlabel('x [Angstrom]');
ylabel('E, U [eV]');

%% ANIMATION
while 1
t = t+dt;
psi = wavepacket_synthes(k1,k2,B1,S2,w,p,C,x_step,x,t);
set(line_psi_psi,'ydata',psi.*conj(psi));
set(line_psi,'ydata',U0/q_element+real(psi));
pause(0.001);
end

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

function y = U_step(U0,x_step,x)
%% potential energy
lenx = length(x);
y = zeros(1,lenx);
for i = 1:lenx
if(x(i)<=x_step)
y(i) = 0;
else
y(i) = U0;
end
end

function y = Energy(E,x)
%% particle's energy
y = E*ones(size(x));

function psi = de_Broglie_step(k1,k2,B1,S2,w,x_step,x,t)
%% de Broglie's wave on step barrier
global j
psi = zeros(size(x));

for i = 1:length(x)
if(x(i)<=x_step)
psi(i) = exp(j*k1*x(i)-j*w*t)+B1*exp(-j*k1*x(i)-j*w*t);
else
psi(i) = S2*exp(j*k2*x(i)-j*w*t);
end
end

function [p,k1,k2,B1,S2,w,C] = packet_analys(E,U0,x0,delta_x,x_step,Np)
%% Calculate koefficients of Fourier's series
global j h m dp
x = linspace(-200e-10,200e-10,5000);

p0 = sqrt(2*m*E);
psi0 = wavepacket_Gauss(E,x0,delta_x,x);
p = zeros(1,2*Np+1);
k1 = zeros(1,2*Np+1);
k2 = zeros(1,2*Np+1);
B1 = zeros(1,2*Np+1);
S2 = zeros(1,2*Np+1);
w = zeros(1,2*Np+1);
C = zeros(1,2*Np+1);

p(Np+1) = p0;
E0 = p0*p0/2/m;

w(Np+1) = E0/h;
k1(Np+1) = sqrt(2*m*E0)/h;
k2(Np+1) = sqrt(2*m*(E0-U0))/h;
B1(Np+1) = (k1(Np+1)-k2(Np+1))/(k1(Np+1)+k2(Np+1))*exp(2*j*k1(Np+1)*x_step);
S2(Np+1) = 2*k1(Np+1)/(k1(Np+1)+k2(Np+1))*exp(j*(k1(Np+1)-k2(Np+1))*x_step);

C(Np+1) = (1/2/pi/h)...
*trapz(x,conj(de_Broglie_step(k1(Np+1),k2(Np+1),B1(Np+1),S2(Np+1),w(Np+1),x_step,x,0)).*psi0);
for ip = 1:Np
p(Np+1-ip) = p0-ip*dp;
E = p(Np+1-ip)*p(Np+1-ip)/2/m;

w(Np+1-ip) = E/h;
k1(Np+1-ip) = sqrt(2*m*E)/h;
k2(Np+1-ip) = sqrt(2*m*(E-U0))/h;
B1(Np+1-ip) = (k1(Np+1-ip)-k2(Np+1-ip))/(k1(Np+1-ip)+k2(Np+1-ip))*exp(2*j*k1(Np+1-ip)*x_step);
S2(Np+1-ip) = 2*k1(Np+1-ip)/(k1(Np+1-ip)+k2(Np+1-ip))*exp(j*(k1(Np+1-ip)-k2(Np+1-ip))*x_step);

C(Np+1-ip) = (1/2/pi/h)...
*trapz(x,conj(de_Broglie_step(k1(Np+1-ip),k2(Np+1-ip),B1(Np+1-ip),S2(Np+1-ip),w(Np+1-ip),x_step,x,0)).*psi0);

p(Np+1+ip) = p0+ip*dp;
E = p(Np+1+ip)*p(Np+1+ip)/2/m;

w(Np+1+ip) = E/h;
k1(Np+1+ip) = sqrt(2*m*E)/h;
k2(Np+1+ip) = sqrt(2*m*(E-U0))/h;
B1(Np+1+ip) = (k1(Np+1+ip)-k2(Np+1+ip))/(k1(Np+1+ip)+k2(Np+1+ip))*exp(2*j*k1(Np+1+ip)*x_step);
S2(Np+1+ip) = 2*k1(Np+1+ip)/(k1(Np+1+ip)+k2(Np+1+ip))*exp(j*(k1(Np+1+ip)-k2(Np+1+ip))*x_step);

C(Np+1+ip) = (1/2/pi/h)...
*trapz(x,conj(de_Broglie_step(k1(Np+1+ip),k2(Np+1+ip),B1(Np+1+ip),S2(Np+1+ip),w(Np+1+ip),x_step,x,0)).*psi0);
end

function psi = wavepacket_synthes(k1,k2,B1,S2,w,p,C,x_step,x,t)
%% Synthesys initial wave packet from koefficients of Fourier's series
global j
lenp = length(p);
lenx = length(x);
psi = zeros(size(x));
for ik = 1:lenp-1
for i = 1:lenx
if(x(i)<=x_step)
psi(i) = psi(i)...
+C(ik)*(exp(j*(k1(ik)*x(i)-w(ik)*t))...
+B1(ik)*exp(j*(-k1(ik)*x(i)-w(ik)*t)))...
*(p(ik+1)-p(ik));
else
psi(i) = psi(i)...
+C(ik)*S2(ik)*exp(j*(k2(ik)*x(i)-w(ik)*t))...
*(p(ik+1)-p(ik));
end
end
end