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.

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

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:

emulate!(prob)

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

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

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.

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
    adaptive: false
    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.

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

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

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.