Wolfram Language Paclet Repository
Community-contributed installable additions to the Wolfram Language
Time Evolution
Time Evolution
solves the equation of motion of a quantum system. Which equation that is depends on what you hand it: the Schrodinger equation for a closed system, the Lindblad master equation once jump operators are present, the Heisenberg equation when the thing being evolved is an operator, and the Kossakowski equation when the jump operators are correlated.
The call is the same shape in every case - a Hamiltonian, optionally some dissipation, an initial condition, and a time - and the result is always an object that depends on time. Supply a time to it and you get an ordinary state or operator back, which answers every property the framework knows about.
Two choices control everything else. A bare symbol for the time solves the equation exactly, through ; a interval solves it numerically, through . And an initial condition of evolves the propagator instead of any particular state, which can then be applied to whatever state you like.
{t,t0,t1}
In[1]:=
<<Wolfram`QuantumFramework`
A Driven Qubit
A Driven Qubit
A Hamiltonian drives a qubit between and . Evolving over one full period:
PauliX
|0〉
|1〉
|0〉
In[2]:=
psi=
["PauliX"],
["0"],{t,0,2Pi}
 | QuantumEvolve |
 | QuantumOperator |
| QuantumState |
The result is a state as a function of time, and is recorded on it as a parameter:
t
In[3]:=
psi["Parameters"]
Supplying a time gives the state at that time - an ordinary with nothing left symbolic:
In[4]:=
psi[1.]
From there every state property is available as usual:
In[5]:=
psi[1.]["ProbabilitiesList"]
Sampling the interval shows the population sloshing between the two levels and back:
In[6]:=
Table[Round[psi[t]["ProbabilitiesList"],0.0001],{t,0,3}]
Plotting it gives the Rabi oscillation:
In[7]:=
Plot[psi[t]["ProbabilitiesList"]〚2〛,{t,0,2Pi},AxesLabel{"t","population of |1〉"}]
What the Solution Actually Is
What the Solution Actually Is
A numerically evolved state is not a table of samples and it is not a symbolic expression. Its amplitudes are the solver's own , held as the state's array and applied to the time parameter:
In[8]:=
psi["State"]
That single object carries the whole two-component solution, so the state still knows its own shape without interpolating anything:
In[9]:=
psi["Dimensions"]
Nothing is computed until a time is supplied; is what triggers an interpolation. This matters for large systems, where the alternative - one scalar interpolation per amplitude, each carrying its own copy of the solver's time grid - costs a copy of that grid per dimension. Asking for it explicitly with shows what is being avoided:
psi[1.]
"Expand"->True
In[10]:=
expanded=
["PauliX"],
["0"],{t,0,2Pi},"Expand"True;Normal[expanded["State"]]
 | QuantumEvolve |
 | QuantumOperator |
| QuantumState |
The two forms give the same numbers:
In[11]:=
expanded[1.]["ProbabilitiesList"]
Exact Solutions
Exact Solutions
Passing a bare symbol instead of an interval solves the equation symbolically rather than numerically:
In[12]:=
exact=
["PauliX"],
["0"],t;Simplify[Normal[exact["StateVector"]]]
| QuantumEvolve |
| QuantumOperator |
| QuantumState |
That is the textbook solution, and its populations are what the plot above traced:
In[13]:=
Simplify[exact["ProbabilitiesList"],Assumptionst∈Reals]
The initial state may be symbolic too, in which case so is the answer:
In[14]:=
SimplifyNormal
["X"],
[{α,β}],t["StateVector"]
| QuantumEvolve |
 | QuantumOperator |
| QuantumState |
Omit the initial state and the time entirely and the register state is used, with the formal symbol as the time parameter. This is the shortest thing that says anything:
|0〉
t
In[15]:=
 | QuantumEvolve |
