# Emulation of the Quantum Dynamics

After we create the Rydberg Hamiltonian and the register for storing the quantum information, we can simulate the quantum many-body dynamics. The coherent dynamics of the system is governed by the Schrödinger equation. The emulation interface of Bloqade is designed in a define-and-run style. With Bloqade, we have two major types of emulation:

• ODE-solver based emulation for most of the problems.
• Krylov based emulation for piecewise constant problems or QAOA-like problem.

### Define the ODE Emulation Problem

ODE solvers are the major backend we use for most of the exact quantum dynamics simulation. The ODE solvers for Bloqade are powered by the DifferentialEquations.jl package.

Bloqade provides a special problem type SchrodingerProblem that supports most of the integrator interface of DifferentialEquations, and most of the solver options. Here, we introduce some common use cases of the integrator and solver options. For more advanced usage of the solvers, please refer to the above link.

BloqadeODE.SchrodingerProblemType
struct SchrodingerProblem
SchrodingerProblem(reg, tspan, hamiltonian; kw...)

Define a Schrodinger equation problem that uses ODE solver from OrdinaryDiffEq to solve the dynamics.

Arguments

• register: required, the evolution problem register, can be a SubspaceArrayReg or an ArrayReg from Yao.
• tspan: required, a (start, stop) tuple or a single number t, the single value form t is equivalent to (zero(t), t).
• hamiltonian: required, the evolution hamiltonian, can be created via rydberg_h.

Common Keyword Arguments

• algo: optional, algorithm to use, this only works for the emulate! interface. for solve or integrator interface, one will need to specify the algorithm explicitly.
• progress: print progress bar or not, this may effect the performance when problem scale is small, default is true.
• progress_steps: steps to update the progress bar, default is 5.
• reltol: relative tolerance, default is 1e-8.
• abstol: absolute tolerance, default is 1e-8.

Further References

For more ODE options, please refer to Common Solver Options. The SchrodingerProblem type supports most of the standard DiffEq problem interface.

source

## Run ODE-based Emulation

To run the emulation, you need to define the exact evolution and solver you would like to run with via BloqadeODE.SchrodingerProblem, and then feed the corresponding object to the emulate! function:

For example, we can simulate the quantum dynamics of a time-dependent Hamiltonian with the following:

using Bloqade
atoms = generate_sites(SquareLattice(), 3, 3; scale=5.1);
clocks = [0.0, 0.1, 0.2, 0.3, 0.4];
wf = piecewise_constant(;clocks, values=2π*[1.0, 2.0, 3.0, 4.0]);
h = rydberg_h(atoms; Δ=2π*2.0, Ω=wf); # create the Hamiltonian
reg = zero_state(length(atoms)); # create fullspace register
ev = SchrodingerProblem(reg, 0.3, h) # the second input is the total time
emulate!(ev)
SchrodingerProblem:
register info:
type: ArrayReg{2, ComplexF64, Matrix{ComplexF64}}
storage size: 8 bytes

time span (μs): (0.0, 0.3)

equation:
storage size: 76.016 KiB
expression:
nqubits: 9
+
├─ [+] ∑ 2π ⋅ 8.627e5.0/|x_i-x_j|^6 n_i n_j
├─ [+] Ω(t) ⋅ ∑ σ^x_i
└─ [-] 2π ⋅ 2.0 ⋅ ∑ n_i

algorithm: DP8(stage_limiter! = trivial_limiter!, step_limiter! = trivial_limiter!, thread = static(false))
options:
save_everystep: false
save_start: false
save_on: false
dense: false
reltol: 1.0e-10
abstol: 1.0e-10


With emulate!, the quantum state stored in reg has been updated to the state after the time evolution.

If you want to do operations during the real-time evolution, such as measuring observables, you can instead using the integrator interface with for loop and with TimeChoiceIterator on your desired clocks, e.g.:

integrator = init(ev, Vern8())
for _ in TimeChoiceIterator(integrator, [0.1, 0.25])
ev.reg # state at selected time
@show measure(ev.reg)[] # measure the state at each time
end
(measure(ev.reg))[] = 000000000 ₍₂₎
(measure(ev.reg))[] = 101000001 ₍₂₎

