Let’s define a two-level system by the time-independent Hamiltonian
H
g
=
ϵ
2
σ
z
+
Δ
2
σ
x
where ϵ represents the bias (tilt) between wells, and Δ is the tunnel splitting. A classical drive in the lab frame is given by
H
drive
=
1
2
V(t)Cos(
ω
d
t+ϕ)
σ
x
with
ω
d
the carrier (drive) frequency, and ϕ the carrier phase.​The Rabi drive can be written as
Ω(t)=
1
2
V(t)
ϵ
2
Δ
+
2
ϵ
. For a π-pulse, we set it as
∫Ω(t)t=π
; thus for a sequence of square pulses with duration τ, one finds
V(t)=
2π
τ
ϵ
2
Δ
+
2
ϵ
∑
k
Pulse
t-
t
k
τ
with the pulse
Pulse(x)=

1
0<=x<1
0
else
. The overall Hamiltonian is given by
H=
H
g
+
H
drive
Define the pulse:
In[]:=
ClearAll[pulse,t];​​pulse[τ_,t_,tkList_]:=Sum
1.
0<=t-tk<τ
0
True
,{tk,tkList}
Visualize a series of pulses:
In[]:=
Plot[pulse[.3,t,Range[5]-.3],{t,0,6},Exclusions->None,AspectRatio->1/2]
Out[]=
The above behavior assumes coherent tunneling, i.e., no noise. In practice, environmental interactions introduce dephasing, relaxation, and other channels, so the dynamics deviate from the Schrödinger equation and the qubit becomes an open system. A common model for such open-system dynamics is the Lindblad master equation:
∂
t
ρ=-[H,ρ]+
∑
k
[
L
k
][ρ]
with
[
L
k
][ρ]=
L
k
ρ
†
L
k
-
1
2

†
L
k
L
k
,ρ
with
L
k
the corresponding Lindblad/jump operators.
Define the Lindblad terms:
In[]:=
ClearAll[,ℒ]​​[L_,ρ_]:=L.ρ.ConjugateTranspose[L]-
1
2
(ConjugateTranspose[L].L.ρ+ρ.ConjugateTranspose[L].L)​​ℒ[ρ_,h_,Lk_List]:=-(h.ρ-ρ.h)+Sum[[L,ρ],{L,Lk}]
Assign values to parameters and set DEs and solve it
In[]:=
(*Timeoffirstpulse*)t0=1;​​(*τdurationofapulse*)τ=.1;​​(*Gapbetweentwopulses*)gap=1;​​(*timesatwhichapulseisapplied*)tk=Table[t0+τ/2+k,{k,0,gapFloor[T/gap],gap}];​​(*
ω
d
/
ω
p
with
ω
p
=
2
ϵ
+
2
Δ
*)ratio=1;​​(*ϕcarrierphase*)ϕ=0;​​(*Δtunnelsplitting*)Δ=2π5;​​(*ϵbias(tilt)betweenwells*)ϵ=2π10;​​(*γrelaxationrate*)γ=.75;​​(*
Γ
ϕ
dephasingrate*)Γϕ=1;​​(*overallevolutiontime*)T=4;​​(*carrierdrive*)ωd=ratio
2
ϵ
+
2
Δ
;​​​​(*Hamiltonian*)h=
ϵ
2
PauliMatrix[3]+
Δ
2
PauliMatrix[2]+
1
2
2π
τ
ϵ
2
Δ
+
2
ϵ
pulse[τ,t,tk]Cos[ωdt+ϕ]PauliMatrix[1];​​(*Lindbladoperatorswithrates*)Lk=
γ
{{0,0},{0,1}},
Γϕ
2
PauliMatrix[3];​​(*initialdenistymatric*)ρ0=With{θ=0,ϕ=0},
2
Cos
θ
2

,
1
2
-ϕ

Sin[θ],
1
2
ϕ

Sin[θ],
2
Sin
θ
2

;​​(*SolutionofLindbladmasterequation*)ρt=NDSolveValue[{
′
ρ
[t]==ℒ[ρ[t],h,Lk],ρ[0]==ρ0},ρ,{t,0,T}];​​(*Blochvector*)vector=Re@Table[Tr[ρt[t].PauliMatrix[j]],{j,3}];​​​​Manipulate​​Row@​​Column"H=
ϵ
2
σ
z
+
Δ
2
σ
x
+
1
2
V(t)Cos(
ω
d
t+ϕ)
σ
x
","
∂
t
ρ=-[H,ρ]+γ​
σ
-
ρσ
+
-
1
2
​{​
σ
+
σ
-
​,ρ}​+
Γ
ϕ
2
​(​
σ
z
ρσ
z
-ρ)","withtherelaxationrateγ=1
T
1
andthepuredephasingrate
Γ
ϕ
=1
T
ϕ
.",​​​​Plot[{vector[[1]],vector[[2]],vector[[3]],pulse[τ,t,tk]},{t,0,tf},PlotRange{-1,1},PlotLegendsPlaced[{"〈
σ
x
〉","〈
σ
y
〉","〈
σ
z
〉","π-pulses"},Below],ImageSizeMedium,AxesLabel{"Time",""},GridLinesAutomatic,Exclusions->None,PlotPoints->100],​​​​Show
Graphics3D[
]
,ParametricPlot3D[vector,{t,0,tf}],ImageSize->Medium,{tf,.1T,T},SaveDefinitions:>True
Out[]=
​
tf
H=
ϵ
2
σ
z
+
Δ
2
σ
x
+
1
2
V(t)Cos(
ω
d
t+ϕ)
σ
x
∂
t
ρ=-[H,ρ]+γ​
σ
-
ρσ
+
-
1
2
​{​
σ
+
σ
-
​,ρ}​+
Γ
ϕ
2
​(​
σ
z
ρσ
z
-ρ)
withtherelaxationrateγ=1/
T
1
andthepuredephasingrate
Γ
ϕ
=1
T
ϕ
.
〈
σ
x
〉
〈
σ
y
〉
〈
σ
z
〉
π-pulses
​​