Interactively Simulating Open Quantum Systems

Dr Harini Hapuarachchi, RMIT University

Join the Conversation #AustralasiaWTC2024

Open Quantum Systems - A Quick Overview

Open quantum systems are real-world, imperfect quantum systems irreversibly interacting with their environments. 



Further details:

Interactive Simulation using ParametricNDSolve & Manipulate

(***Function to output the Hamiltonian matrix, where δ = (Subscript[(Subscript[ω, 0])-ω, d]) ***)
getH[δ_,Ω0_, ϕ0_]:=-{{(δ/2), Ω0*Cos[ϕ0]-I*Ω0*Sin[ϕ0]},{Ω0*Cos[ϕ0]+I*Ω0*Sin[ϕ0], -(δ/2)}};


(***Functions to obtain the RHS of the QME***)
CT[M_]:=ConjugateTranspose[M];
getRhoDot[δ_,Ω0_, ϕ0_,γd_,ρ_]:=Module[{H, decayTerm}(* local variables of Module[] *),
H = getH[δ,Ω0, ϕ0];                          
Lk={{0,1},{0,0}}; 
decayTerm = γd*(Lk.ρ.CT[Lk]-1/2*((CT[Lk].Lk).ρ+ρ.(CT[Lk].Lk)));                          
rhoDot = -I*(H.ρ-ρ.H)+decayTerm]                                                       


(*** Define static parameters and the initial density matrix ***)
δ=0.;                                        (*Detuning δ = (Subscript[(Subscript[ω, 0])-ω, d])*)
ϕ0=0.;                                     (*Phase of the laser field*)
ρ0={{1,0},{0,0}};      (*Initial density matrix*)


