Issue in getting correct probabilities for Landau Zener Problem

I am trying to solve the Landau Zener problem where i have the equations

I have written a code using Finite Difference as given below. The result is wrong. The correct result is

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')

[Torsten]

Why don't you use MATLAB's ode45 or ode15s ?

[Aryaman]

I have been asked to do it using finite difference

Answer 1

[David Goodmanson]

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).

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.

[Aryaman]

Thank you for your response. I rewrote the equation and rewrote the the time step.

clc;
close all;

%-----------------------------------------------------------------------%

%N=1000; % no of time grid points
%t=linspace(0,0.25e-6,N); % time array in micro seconds
%dt=t(2)-t(1); % step size in time
v = 0.1;                    % sweep velocity
delta =1;                % coupling parameter between the ground and excited states
t_max = 20;
dt = 0.000001;                 % 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=2:N % finite difference loop for advancing in time
    cg(i+1)=(dt*v*t(i)/1i*hbar)*cg(i-1)+cg(i-1)+ dt*b*ce(i)/1i*hbar;
    ce(i+1)=-(dt*v*t(i)/1i*hbar)*ce(i-1)+ce(i-1) + 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('|cg|^2', '|ce|^2');

[Torsten]

I reduced dt and divided by hbar instead of multiplying.

The results from Explicit Euler and ode45 are almost equal.

clc;
close all;

%-----------------------------------------------------------------------%

%N=1000; % no of time grid points
%t=linspace(0,0.25e-6,N); % time array in micro seconds
%dt=t(2)-t(1); % step size in time
v = 0.001;                    % sweep velocity
delta =1;                % coupling parameter between the ground and excited states
t_max = 10;
dt = 0.001;                 % 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
    cg(i+1)=(dt*v*t(i)/(1i*hbar))*cg(i)+cg(i)+ dt*b*ce(i)/(1i*hbar);
    ce(i+1)=-(dt*v*t(i)/(1i*hbar))*ce(i)+ce(i) + dt*b*cg(i)/(1i*hbar);
   
end

% plotting the probabilities in ground and excited states
hold on
plot(t,abs(cg(1:N).*cg(1:N)),'linewidth',2)
plot(t,abs(ce(1:N).*ce(1:N)),'r','linewidth',2)
xlabel('Time');
ylabel('Probability');
legend('|cg|^2', '|ce|^2');

%professional solver to compare results
f = @(t,x) -1i*[0.001/1*t 0.5/1*1;0.5/1*1 -0.001/1*t]*x;
[t y] = ode45(f,[-10 10],[1 0]);
plot(t,abs(y).^2)
grid on
hold off

[Aryaman]

Thanks

[Torsten]

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.