Issue in getting correct probabilities for Landau Zener Problem
6 views (last 30 days)
Show older comments
Aryaman
on 9 Apr 2024
Edited: David Goodmanson
on 20 Apr 2024
I am trying to solve the Landau Zener problem where i have the equations
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1664956/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1664961/image.png)
I have written a code using Finite Difference as given below. The result is wrong. The correct result is ![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1664966/image.jpeg)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1664966/image.jpeg)
Can anybody help find the error?
clc;
close all;
%-----------------------------------------------------------------------%
v = 0.1; % sweep velocity
delta =4; % coupling parameter between the ground and excited states
t_max = 1000;
dt = 0.1; % time step
t = -t_max:dt:t_max; % defining time range with step size dt
N = length(t);
b=delta/2;
hbar=1;
%---- initial conditions-------------
cg(1:N+1)=1; % initial ground state
ce(1:N+1)=0; % initial excited state
for i=1:N % finite difference loop for advancing in time
a=v*t(i)/2;
cg(i+1)=(dt*a*t(i)/1i*hbar)*cg(i)+cg(i)+ dt*b*ce(i)/1i*hbar;
ce(i+1)=-(dt*a*t(i)/1i*hbar)*ce(i)+ce(i) + dt*b*cg(i)/1i*hbar;
end
% plotting the probabilities in ground and excited states
plot(t,abs(cg(1:N).*cg(1:N)),'linewidth',2)
hold on
plot(t,abs(ce(1:N).*ce(1:N)),'r','linewidth',2)
xlabel('Time');
ylabel('Probability');
legend('|c1|^2', '|c2|^2')
2 Comments
Accepted Answer
David Goodmanson
on 9 Apr 2024
Edited: David Goodmanson
on 20 Apr 2024
Hello Aryaman,
I believe your step sizes are too large, but the first issue is, in the for loop you define a=v*t(i)/2; but then in the next line you multiply a by t(i) again, so now the Hamiltonian has factors of t^2. That's not helping.
It's good to program up a solver on your own, but in this case the solutions for large positive and negative t are going like (for the posted equations at the very top of the question)
e^(i*(a/2)*t^2) and e^(-i*(a/2)*t^2)
which oscillate increasigly quickly the further out you go. In this case it's probably best, at least initially, to let a Matlab ode solver pick dt, which it does dynamically. In the following code I arbitrarily picked factors of 1 in the Hamiltonian:
f = @(t,x) -1i*[t 1;1 -t]*x;
[t y] = ode45(f,[-20 20],[1 0]);
fig(1)
plot(t,abs(y).^2)
grid on
which produces the plot below for each of cg and ce. I imagine you could reproduce this reasonably well with the finite difference loop, although the Euler forward difference method you are implementing means there will need to be a lot of points.
From the posted equations it appears that there are two independent parameters, a and b (or three if b is complex) , but if the problem is rescaled there is really only one (see below).
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1665041/image.png)
The posted equations are, with complex b represented as b e^( i phi) and b positive real,
i hbar d/dt c1 = at c1 + b e^(i phi) c2
i hbar d/dt c2 = b (e^-i phi)* c1 -at c2
Absorbing the phase of b into c2 with c2 --> e^(-i phi) c2 results in
i hbar d/dt c1 = at c1 + b c2
i hbar d/dt c2 = b c1 -at c2
(this phase must be kept track of in order to readjust c2 after the solution is found!).
The far more important result comes from rescaling to go to dimensionless time. There are two physical units here, time and energy, and three parameters, hbar ~ Et, a ~ E/t and b ~ E. Let t = t0*u, where t0 = sqrt(hbar/a) and u is dimensionless. Plugging that in results in
i d/du c1 = u c1 + B c2
i d/du c2 = B c1 - u c2
where B = b/sqrt(hbar a)
and B is dimensionless. This agrees with the Buckham pi theorem: 2 units, 3 variables implies 1 dimensionless constant. (You are free to do some implicit scaling as well by choosing hbar = 1 as you did). B is similar to the Reynolds number in that it determines the characteristic behavior of the solution. A function of B determines the probabilities as u --> infinity.
4 Comments
Torsten
on 10 Apr 2024
If hbar = h' and h*h' = 1, then you were right multiplying your equations with hbar.
If hbar = h, you were wrong.
I interpreted h = hbar.
More Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!