Time Evolution Solutions
Solution
QuantumToolbox
utilizes the powerful DifferentialEquation.jl
to simulate different kinds of quantum system dynamics. Thus, we will first look at the data structure used for returning the solution (sol
) from DifferentialEquation.jl
. The solution stores all the crucial data needed for analyzing and plotting the results of a simulation. A generic structure TimeEvolutionSol
contains the following properties for storing simulation data:
Fields (Attributes) | Description |
---|---|
sol.times | The list of time points at which the expectation values are calculated during the evolution. |
sol.times_states | The list of time points at which the states are stored during the evolution. |
sol.states | The list of result states corresponding to each time point in sol.times_states . |
sol.expect | The expectation values corresponding to each time point in sol.times . |
sol.alg | The algorithm which is used during the solving process. |
sol.abstol | The absolute tolerance which is used during the solving process. |
sol.reltol | The relative tolerance which is used during the solving process. |
sol.retcode (or sol.converged ) | The returned status from the solver. |
Accessing data in solutions
To understand how to access the data in solution, we will use an example as a guide, although we do not worry about the simulation details at this stage. The Schrödinger equation solver (sesolve
) used in this example returns TimeEvolutionSol
:
H = 0.5 * sigmay()
ψ0 = basis(2, 0)
e_ops = [
proj(basis(2, 0)),
proj(basis(2, 1)),
basis(2, 0) * basis(2, 1)'
]
tlist = LinRange(0, 10, 100)
sol = sesolve(H, ψ0, tlist, e_ops = e_ops, progress_bar = Val(false))
To see what is contained inside the solution, we can use the print
function:
print(sol)
Solution of time evolution
(return code: Success)
--------------------------
num_states = 1
num_expect = 3
ODE alg.: OrdinaryDiffEqTsit5.Tsit5{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}(OrdinaryDiffEqCore.trivial_limiter!, OrdinaryDiffEqCore.trivial_limiter!, static(false))
abstol = 1.0e-8
reltol = 1.0e-6
It tells us the number of expectation values are computed and the number of states are stored. Now we have all the information needed to analyze the simulation results. To access the data for the three expectation values, one can do:
expt1 = real(sol.expect[1,:])
expt2 = real(sol.expect[2,:])
expt3 = real(sol.expect[3,:])
Recall that Julia
uses Fortran
-style indexing that begins with one (i.e., [1,:]
represents the 1-st observable, where :
represents all values corresponding to tlist
).
Together with the list of time points at which these expectation values are calculated:
times = sol.times
we can plot the resulting expectation values:
# plot by CairoMakie.jl
fig = Figure(size = (500, 350))
ax = Axis(fig[1, 1], xlabel = L"t")
lines!(ax, times, expt1, label = L"\langle 0 | \rho(t) | 0 \rangle")
lines!(ax, times, expt2, label = L"\langle 1 | \rho(t) | 1 \rangle")
lines!(ax, times, expt3, label = L"\langle 0 | \rho(t) | 1 \rangle")
ylims!(ax, (-0.5, 1.0))
axislegend(ax, position = :lb)
fig
State vectors, or density matrices, are accessed in a similar manner:
sol.states
1-element Vector{QuantumObject{Ket, Dimensions{1, Tuple{Space}}, Vector{ComplexF64}}}:
Quantum Object: type=Ket() dims=[2] size=(2,)
2-element Vector{ComplexF64}:
0.2836622138778657 + 0.0im
-0.9589242324580376 + 0.0im
Together with the list of time points at which these states are stored:
times = sol.times_states
Here, the solution contains only one (final) state. Because the states
will be saved depend on the keyword argument saveat
in kwargs
. If e_ops
is empty, the default value of saveat=tlist
(saving the states corresponding to tlist
), otherwise, saveat=[tlist[end]]
(only save the final state). One can also specify e_ops
and saveat
separately.
Some other solvers can have other output.
Multiple trajectories solution
The solutions are different for solvers which compute multiple trajectories, such as the TimeEvolutionMCSol
(Monte Carlo) or the TimeEvolutionStochasticSol
(stochastic methods). The storage of expectation values and states depends on the keyword argument keep_runs_results
, which determines whether the results of all trajectories are stored in the solution.
When the keyword argument keep_runs_results
is passed as Val(false)
to a multi-trajectory solver, the states
and expect
fields store only the average results (averaged over all trajectories). The results can be accessed by the following index-order:
sol.states[time_idx]
sol.expect[e_op,time_idx]
For example:
tlist = LinRange(0, 1, 11)
c_ops = (destroy(2),)
e_ops = (num(2),)
sol_mc1 = mcsolve(H, ψ0, tlist, c_ops, e_ops=e_ops, ntraj=25, keep_runs_results=Val(false), progress_bar=Val(false))
size(sol_mc1.expect)
(1, 11)
If the keyword argument keep_runs_results = Val(true)
, the results for each trajectory and each time are stored, and the index-order of the elements in fields states
and expect
are:
sol.states[trajectory,time_idx]
sol.expect[e_op,trajectory,time_idx]
For example:
sol_mc2 = mcsolve(H, ψ0, tlist, c_ops, e_ops=e_ops, ntraj=25, keep_runs_results=Val(true), progress_bar=Val(false))
size(sol_mc2.expect)
(1, 25, 11)
We also provide the following functions for statistical analysis of multi-trajectory sol
utions:
Functions | Description |
---|---|
average_states(sol) | Return the trajectory-averaged result states (as density Operator ) at each time point. |
average_expect(sol) | Return the trajectory-averaged expectation values at each time point. |
std_expect(sol) | Return the trajectory-wise standard deviation of the expectation values at each time point. |