You can use any function on the reg object. For calculating observables, please see the Registers and Observables section.

Tip

Remember to make sure your operation does not mutate your state so that it won't affect the evolution itself, since the entire time evolution is simulated by mutating the state vector stored in reg. Avoid using any function that has a ! in its name (indicating that it mutates its arguments per Julia convention) on the register info.reg unless you are certain about what you are doing.

## Choose an ODE Solver

One of the most powerful tools of the Julia ecosystem is the DifferentialEquations ecosystem that implements many different solvers. These solvers have different advantages and trade-offs. Since simulating a quantum many-body Schrödinger equation has some special properties compared to a general ODE problem, we will discuss some general heuristics in this section on how to choose a good ODE solver and how to check if your simulation converges. Because the stiffness of the many-body Schrödinger equation is unknown, we will opt instead to use non-stiff problem algorithms or auto-switching algorithms.

For most of the cases, one can use the VCABM solver for a large system simulation. However, this method requires more memory, which can be a bottleneck with GPUs.

The Vern family is another set of solvers that is good for the many-body Schrödinger equation, such as Vern6, Vern7, and Vern8. They also have relatively good memory usage when utilizing GPUs.

For a more detailed list of solvers, please refer to DifferentialEquations:Full List of Solvers. For more detailed explanation on ODE solvers, please refer to DifferentialEquations:Recommended Methods.

Default Solver and Auto-Switching

Bloqade uses the Dormand-Prince solver (DP8) by default for all SchrodingerProblem instances. For dealing with potentially stiff equations one can take advantage of DifferentialEquations auto-switching algorithms where a default solver can be supplied but upon encountering stiffness can switch to a specific solver. A recommended combination for Bloqade is AutoVern9(Rodas5) although other combinations are possible (see DifferentialEquations: Pre-Built Stiffness Detecting and Auto-Switching Algorithms)

If you are familiar with MATLAB or Python, you may wish to compare the same methods that you use in MATLAB or Python; you can find the corresponding solvers in Julia in DifferentialEquations:Translations from MATLAB/Python/R.

## Adaptive Steps in ODE Solvers

Our ODE solvers use adaptive steps by default. It provides a significant speedup compared to standard fixed-step methods (see our benchmark here). However, if one expects to retrieve the results during the time evolution, e.g., plotting Rydberg densities with the evolution time, fixed-step methods are sometimes preferred.

More specifically, when the adaptive steps are turned on the time steps can become large. However, if one is interested in measuring some observables in smaller time steps, then the adaptive step method will not produce accurate results for the finer time step, instead outputting results at the specific adaptive steps. In this situation, it's better to use fixed-step methods at the times where the observables are measured.

One can use the code below to turn off the adaptive steps when setting up the SchrodingerProblem:

atoms = generate_sites(SquareLattice(), 3, 3; scale=5.1);
h = rydberg_h(atoms; Δ=2π*2.0, Ω= 2π*1.0); # create the Hamiltonian
reg = zero_state(length(atoms));
prob = SchrodingerProblem(reg, 3.0, h, adaptive = false, dt=1e-3);
SchrodingerProblem:
register info:
type: ArrayReg{2, ComplexF64, Matrix{ComplexF64}}
storage size: 8 bytes

time span (μs): (0.0, 3.0)

equation:
storage size: 83.992 KiB
expression:
nqubits: 9
+
├─ [+] ∑ 2π ⋅ 8.627e5.0/|x_i-x_j|^6 n_i n_j
├─ [+] 2π ⋅ 0.5 ⋅ ∑ σ^x_i
└─ [-] 2π ⋅ 2.0 ⋅ ∑ n_i

algorithm: DP8(stage_limiter! = trivial_limiter!, step_limiter! = trivial_limiter!, thread = static(false))
options:
dt: 0.001
save_everystep: false
save_start: false
save_on: false
dense: false
reltol: 1.0e-10
abstol: 1.0e-10