| QuantumOperator |
The numerical solution agrees with the closed form to the solver's tolerance:
In[16]:=
Max[Abs[psi[1.]["ProbabilitiesList"]-{Cos[1.]^2,Sin[1.]^2}]]
Time-Dependent Hamiltonians
Time-Dependent Hamiltonians
Nothing requires the Hamiltonian to be constant. A symbol in it that matches the evolution parameter simply makes the generator time dependent. This is a spin with splitting driven transversally at the same frequency - a resonant drive:
ω
0
In[17]:=
ω0=2Pi;ωp=Pi/10;tf=2Pi/ωp;driven=ω0
["JZ"]+2ωpCos[ω0t]
["JX"];
| QuantumOperator |
| QuantumOperator |
With no initial state given, the register state is used, so the Hamiltonian and the time interval are the whole call:
In[18]:=
rabi=
[driven,{t,0,tf}]
| QuantumEvolve |
The Bloch vector traces the slow Rabi nutation underneath the fast Larmor precession. Over the interval it makes one full circuit from the north pole and back:
2π/
ω
p
In[19]:=
Plot[Evaluate[Re[rabi[t]["BlochVector"]]],{t,0,tf},PlotRangeAll,PlotLegends{Subscript[r,x],Subscript[r,y],Subscript[r,z]}]
Halfway through, the spin has been driven to the south pole - a pulse:
π
In[20]:=
Round[Chop[Re[rabi[tf/2]["BlochVector"]]],0.0001]
The evolution is unitary throughout, so the state stays pure and the Bloch vector stays on the sphere:
In[21]:=
rabi[tf/2]["Purity"]
Observables and Trajectories
Observables and Trajectories
An evolved state answers every ordinary state property once a time is given, so an observable is measured as it would be on any state:
In[22]:=
| QuantumMeasurement |
| QuantumMeasurementOperator |
| QuantumOperator |
For a qubit the Bloch vector gives the motion directly. Under the state precesses about the axis, so the component stays zero:
PauliX
x
x
In[23]:=
Chop[psi[1.]["BlochCartesianCoordinates"]]
Drawing that as a curve on the Bloch sphere shows the precession as the great circle it is:
In[24]:=
Show
["UniformMixture"]["BlochPlot"],ParametricPlot3D[Evaluate[Re[psi[t]["BlochVector"]]],{t,0,2Pi}]
| QuantumState |
Open Systems: the Lindblad Equation
Open Systems: the Lindblad Equation
Jump operators turn the Schrodinger equation into the Lindblad master equation,
∂
t
∑
i
γ
i
L
i
†
L
i
1
2
†
L
i
L
i
with the density matrix, the jump operators and their rates. Passing them as is enough to switch equations. Here a Hamiltonian with dephasing at rate 0.6:
ρ
L
i
γ
i
Ls->rates
PauliZ
PauliX
In[25]:=
noisy=
["PauliZ"],
["PauliX"]{0.6},
["0"],{t,0,3}
| QuantumEvolve |
| QuantumOperator |
| QuantumOperator |
| QuantumState |
The state is a density matrix now rather than a vector, because dissipation does not preserve purity. Purity falls from 1 toward the maximally mixed value of 1/2:
In[26]:=
Table[noisy[t]["Purity"],{t,0,3}]
The populations equalize as it goes:
In[27]:=
noisy[3.]["ProbabilitiesList"]
Relaxing to a Thermal State
Relaxing to a Thermal State
Dephasing drives a qubit to the maximally mixed state because it has no preferred direction. Amplitude damping does have one, and at finite temperature the two directions compete. With a mean occupation of the bath, decay happens at rate and absorption at rate :
n
γ(n+1)
γn
In[28]:=
Ω=50;γ=1;n=3;H=Ω/2
["Z"];σm=
[{{0,1},{0,0}}];Ls=σm,;γs=γ{n+1,n};
| QuantumOperator |
| QuantumOperator |
†
σm
Starting from a pure state tilted off the axis:
In[29]:=
ρ0=
[{Cos[Pi/8],Exp[IPi/4]Sin[Pi/8]}]
| QuantumState |
In[30]:=
ρt=
[H,Lsγs,ρ0,{t,0,1}]
| QuantumEvolve |
By the populations have all but reached the thermal ratio , which for is :
t=1
{(n+1)/(2n+1),n/(2n+1)}
n=3
{4/7,3/7}
In[31]:=
ρt[1.]["ProbabilitiesList"]
In[32]:=
N[{n+1,n}/(2n+1)]
The trajectory spirals in from the surface of the Bloch sphere toward that point, the fast precession winding around the slow decay:
PauliZ
In[33]:=
Show
["UniformMixture"]["BlochPlot"],ParametricPlot3D[Evaluate[Re[ρt[t]["BlochVector"]]],{t,0,1}]
| QuantumState |
The rates can equally be folded into the operators as , in which case the jump operators are given as a plain list and the result is identical:
γ
i
L
i
In[34]:=
ρt2=
[H,Sqrt[γs]Ls,ρ0,{t,0,1}];Max[Abs[Normal[ρt2[1.]["DensityMatrix"]]-Normal[ρt[1.]["DensityMatrix"]]]]
| QuantumEvolve |
The Liouvillian and the Propagator
The Liouvillian and the Propagator
It acts on a density matrix flattened into a vector, so for a qubit its matrix is 4x4 rather than 2x2:
Because it is the map and not a trajectory, it can be applied to any state without solving again:
Solving the same problem directly as a Lindblad equation gives the same answer, to the solver's tolerance:
The propagator also exists for closed systems, where it is the familiar unitary. Solved exactly it is a closed form:
The Heisenberg Picture
The Heisenberg Picture
In the Heisenberg picture the observable moves and the state does not:
They agree to machine precision across the interval, since both read the same solved propagator:
The Kossakowski Equation
The Kossakowski Equation
The Lindblad equation assumes the jump operators are uncorrelated. Dropping that assumption gives the Kossakowski equation,
Two atoms decaying into a shared field is the standard example. The Hamiltonian is just their bare energies:
Each atom decays through its own lowering operator:
Evolve three initial states with one excitation between them - a product state, the singlet, and the triplet. When the Hamiltonian is already a superoperator, the initial state comes directly after it:
Count the excitations left at each time with the total-number operator:
At the end of the interval they are far apart:
Watching for Events
Watching for Events
The crossings really are crossings - evaluating the returned solution at them gives one half back:
Shading the plot between them marks out the intervals where the qubit is mostly excited:
The Equations Themselves
The Equations Themselves
Details
Details
Possible Issues
Possible Issues
Ask for amplitudes before supplying a time and the answer is a function of the parameter rather than numbers, which is rarely what is wanted from a numerical solve. Supply the time first:
If the exact solver cannot find a closed form, the output is the unevaluated solver call it gave up on, along with a message. Here a Rabi problem with a general detuning, which has no elementary solution:
Numbers in the Hamiltonian make the same problem solvable numerically, so the fallback when a closed form does not exist is to give the parameters values and pass an interval instead.