(********** Preparing the inputs for ParametricNDsolve **********)
ρ[t_]:={{ρ11[t],ρ12[t]},{ρ21[t],ρ22[t]}};   (* Assemble rho using PNDSolve's variables*)

(* Build matrix equations and convert to scalar equation sets*)
scalarEqns = Thread[Flatten/@(D[ρ[t],t]==getRhoDot[δ,Ω0,ϕ0,γd,ρ[t]])];
initConditions = Thread[Flatten/@(ρ[0]==ρ0)];
(**************************************************************)


(***Solve the differential equations using ParametricNDSolve***)
sol=ParametricNDSolve[        (* sol holds the solutions as parametric interpolation functions*)
{scalarEqns,initConditions},(*Differential equations and initial conditions*)
{ρ11,ρ12,ρ21,ρ22},                   (* Dependent variables *)
{t,0.,T},                                          (* Independent variable and its range *)
{Ω0,γd,T}];                                     (* Parameters *)


(************* Alternative method for time evolution **************)
(* Instead of ParametricNDSolve, QuantumEvolve within the Wolfram Quantum Framework 
could also be used for the time evolution, given the Liouvillian superoperator.*)
(***************************************************************)


(*** Creating dynamic plots of the solutions***)
Manipulate[  (*** The code within Manipulate runs every time a dynamic param. changes ***)

(*** Plotting the solutions of ParametricNDsolve ***)
Plot[Evaluate[{ρ11[Ω0, γd, T][t], ρ22[Ω0, γd, T][t]}/.sol],{t,0,T},
PlotLegends->Placed[{"(Subscript[ρ, gg])","(Subscript[ρ, ee])"}, Right],Frame->True,FrameLabel->{"Time","Populations"}, 
LabelStyle->Directive[Black,15.0],PlotRange->All, ImageSize->{300,180}],

(*** Optional text in Manipulate Console ***)
Style["Two-Level System (TLS) Simulator\n",Bold,18],

(*** Configuring the dynamic parameter sliders ***)
(*** Format: {{parameter,defaultValue,displayText,},sliderMin,sliderMax,sliderIncrement,...}*)
{{Ω0, 1., "Laser coupling ((Subscript[Ω, 0])):"}, 0.1, 2.0, 0.1, Appearance-> "Labeled", LabelStyle->{15}},  
{{γd,0., "Decay rate ((Subscript[γ, d])):"}, 0., 2.0, 0.1, Appearance->"Labeled", LabelStyle->{15}},
{{T,5π, "Evolution time:"}, 0., 10π, π/10, Appearance->"Labeled", LabelStyle->{15}},

(***Manipulate console properties***)
AppearanceElements->"ResetButton",
ControlPlacement->Left, 
LabelStyle->{15.0}
]

Fitting Parameters for Experiments

(* IMPORTANT: LETTER 'S' HAS BEEN APPENDED TO PARAMETER NAMES THAT COULD CONFLICT WITH OTHER SIMULATIONS IN THE SLIDE SHOW*)

(*Clear["Global`*"]*)(* Clear all definitions and values associated with symbols in the Global context *)

(*Import and paste below*) (*Ref: https://mathematica.stackexchange.com/questions/204745/how-can-i-import-a-csv-file-into-mathematica*)
(*SetDirectory[NotebookDirectory[]]
Import["ExpDemoData.csv"]*)

ExpDemoData={{0.075265`,0.032678`},{0.41129`,0.084539`},{0.5324`,0.16686`},{0.71379`,0.3067`},{0.69571`,0.46911`},{1.0156`,0.56721`},{0.95265`,0.66872`},{1.2106`,0.81194`},{1.5773`,0.85703`},{1.9611`,0.80287`},{1.9776`,0.73295`},{1.9942`,0.65851`},{2.3167`,0.6021`},{2.5021`,0.50623`},{2.3041`,0.44533`},{2.4279`,0.36976`},{2.7502`,0.32463`},{2.8275`,0.28403`},{2.8127`,0.25245`},{3.0886`,0.23664`},{3.3031`,0.23438`},{3.394`,0.28851`},{3.7605`,0.35052`},{3.7593`,0.41931`},{3.9111`,0.50051`},{4.1249`,0.5332`},{4.2766`,0.61891`},{4.4442`,0.67078`},{4.5362`,0.66175`},{4.8422`,0.67865`},{4.9652`,0.65045`},{5.3182`,0.60531`},{5.3347`,0.53426`},{5.5343`,0.50718`},{5.5195`,0.47448`},{5.8727`,0.42032`},{6.0879`,0.3673`},{6.3333`,0.34924`},{6.7613`,0.3932`},{6.7606`,0.43268`},{7.0205`,0.46199`},{7.3558`,0.55445`},{7.7226`,0.60067`},{8.0751`,0.5826`},{8.5661`,0.52506`},{8.7509`,0.46414`},{9.2264`,0.42126`},{9.7775`,0.43025`},{10.205`,0.49113`},{10.648`,0.54524`},{11.276`,0.54069`},{11.522`,0.50459`},{12.243`,0.4538`},{12.962`,0.46278`},{13.635`,0.51011`},{14.416`,0.52134`},{14.707`,0.4999`},{15.443`,0.46941`},{16.101`,0.48854`},{16.729`,0.49753`},{17.295`,0.51554`},{18.169`,0.47715`},{18.674`,0.49178`},{19.93`,0.48833`}};
                                                

(* Static parameters and initial density matrix *)
δS=0.; ϕ0S=0.;  
ρ0S={{1,0},{0,0}};

(* Function to assemble the density matrix using symbolic matrix elements *)
ρS[tS_]:={{ρ11S[tS],ρ12S[tS]},{ρ21S[tS],ρ22S[tS]}};

(* We then compute matrix differential equations and compare element-wise *)
scalarEqnsS = Thread[Flatten/@(D[ρS[tS],tS]==getRhoDot[δS,Ω0S,ϕ0S,γdS,ρS[tS]])];
initConditionsS = Thread[Flatten/@(ρS[0]==ρ0S)];

(*Solve the differential equations [1]]*)
solS=ParametricNDSolve[
{scalarEqnsS,initConditionsS},
{ρ11S,ρ12S,ρ21S,ρ22S},(* Dependent variables *)
{tS,0.,TS},  (* Independent variable and range *)
{Ω0S,γdS,TS}];(* Parameters *)

(* Initialise local dynamic variables within the scope of DynamicModule to avoid conflicts with other parts of the notebook*)
DynamicModule[{TS=5 π,Ω0S=1.,γdS=0.}, 
Manipulate[ (* Code within Manipulate[] runs every a dynamic parameter changes*)

(* Plotting the ParametricNDsolve solutions *)
PopnPlot = Plot[Evaluate[{ρ22S[Ω0S, γdS, TS][tS]}/.solS],{tS,0,TS},PlotStyle->{Red,Thick},
Frame->True,FrameLabel->{"Time","Counts ∝ (Subscript[ρ, ee])"}, PlotLegends->Placed[{"Simulation"}, Bottom],LabelStyle->Directive[Black,15.0],
PlotRange->{{0,15},{0,1}}, ImageSize->{300,180}];
ExpPlot = ListPlot[ExpDemoData, PlotStyle->{Black,PointSize[0.03]}, PlotLegends->Placed[{"Exp."}, Bottom], LabelStyle->Directive[Black,15.0]];
Show[PopnPlot, ExpPlot],

(* Optional text in Manipulate Console *)
Style["Two-Level System (TLS) Simulator\n",Bold,18],

(* Dynamic parameters *)
{{Ω0S, 1., "Laser coupling ((Subscript[Ω, 0])):"}, 0.1, 2.0, 0.05, Appearance-> "Labeled", LabelStyle->{15}},
{{γdS,0., "Decay rate ((Subscript[γ, d])):"}, 0., 0.5, 0.1, Appearance->"Labeled", LabelStyle->{15}},
{{TS,5π, "Evolution time:"}, 0., 10π, π/10, Appearance->"Labeled", LabelStyle->{15}},

(* Manipulate console properties *)
AppearanceElements->"ResetButton",
ControlPlacement->Left, 
LabelStyle->{15.0}
]]

Visualising Qubit Dynamics

(* IMPORTANT: LETTER 'S' HAS BEEN APPENDED TO PARAMETER NAMES THAT COULD CONFLICT WITH OTHER SIMULATIONS IN THE SLIDE SHOW*)

 (*Clear["Global`*"]*) (* Clear all definitions and values associated with symbols in the Global context *)

(* ****** Preparing the functions needed ******* *)

(* Functions to get Hamiltonian and Lindblad terms in angular frequency units *)
getHB[δB_,Ω0B_, ϕ0B_]:=-(δB/2)*PauliMatrix[3]-Ω0B*Cos[ϕ0B]*PauliMatrix[1]-Ω0B*Sin[ϕ0B]*PauliMatrix[2];
getLindTermB[γB_,σB_,ρB_]:=γB*(σB.ρB.ConjugateTranspose[σB]-1/2*((ConjugateTranspose[σB].σB).ρB+ρB.(ConjugateTranspose[σB].σB)));

(* Function to output the time derivative of density matrix ρ, as per the Lindblad master equation *)
getRhoDotB[δB_,Ω0B_, ϕ0B_,ρB_,γ1B_(*, γ2_*)]:=Module[{HB, decayTermB, dephTermB},(* These variables are treated as local by Module *)
HB = getHB[δB,Ω0B, ϕ0B];                                                                                                          (* Hamiltonian in angular frequency units *)
σ1B={{0,1},{0,0}}; decayTermB = getLindTermB[γ1B,σ1B,ρB];                              (* Decay term of the Lindblad master equation *)
(*σ2=-PauliMatrix[3]; dephTerm = getLindTermB[γ2,σ2,ρ];     *)                           (* Dephasing term of the Lindblad master equation *)
rhoDotB = -I*(HB.ρB-ρB.HB)+decayTermB(*+dephTerm*)]                                                        (* Output of the function *)

(* Function to compute Bloch vector given a density matrix*)
blochVectorFromρ[ρB_]:={Tr[ρB.PauliMatrix[1]],Tr[ρB.PauliMatrix[2]],Tr[ρB.PauliMatrix[3]]};

(* Function to compute the density matrix given a Bloch vector*)
ρFromBlochVector[{rx_, ry_, rz_}]:=(IdentityMatrix[2]+rx*PauliMatrix[1]+ry*PauliMatrix[2]+rz*PauliMatrix[3])/2;
BlochVectorFromInitPosStr[InitPosStr_]:=
If[InitPosStr=="|0⟩", r={0,0,1} ,
If[InitPosStr=="|1⟩", r={0,0,-1}, 
If[InitPosStr=="|+⟩", r={1,0,0},
If[InitPosStr=="|-⟩", r={-1,0,0},
If[InitPosStr=="|i⟩", r={0,1,0},
If[InitPosStr=="|-i⟩", r={0,-1,0},Nothing]]]]]];




(* ************************* Solving the Lindblad Master Equation with ParametricNDSolve ************************* *)

(* Function to compute the density matrix using symbolic matrix elements (which will be the independent variables of ParametricNDSolve) *)
ρB[tB_]:={{ρ11B[tB],ρ12B[tB]},{ρ21B[tB],ρ22B[tB]}};

(*Solve the differential equations [1]]*)
solB=ParametricNDSolve[{
(* Coupled differential equations to solve are obtained by element-wise comparison (Flattening) the matrix equation [2] *)
scalarEqnsB = Thread[Flatten/@(D[ρB[tB],tB]==getRhoDotB[δB,Ω0B,ϕ0B, ρB[tB],γ1B(*,γ2B*)])],       
ρ11B[0] == ρ11initB,ρ12B[0]==ρ12initB,ρ21B[0]==ρ21initB,ρ22B[0]==ρ22initB},(* Initial conditions (note these are set to parameters) *)
{ρ11B,ρ12B,ρ21B,ρ22B},                                                                                                                       (* Dependent variables we solve for *)
{tB,0.,TB},                                                                                                                                              (* Independent variable and its range *)
{δB, Ω0B,ϕ0B, γ1B, (*γ2B, *)TB, ρ11initB,ρ12initB,ρ21initB,ρ22initB}];                              (* Parameters the solution will depend on *)

