Quit[]

adiabatic problem 1: standard LZ probability.

The goal of this section is to solve the standard LZ problem, i.e. a time-dependent two level system with Hamiltonian

H={{-α*t,β},{β,α*t}}

solvesys does this numerical using NDsolve, subject to the initial condition that the state is initially all in channel 1. 
LZ gives the LZ probability of a diabatic transition. 
The other functions (markov, interino, infunc) are used to get the Markov approximation used in Glasbrenner and Schleich. 
Runall runs all of these and makes plots of Re, Im parts of the two channel solutions (from NDsolve), the probabilities as a function of time vs Landau Zener, and a comparison of the same with the Markov approximation

solvesys[α_,β_,t0_,tf_]:={ψ1[t],ψ2[t]}/.NDSolve[{I*D[{ψ1[t],ψ2[t]},t]=={{-α*t,β},{β,α*t}}.{ψ1[t],ψ2[t]},ψ1[t0]==1,ψ2[t0]==0.0},{ψ1[t],ψ2[t]},{t,t0,tf}][[1]]

plotpots[α_,β_,t0_,tf_]:=Module[{t,ev},
ev=Eigenvalues[{{-α*t,β},{β,α*t}}];
Plot[{α*t,-α*t,ev[[1]],ev[[2]]},{t,t0,tf}]
]

LZ[α_,β_]:=Exp[-2*Pi*β^2/2/α]

infunc[α_,β_,t_,t0_]:=Integrate[Exp[I*α*tp^2/2],{tp,t0,t}]

markov[α_,β_,t0_,tf_]:=Module[{bhz,bhzi},
bhz=Table[{t,infunc[α,β,t,t0]},{t,t0,tf,(tf-t0)/10000}];
bhzi=Interpolation[bhz];
a[t]/.NDSolve[{D[a[t],t]==-β^2*Exp[-I*α*t^2/2]*bhzi[t]*a[t],a[t0+1]==1},a[t],{t,t0+1,tf}][[1]]
]

Integrate[Exp[I*α*q^2],{q,t0,p}]

interino[α_,β_,t0_,t_]:=NIntegrate[Exp[-I*α*p^2]*((((((-1))^(1/4)) (Sqrt[π]) (-Erf[(((-1))^(3/4)) p (Sqrt[α])]+Erf[(((-1))^(3/4)) t0 (Sqrt[α])]))/(2 (Sqrt[α])))),{p,t0,t}]

