# 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 DiffEq.jl package.

Bloqade provides a special problem type `SchrodingerProblem`

that supports most of the integrator interface of `DiffEq`

, 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.SchrodingerProblem`

— Type```
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.

## 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:

`BloqadeExpr.emulate!`

— Function`emulate!(prob)`

Run emulation of a given problem.

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

```
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/|r_i-r_j|^6 n_i n_j
├─ [+] Ω(t) ⋅ ∑ σ^x_i
└─ [-] 2π ⋅ 2.0 ⋅ ∑ n_i
algorithm: CompositeAlgorithm{Tuple{Vern9, Rodas5{0, false, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}}, AutoSwitch{Vern9, Rodas5{0, false, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, Rational{Int64}, Int64}}((Vern9(true), Rodas5{0, false, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}(nothing, OrdinaryDiffEq.DEFAULT_PRECS)), AutoSwitch{Vern9, Rodas5{0, false, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, Rational{Int64}, Int64}(Vern9(true), Rodas5{0, false, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}(nothing, OrdinaryDiffEq.DEFAULT_PRECS), 10, 3, 9//10, 9//10, 2, false, 5))
options:
save_everystep: false
save_start: false
save_on: false
dense: false
```

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))[] = 100000000 ₍₂₎
```

You can use any function on the `reg`

object. For calculating observables, please see the Registers and Observables section.

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`

. Thus, do not use any function that has a `!`

in its name on the register `info.reg`

unless you are certain about what you are doing.

## Choose an ODE Solver

One of the most powerful tool of the Julia ecosystem is the DiffEq 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 many-body Schrödinger equation's stiffness is unknown, we will not be using stiff problem solvers, but instead using 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 many-body Schrödinger equation, such as `Vern6`

, `Vern7`

, and `Vern8`

. They also have relatively good memory usage when utilize GPUs.

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

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 DiffEq:Translation 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 might be large, but 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, but instead output results at the specific adaptive steps. In this situation, it's better to use fixed-step methods at the clocks 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/|r_i-r_j|^6 n_i n_j
├─ [+] 2π ⋅ 0.5 ⋅ ∑ σ^x_i
└─ [-] 2π ⋅ 2.0 ⋅ ∑ n_i
algorithm: CompositeAlgorithm{Tuple{Vern9, Rodas5{0, false, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}}, AutoSwitch{Vern9, Rodas5{0, false, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, Rational{Int64}, Int64}}((Vern9(true), Rodas5{0, false, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}(nothing, OrdinaryDiffEq.DEFAULT_PRECS)), AutoSwitch{Vern9, Rodas5{0, false, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, Rational{Int64}, Int64}(Vern9(true), Rodas5{0, false, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}(nothing, OrdinaryDiffEq.DEFAULT_PRECS), 10, 3, 9//10, 9//10, 2, false, 5))
options:
dt: 0.001
adaptive: false
save_everystep: false
save_start: false
save_on: false
dense: false
```

Here, we've specified the fixed time step as `dt = 1e-3`

. If one only expects the final state of the evolution, or the intervals between each chosen clock is 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 `KrylovEvolution`

object.

`BloqadeKrylov.KrylovEvolution`

— Type```
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/|r_i-r_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
```

## 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 it's 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`

.