A generative AI envisions a quantum computer based on neural atom architecture, which leverages two different species of atoms, precisely entrapped by laser technology.

Bloqade is a cutting-edge package, written in Julia, dedicated to the rapid dynamics simulation of neural atom systems, specifically the Rydberg atom system. This software was launched by QuEra Computing Inc. in the last year and has since garnered favorable reviews for its superior performance. It boasts the support of GPU acceleration, resulting in significantly faster calculations.

Benchmark of Bloqade (Julia) with popular quantum simulator QuTip (Python). Credit: Bloqade Manual

However, what if one wishes to utilize a Python-based optimizer to manage pulse control or if the user is more comfortable with Python? How can they harness the power of Bloqade? The answer is more straightforward than one might expect, thanks to the assistance of JuliaCall.

First things first: you need to make sure JuliaCall is set up and ready to go in your current Python environment. Also, you should have your project.toml handy in the current directory. All the packages you need for your project - as listed in the project.toml file - should be installed and up-to-date. And, of course, we need to double-check that Bloqade is installed.

Let’s see how to use JuliaCall:

from juliacall import Main as jl
import numpy as np
from juliacall import convert as jlconvert
import matplotlib.pyplot as plt

First, we want to activate the Julia environment in the project.toml

jl.Pkg.activate(".")

To import Julia packages, one can use jl.seval("String") function.

jl.seval("using Bloqade")
jl.seval("using Yao")

Julia supports a variety of data types. To facilitate interaction between Python and Julia, it’s necessary to appropriately convert Python data to Julia data types. This can be efficiently accomplished using juliaconvert. For example, we can transfer numpy arrays to julia arrays:

jlconvert(jl.Array, np.array([1,2,3]))
3-element Vector{Int64}:
 1
 2
 3
jlconvert(jl.Array, np.array([[1,2,3],[3.4,5.6,7.8]]))
2×3 Matrix{Float64}:
 1.0  2.0  3.0
 3.4  5.6  7.8
jlconvert(jl.Array, np.array([[1j,2,3],[3.4,5.6,7.8]]))
2×3 Matrix{ComplexF64}:
 0.0+1.0im  2.0+0.0im  3.0+0.0im
 3.4+0.0im  5.6+0.0im  7.8+0.0im
jlconvert(jl.Symbol,"a")
:a

Alright, with the necessary preparations in place, we are now equipped to delve into the use of Bloqade. As an illustrative example, we will use a scenario detailed in the manual: the adiabatic preparation of the Z2-ordered state in a one-dimensional atom chain.

nsites = 9
atoms = jl.Bloqade.generate_sites(jl.Bloqade.ChainLattice(), nsites, scale = 5.72)

We can see that we generated a chain of nine atoms:

atoms

The configuration of 9 atoms

Now let’s use a piecewise linear function to schedule the pulse of Rabi frenquency $\Omega$ and detuning $\Delta$:

total_time = 3.0
Ω_max = 2*np.pi * 4
clocks = jlconvert(jl.Array, np.array([0.0, 0.1, 2.1, 2.2, total_time]))
values = jlconvert(jl.Array, np.array([0.0, Ω_max, Ω_max, 0, 0]))
Ω = jl.Bloqade.piecewise_linear(clocks = clocks, values = values)
Ω
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀Waveform{_, Float64}⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ 
                    ┌────────────────────────────────────────┐ 
                  5 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
                    │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
                    │⠀⣀⣀⣀⢀⣀⣀⣀⣀⣀⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
                    │⠀⡏⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
                    │⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
                    │⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
                    │⢸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
   value / 2π (MHz) │⢸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
                    │⢸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
                    │⢸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
                    │⡜⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
                    │⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
                    │⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
                    │⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
                  0 │⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣇⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀│ 
                    └────────────────────────────────────────┘ 
                    ⠀0⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀clock (μs)⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀3⠀
U1 = -2*np.pi * 10
U2 = 2*np.pi * 10
clocks = jlconvert(jl.Array, np.array([0.0, 0.6, 2.1, total_time]))
values = jlconvert(jl.Array, np.array([U1, U1, U2, U2]))
Δ = jl.Bloqade.piecewise_linear(clocks = clocks, values = values)
Δ
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀Waveform{_, Float64}⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ 
                    ┌────────────────────────────────────────┐ 
                 20 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
                    │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
                    │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
                    │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀│ 
                    │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⡤⠞⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
                    │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣠⠖⠉⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
                    │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⡴⠋⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
   value / 2π (MHz) │⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⢤⡤⠾⠥⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤│ 
                    │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣠⠖⠉⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
                    │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⡴⠋⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
                    │⠀⠀⠀⠀⠀⠀⠀⠀⢀⡤⠞⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
                    │⠒⠒⠒⠒⠒⠒⠒⠒⠉⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
                    │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
                    │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
                -20 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
                    └────────────────────────────────────────┘ 
                    ⠀0⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀clock (μs)⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀3⠀

The time-dependent Hamiltonian is generated by:

h = jl.Bloqade.rydberg_h(atoms,Δ=Δ, Ω=Ω)
reg = jl.Bloqade.zero_state(9)
prob = jl.Bloqade.SchrodingerProblem(reg, total_time, h)
integrator = jl.Bloqade.init(prob, jl.Bloqade.Vern8())

Now, we can calculate some expectation values. For example, we want to get the Pauli Z expectation value of the first atom. To achieve this, we need to use some functions provided by Yao, since Bloqade is a package based on Yao, and the wavefunction register is defined by Yao:

(comment: In Yao, the operators are defined with YaoBlocks, and it is usually easy to write with Julia.Pair. However, we believe there is a bug in JuliaCall that will generate a tuple instead of Julia.Pair, so we use the following solution)

obs_Z = jl.Yao.PutBlock(9, jl.Yao.Z,(1,))
obsZ_lists = []
for _ in jl.Bloqade.TimeChoiceIterator(integrator,jl.LinRange(0,total_time,100)):
    obsZ = jl.Yao.expect(obs_Z, reg)
    obsZ_lists.append(obsZ)
plt.plot(np.linspace(0,total_time,100),obsZ_lists)

The expectation of Pauli Z of the first qubit as a function of time.

Now let’s get the Rydberg density of all the atoms during the evolution:

reg = jl.Bloqade.zero_state(9)
prob = jl.Bloqade.SchrodingerProblem(reg, total_time, h)
integrator = jl.Bloqade.init(prob, jl.Bloqade.Vern8())
densities = []
for _ in jl.Bloqade.TimeChoiceIterator(integrator, jl.LinRange(0,total_time,100)):
    densities.append(jl.Bloqade.rydberg_density(reg))
densities = np.array(densities)
fig, ax = plt.subplots(figsize = (10, 4))
shw = ax.imshow(np.transpose(densities), interpolation = "nearest", aspect = "auto", extent = [0, total_time, 0.5, nsites + 0.5])
ax.set_xlabel("time (μs)")
ax.set_ylabel("site")
ax.set_xticks(np.linspace(0,total_time,10))
ax.set_yticks(np.arange(nsites)+1)
bar = fig.colorbar(shw)

Rydberg density of 9 atoms as a function of time.

Now, we have successfully reproduced the figure of the density adiabatic evolution of 9 atoms!

Acknowledgment: We thank the authors of Bloqade, especially Dr. Chen Zhao, for their assistance in preparing this tutorial.