runall[α_,β_,t0_,tf_,t0p_,tfp_]:=Module[{p1,p2,bhz,ffa,tt,tti},
{p1,p2}=solvesys[α,β,t0,tf];
Print[plotpots[α,β,t0p/10,tfp/10]];
Print[Plot[{Re[p1],Im[p1]},{t,t0p,tfp}]];
Print[Plot[{Re[p2],Im[p2]},{t,t0p,tfp}]];
Print[Plot[{Abs[p1]^2,Abs[p2]^2,LZ[α,β],1-LZ[α,β],Abs[p1]^2+Abs[p2]^2},{t,t0p,tfp},PlotRange->All]];
(*Print[ParametricPlot[{Re[p1],Im[p1]},{t,t0,0},PlotRange->{{0,1},{-1,0}}]]*);
tt=ParallelTable[{t,Exp[-β^2*interino[α,β,t0p,t]//N]},{t,t0p,tfp,(tfp-t0p)/1000}];
tti=Interpolation[tt];
Print[Plot[{Abs[p1]^2,Abs[p2]^2,LZ[α,β],1-LZ[α,β],Abs[p1]^2+Abs[p2]^2,Abs[tti[t]]^2},{t,t0p,tfp},PlotRange->All]];
]

runall[1,0.5,-100,100,-20,20]

adiabatic problem 2: a problem with two crossings

Now we are going to tackle a harder problem involving the diabatic Hamiltonian

H={{0,a},{a,c-b*(r0-R)^2}}

There are two avoided crossings present with identical LZ parameters. This code analyzes this problem using LZ formalism, trying to include also Stueckelberg oscillations.

(*parameters*)

{A,B,CC,R0}={0.1,1,1,10}

H[R_,a_,b_,c_,r0_]:={{0,a},{a,c-b*(r0-R)^2}};

(*dH is used for the Feynman Hellman theorem*)

dH[R_,a_,b_,c_,r0_]=D[H[R1,a,b,c,r0],R1]/.R1->R;

(*numerically calculate the adiabatic PECs and adiabatic states*)

{ev1,ev2}=Eigenvalues[H[R,a,b,c,r0]];
numevec=Table[{R1,Eigensystem[H[R1,A,B,CC,R0]-IdentityMatrix[2]*10]},{R1,7,12,0.0021}];
r=numevec[[All,1]];
evec[j_,c_]:=Interpolation[Table[{r[[i]],numevec[[i,2,2,j,c]]},{i,1,Length[numevec]}]]
evm=Table[Table[evec[j,c],{j,1,2}],{c,1,2}];

(*plot adiabatic and diabatic curves*)

Show[{ListLinePlot[{Table[{r[[i]],numevec[[i,2,1,2]]+10},{i,1,Length[numevec]}],Table[{r[[i]],numevec[[i,2,1,1]]+10},{i,1,Length[numevec]}]}],Plot[{H[R,A,B,CC,R0][[1,1]],H[R,A,B,CC,R0][[2,2]]},{R,7,12},PlotRange->All,PlotStyle->{{Black,Dashed}}]}]

(*the crossing points are given by (for p = -1 and p = +1)*)

Rm[a_,b_,c_,r0_,p_]:=((p*(Sqrt[c])+(Sqrt[b]) r0)/(Sqrt[b]))

Rm[A,B,CC,R0,-1]

Rm[A,B,CC,R0,1]

(*the gap between adiabatic energy levels is*)

del[a1_,b1_,c1_,r01_,R1_]:=(ev1-ev2)/.{a->a1,b->b1,c->c1,r0->r01,R->R1}

del[A,B,CC,R0,Rm[A,B,CC,R0,-1]]

Δ=del[A,B,CC,R0,Rm[A,B,CC,R0,1]]

(*the difference in slope of the two potential curves at the crossing is*)
dSlope=dH[Rm[A,B,CC,R0,1],A,B,CC,R0]

(*compute the P matrix two ways: 1 using the Feynman Hellman theorem, and 2 using the definition of the P matrix*)

Pmat1=Table[{r[[i]],numevec[[i,2,2,1]].dH[r[[i]],A,B,CC,R0].numevec[[i,2,2,2]]/(numevec[[i,2,1,2]]-numevec[[i,2,1,1]])},{i,1,Length[r]}];

Pmat2=Table[{r[[i]],numevec[[i,2,2,1]].numevec[[i+1,2,2,2]]/(r[[i+1]]-r[[i]])},{i,1,Length[r]-1}];

Pmat1i=Interpolation[Pmat1];

(*they agree, and additionally. they even match the Lorentzian formula (Eq. 12) from CW Clark Phys. Lett. 70, 1979 very nicely*)

Show[{ListLinePlot[{Pmat1,Pmat2},PlotRange->All,PlotStyle->{{Black,Thick},{Red,Dashed}}],Plot[1/2*(dSlope[[2,2]]/Δ)/(1+(dSlope[[2,2]]/Δ*(R-Rm[A,B,CC,R0,1]))^2),{R,7,12},PlotRange->All,PlotStyle->Green]}]

(*Landau Zener probability two ways: from Clark's Lorentzian formula and from the definition using diabataic potentials*)

LZP[en_,V_,a_,b_,c_,r0_]:=
Module[{v,V11t,V22t,P1,P2,Rmm},
v=Sqrt[2*(en-V)];
Rmm=Rm[a,b,c,r0,-1];
V11t=0;
V22t=(2 b (-Rmm+r0));
P1=Exp[-2*Pi*a^2/v/Abs[(V11t-V22t)]];
P2=Exp[-2*Pi*del[a,b,c,r0,Rmm]/8/v/Pmat1i[Rmm]];
Return[{P1,P2}]
]

LZP[2.5,0,A,B,CC,R0]

Now, let’s get started solving the following question: given a particle incident with energy E from the left, originally traveling on the flat yellow potential, what is the probability that the particle exits on the right in the lower (blue) channel?

ev2s[a1_,b1_,c1_,r01_,R1_]:=(ev2)/.{a->a1,b->b1,c->c1,r0->r01,R->R1}
ev1s[a1_,b1_,c1_,r01_,R1_]:=(ev1)/.{a->a1,b->b1,c->c1,r0->r01,R->R1}

Plot[{ev1s[A,B,CC,R0,R],ev2s[A,B,CC,R0,R]},{R,8,12}]

(*A first step is to incoherently add the two LZ probabilities. The particle can first do a diabatic crossing, then stay adiabatic, or it can do the reverse. THe probability is just the sum of these two options*)

incoherentProb[en_,V_,a_,b_,c_,r0_]:=2*LZP[en,V,a,b,c,r0][[1]]*(1-LZP[en,V,a,b,c,r0][[1]])

(*A more accurate calculation includes the phase accumulated along the two different pathways coherently. For this, we need first the phases accumulated*)

te1[en_]:=NIntegrate[Sqrt[2*(en-ev1s[A,B,CC,R0,R])],{R,Rm[A,B,CC,R0,-1],Rm[A,B,CC,R0,1]}]
te2[en_]:=NIntegrate[Sqrt[2*(en-ev2s[A,B,CC,R0,R])],{R,Rm[A,B,CC,R0,-1],Rm[A,B,CC,R0,1]}]

ListLinePlot[{Table[{en,te1[en]},{en,1.01,6,0.1}],Table[{en,te2[en]},{en,1.01,6,0.1}]}]

(*Now we have three options for including the coherence of these two pathways in the probability. 
Option 1 (T1 below) computes the probability the opposite way: 1 - P_blue curve. To get that probability, we add path one (adiabatic + adiabatic) with phase accumulated along the upper potential, to path two (diabatic + diabatic) with phase accumulated along the lower potential.
Note that this requires subtracting rather than adding the amplitudes. I am not sure how to justify this additional phase just yet. 

Option 2 (T2 below) assumes that the particle takes path one (diabatic then adiabatic) with phase accumulated along the lower potential, and that it takes path two (adiabatic then diabatic) with phase accumulated along the upper potential. 

Option 3 (T3) assumes the same as T1, but including now also the Stokes phase (for now I have guessed the signs in the exponentials - other options seem to give worse results!)
*)

coherentProb[en_,V_,a_,b_,c_,r0_,R1_,R2_]:=
Module[{T1,int1,int2,T2,T,v,Rmm,δ,P2,ϕ,T3},
v=Sqrt[2*(en-V)];
Rmm=Rm[a,b,c,r0,-1];
δ=del[a,b,c,r0,Rmm]/8/v/Pmat1i[Rmm];
P2=Exp[-2*Pi*δ];
ϕ=Pi/4+Arg[Gamma[1-I*δ]]+δ*(Log[δ]-1);
int1=NIntegrate[Sqrt[2*(en-ev1s[a,b,c,r0,R])],{R,R1,R2}];
int2=NIntegrate[Sqrt[2*(en-ev2s[a,b,c,r0,R])],{R,R1,R2}];
T1=Sqrt[LZP[en,V,a,b,c,r0][[1]]]Exp[I*int1]Sqrt[LZP[en,V,a,b,c,r0][[1]]]-Sqrt[1-LZP[en,V,a,b,c,r0][[1]]]*Exp[I*int2]*Sqrt[1-LZP[en,V,a,b,c,r0][[1]]];
T2=Sqrt[LZP[en,V,a,b,c,r0][[1]]]Exp[I*int1]Sqrt[1-LZP[en,V,a,b,c,r0][[1]]]+Sqrt[1-LZP[en,V,a,b,c,r0][[1]]]*Exp[I*int2]*Sqrt[LZP[en,V,a,b,c,r0][[1]]];
T3=Sqrt[LZP[en,V,a,b,c,r0][[1]]]Exp[I*int1+I*ϕ]Sqrt[1-LZP[en,V,a,b,c,r0][[1]]]-Sqrt[1-LZP[en,V,a,b,c,r0][[1]]]*Exp[I*int2-I*ϕ]*Sqrt[LZP[en,V,a,b,c,r0][[1]]];
Return[{1-Abs[T1]^2,Abs[T2]^2,Abs[T3]^2}]
]

{RL,RR}={Rm[A,B,CC,R0,-1],Rm[A,B,CC,R0,1]}

ListLinePlot[{Table[{en,coherentProb[en,0,A,B,CC,R0,RL,RR][[1]]},{en,1.15,4,0.1}],Table[{en,coherentProb[en,0,A,B,CC,R0,RL,RR][[2]]},{en,1.15,4,0.1}],Table[{en,coherentProb[en,0,A,B,CC,R0,RL,RR][[3]]},{en,1.15,4,0.1}],Table[{en,incoherentProb[en,0,A,B,CC,R0]},{en,1.15,4,0.1}]},PlotStyle->{{Black,Thick},{Red,Dashed},{Blue,Thick},Green}]

(*Finally, we solve the problem numerically using NDsolve. An outgoing wave on the left with unit probability in the upper curve propagates through the crossing region, and we check how much amplitude it has on the right. 
*)

solu[ee_,a_,b_,c_,r0_,R1_,R2_]:={ψ1,ψ2}/.NDSolve[{-1/2*D[{ψ1[R],ψ2[R]},{R,2}]+(H[R,a,b,c,r0]).{ψ1[R],ψ2[R]}==ee*{ψ1[R],ψ2[R]},ψ1[R1]==1/Sqrt[Sqrt[2*ee]]Exp[I*Sqrt[2*ee]*R1],ψ1'[R1]==I*Sqrt[2*ee]*ψ1[R1],ψ2[R1]==0,ψ2'[R1]==0},{ψ1,ψ2},{R,-R1,R2}][[1]]

{f1,f2}=solu[1.1,A,B,CC,R0,0,40]

Plot[{Re[f1[R]],Re[f2[R]]},{R,0,40}]

(*we do this two kind of ad-hoc ways: by finding the amplitude and by averaging the oscillations over a small distance. Both are *)

solvo[R1_,R2_,f_]:=FindMaximum[{Abs[f[R]]^2,30<R<40},{R,35}][[1]]

solvo2[R1_,R2_,f_]:=NIntegrate[Abs[f[R]]^2,{R,35,38}]/3

vals2=Table[{ee,{f1,f2}=solu[ee,A,B,CC,R0,0,40];amp=solvo2[30,40,f1];1-amp*Sqrt[2*ee]},{ee,0.1,6,0.02}];

vals=Table[{ee,{f1,f2}=solu[ee,A,B,CC,R0,0,40];amp=solvo[30,40,f1];1-amp*Sqrt[2*ee]},{ee,0.1,6,0.02}];

(*black is the coherent probability without additional Stoke's phase. Red is the coherent probabilty with that phase. *)

(*bblue is the incoherent probability. Green and orange are the averaged and amplitude-matched numerical probabilities, respectively. *)

Show[{ListLinePlot[{Table[{en,coherentProb[en,0,A,B,CC,R0,9,11][[1]]},{en,1,6,0.1}],Table[{en,coherentProb[en,0,A,B,CC,R0,9,11][[3]]},{en,1,6,0.1}],Table[{en,incoherentProb[en,0,A,B,CC,R0]},{en,0.1,6,0.1}],vals,vals2},PlotStyle->{Black,Red,Blue,Green,Orange}]}]