(* If ParametricNDSolve solutions should be used in more than one place, we define solution variables as follows [3] *)
{ρ11SolB, ρ12SolB, ρ21SolB, ρ22SolB}={ρ11B,ρ12B,ρ21B,ρ22B}/.solB;




(* *********** Drawing the sphere and axes of the Bloch sphere to be displayed within Manipulate ***********  *)
blochSphere=SphericalPlot3D[1,{θ,0,π},{φ,0,2 π},Mesh->{20,20},MeshStyle->Directive[Blue,Thin,Opacity[0.2]],PlotStyle->Directive[Opacity[0.1],LightBlue],Lighting->{"Neutral",Blue},Axes->False];

axes=Graphics3D[{
Black,Arrow[{{0,0,0},{1.8,0,0}}],Text[Style["(OverHat[x])",14],{1.9,0,0}],
Black,Arrow[{{0,0,0},{0,1.8,0}}],Text[Style["(OverHat[y])",14],{0,1.9,0}],
Black,Arrow[{{0,0,0},{0,0,1.8}}],Text[Style["(OverHat[z])",14],{0,0.2,1.7}],
Dashed,Line[{{0,0,0},{-1.0,0,0}}],Dashed,Line[{{0,0,0},{0,-1.0,0}}],Dashed,Line[{{0,0,0},{0,0,-1.0}}],

Text[Style["Bloch vector evolution", Bold, 15],{0.,0.,2.0}],
Red,PointSize[Large],Point[{0.1,0.,-1.7}], 
Blue,PointSize[Medium],Point[{0.1,0.,-2.0}],
Black, Text[Style["\nInitial \nFinal",13],{-.6,0,-2.0}]


}];