Here, we've specified the fixed time step as dt = 1e-3. If one only wants the final state of the evolution or if the intervals between each chosen time are much larger than the maximum step size, then adaptive steps are preferred.

### Define the Krylov Emulation Problem

The Krylov-based method expects time-independent Hamiltonians. One can define such a time evolution via a KrylovEvolution object.

BloqadeKrylov.KrylovEvolutionType
struct KrylovEvolution
KrylovEvolution(reg::AbstractRegister, clocks, h; kw...)

Create a KrylovEvolution object that describes a time evolution using Krylov subspace methods.

Arguments

• reg: a register, should be a subtype of AbstractRegister.
• clocks: the clocks of this time evolution at each step.
• h: a hamiltonian expression.

Keyword Arguments

• progress: show progress bar, default is false.
• progress_name: progress bar name, default is "emulating".
• normalize_step: normalize the state every normalize_step.
• normalize_finally: wether normalize the state in the end of evolution, default is true.
• tol: tolerance of the Krylov method, default is 1e-7

Examples

The following is the simplest way of using KrylovEvolution via emulate!. For more advanced usage, please refer to documentation page Emulation.

julia> using Bloqade

julia> r = zero_state(5)
ArrayReg{2, ComplexF64, Array...}
active qudits: 5/5
nlevel: 2

julia> atoms = [(i, ) for i in 1:5]
5-element Vector{Tuple{Int64}}:
(1,)
(2,)
(3,)
(4,)
(5,)

julia> h = rydberg_h(atoms; Ω=sin)
nqubits: 5
+
├─ [+] ∑ 5.42e6/|x_i-x_j|^6 n_i n_j
└─ [+] Ω(t) ⋅ ∑ σ^x_i

julia> prob = KrylovEvolution(r, 0.0:1e-2:0.1, h);

julia> emulate!(prob); # run the emulation
source

## Run Krylov-based Emulation

We can run the Krylov-based emulation in a similar way using emulate!:

julia> emulate!(KrylovEvolution(reg, clocks, h))KrylovEvolution:
register info:
type: ArrayReg{2, ComplexF64, Matrix{ComplexF64}}
storage size: 8 bytes

clocks
start:0.0μs
last:0.4μs

Hamiltonian
number of dynamic terms: 0
storage size: 40 bytes

Options:
progress=false
progress_step=1
progress_name="emulating"
normalize_step=5
normalize_finally=true
tol=1.0e-7

However, as its name suggests, the Krylov-based emulation is not a standard ODE problem that DiffEq supports. Thus, it does not support the ODE problem interface, but is more like a gate-based interface. For example, the object KrylovEvolution is iterable:

for (step, reg, duration) in KrylovEvolution(reg, clocks, h)
@show step
@show reg
@show duration
println("==========")
end
step = 1
reg = ArrayReg{2, ComplexF64, Array...}
active qubits: 9/9
nlevel: 2
duration = 0.0
==========
step = 2
reg = ArrayReg{2, ComplexF64, Array...}
active qubits: 9/9
nlevel: 2
duration = 0.1
==========
step = 3
reg = ArrayReg{2, ComplexF64, Array...}
active qubits: 9/9
nlevel: 2
duration = 0.2
==========
step = 4
reg = ArrayReg{2, ComplexF64, Array...}
active qubits: 9/9
nlevel: 2
duration = 0.3
==========
step = 5
reg = ArrayReg{2, ComplexF64, Array...}
active qubits: 9/9
nlevel: 2
duration = 0.4
==========

## Krylov vs ODE Solvers

The KrylovEvolution uses the Krylov subspace methods to simulate the time evolution of time-independent operators $\exp(i\Delta t_i H)$, where $\Delta t_i$ is the duration of time-independent Hamiltonian $H$ at time $t$. This method is more efficient when the evolution itself is a discrete evolution, e.g. in QAOA and with piecewise_constant waveforms. As for other cases, ODE solvers are usually more efficient than KrylovEvolution.