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.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 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.

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. 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.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/|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
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 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.