DynamicModule[{InitPosStrB= "|0⟩", TB= π, δB=0.,Ω0B=1., ϕ0B=0.,γ1B=0.(*, γ2B=0.*)},  (* DynamicModule ensures that dynamic variables reside in a local scope *)
Manipulate[(* The code within Manipulate[] below is run every time the user changes a dynamic parameter *)

(* Create initial density matrix ρ0 using user input and define a 3D point to display later *)
InitialBlochVector = BlochVectorFromInitPosStr[InitPosStrB];
ρ0B=ρFromBlochVector[InitialBlochVector];
StartingPos = Graphics3D[{ Red,PointSize[Large],Point[InitialBlochVector], 
If[InitPosStrB =="|0⟩",Text[Style[InitPosStrB,12,Bold,Red],InitialBlochVector+{-0.1,0.12,0.2}],
If[InitPosStrB =="|1⟩",Text[Style[InitPosStrB,12,Bold,Red],InitialBlochVector+{ 0.0,0.0,-0.3}],
If[InitPosStrB =="|+⟩",Text[Style[InitPosStrB,12,Bold,Red],InitialBlochVector+{ 0.2,0.0,-0.1}], 
If[InitPosStrB =="|-⟩",Text[Style[InitPosStrB,12,Bold,Red],InitialBlochVector+{ -0.9,0.0,0.}],
If[InitPosStrB =="|i⟩",Text[Style[InitPosStrB,12,Bold,Red],InitialBlochVector+{ -0.5,0.05,0.}], 
If[InitPosStrB =="|-i⟩",Text[Style[InitPosStrB,12,Bold,Red],InitialBlochVector+{ -0.5,-1.0,0.}],Nothing]]]]]]}];

(* Obtain the density matrix solution at user defined parameters as a function of t, and define a 3D plot of Bloch vector from t=0 to T *)
rhoSolB = {
{ρ11SolB[δB, Ω0B, ϕ0B, γ1B, (*γ2B,*) TB, ρ0B[[1]][[1]],ρ0B[[1]][[2]],ρ0B[[2]][[1]],ρ0B[[2]][[2]]][tB], ρ12SolB[δB,Ω0B, ϕ0B, γ1B, (*γ2B,*) TB, ρ0B[[1]][[1]],ρ0B[[1]][[2]],ρ0B[[2]][[1]],ρ0B[[2]][[2]]][tB]},
{ρ21SolB[δB, Ω0B, ϕ0B, γ1B, (*γ2B,*) TB, ρ0B[[1]][[1]],ρ0B[[1]][[2]],ρ0B[[2]][[1]],ρ0B[[2]][[2]]][tB], ρ22SolB[δB, Ω0B, ϕ0B, γ1B, (*γ2B,*) TB, ρ0B[[1]][[1]],ρ0B[[1]][[2]],ρ0B[[2]][[1]],ρ0B[[2]][[2]]][tB]}};
BlochEvolution = ParametricPlot3D[blochVectorFromρ[rhoSolB], {tB,0,TB}, PlotStyle->Black];

(* Obtain the final density matrix solution and its Blochvector (at T). Define a 3D Blochvector point to display later. *)
ρFinal={
{ρ11SolB[δB, Ω0B, ϕ0B, γ1B, (*γ2B,*) TB, ρ0B[[1]][[1]],ρ0B[[1]][[2]],ρ0B[[2]][[1]],ρ0B[[2]][[2]]][TB],ρ12SolB[δB, Ω0B, ϕ0B, γ1B, (*γ2B,*) TB, ρ0B[[1]][[1]],ρ0B[[1]][[2]],ρ0B[[2]][[1]],ρ0B[[2]][[2]]][TB]},
{ρ21SolB[δB, Ω0B, ϕ0B, γ1B, (*γ2B,*) TB, ρ0B[[1]][[1]],ρ0B[[1]][[2]],ρ0B[[2]][[1]],ρ0B[[2]][[2]]][TB],ρ22SolB[δB, Ω0B, ϕ0B, γ1B, (*γ2B,*) TB, ρ0B[[1]][[1]],ρ0B[[1]][[2]],ρ0B[[2]][[1]],ρ0B[[2]][[2]]][TB]}};
FinalBlochVector = blochVectorFromρ[ρFinal];
EndingPos = Graphics3D[{Blue,PointSize[Medium],Point[FinalBlochVector]}];


(* We display plots within the following Grid *)
Grid[{{

(* Show the 3D plots prepared earlier, using the same view parameters *)
Show[blochSphere,BlochEvolution, axes,StartingPos, EndingPos,
PlotRange->{{-2,2},{-2,2},{-2,2}},
Axes->False,Boxed->False,Background->White, 
ViewVertical->{0,0,1},ViewPoint->{1.0,1.0,.5},
ImageSize->{400,400}]
(*
(* Plotting the density matrix elements in an inner Grid *)
Grid[{{
Plot[
{ρ11SolB[δB, Ω0B, ϕ0B, γ1B, (*γ2B,*) TB, ρ0B[[1]][[1]],ρ0B[[1]][[2]],ρ0B[[2]][[1]],ρ0B[[2]][[2]]][tB],
ρ22SolB[δB, Ω0B, ϕ0B, γ1B, (*γ2B,*) TB, ρ0B[[1]][[1]],ρ0B[[1]][[2]],ρ0B[[2]][[1]],ρ0B[[2]][[2]]][tB]},
{tB,0,TB},PlotLegends->Placed[{"(Subscript[ρ, gg])    ","(Subscript[ρ, ee])    "}, Right],Frame->True,FrameLabel->{"","Populations"}, LabelStyle->Directive[Black,11],PlotRange->{{0,TB},{0., 1.0}}]},

{Plot[
{Re[ρ12SolB[δB, Ω0B, ϕ0B, γ1B, (*γ2B,*) TB, ρ0B[[1]][[1]],ρ0B[[1]][[2]],ρ0B[[2]][[1]],ρ0B[[2]][[2]]][tB]],
Im[ρ12SolB[δB, Ω0B, ϕ0B, γ1B, (*γ2B,*) TB, ρ0B[[1]][[1]],ρ0B[[1]][[2]],ρ0B[[2]][[1]],ρ0B[[2]][[2]]][tB]]},
{tB,0,TB},PlotLegends->Placed[{"Re[(Subscript[ρ, ge])]","Im[(Subscript[ρ, ge])]"},Right],Frame->True,FrameLabel->{"Time","Coherences"}, LabelStyle->Directive[Black,11],PlotRange->{{0,TB},{-0.5, 0.5}}]
}}]
*)

}}],

(*Style["TLS Bloch Trajectory Simulator",Bold, 17], *)
Style["Bloch vector (Overscript[v, →])(t)=Tr[ρ(t)(Overscript[σ, →])]", Bold, 17],
Style["(where (Overscript[σ, →]) = X(OverHat[x]) + Y(OverHat[y]) + Z(OverHat[z]))", 13],

(* User defined dynamic parameters *)
Style["\nControl parameters:",15, Bold],
{{InitPosStrB, "|0⟩", "Initial State |(Subscript[ψ, 0])⟩:"}, {"|0⟩", "|1⟩", "|+⟩", "|-⟩","|i⟩", "|-i⟩"}},
(*Style["(Where (Subscript[ω, 0]) is the natural frequency of TLS)",11],*)
{{δB,0., "Detuning (δ=(Subscript[ω, 0])-(Subscript[ω, d])):"}, -1, 1, 0.1, Appearance-> "Labeled", LabelStyle->{13}},
{{Ω0B,1., "Laser coupling ((Subscript[Ω, 0])):"}, -1, 1, 0.1, Appearance-> "Labeled", LabelStyle->{13}},
{{ϕ0B,0., "Laser phase ((Subscript[ϕ, 0])):"}, 0., 2π, 0.05π, Appearance-> "Labeled", LabelStyle->{13}},
{{TB,2π, "Evolution time:"}, 0.01,10,0.1,  Appearance->"Labeled", LabelStyle->{13}}, 

Style["\nEnvironmental influence(s):",15, Bold],
{{γ1B,0., "Decay rate:"}, 0., 1., 0.1, Appearance->"Labeled", LabelStyle->{13}},
(*{{γ2B,0., "Deph. rate:"}, 0., 1.,0.1, Appearance->"Labeled", LabelStyle->{13}},*)

(*Style["\nInitial pure state (where ρ[t=0]=|(Subscript[ψ, 0])⟩⟨(Subscript[ψ, 0])|) & time:", 11, Bold],*)
(*Style["\nInitial pure state & time:", 15, Bold],*)
  
AppearanceElements->"ResetButton",
ControlPlacement->Left ,
LabelStyle->{13.0}
]]