Edit on GitHub

CZ gate in Rydberg atoms

In order to implement digital quantum computation in Rydberg atoms, ground state $|g\rangle$ and hyperfine state $|h\rangle$ are used as qubit states. Transitions between these two levels $|g\rangle \leftrightarrow |h\rangle$ can be induced with Raman pulses. However, in order to create entanglement between two atoms, another ancillary atomic level has to be introduced and it's called the Rydberg level $|r\rangle$ (essentially a highly excited state). Transitions $|h\rangle \leftrightarrow |r\rangle$ are induced with so-called Rydberg pulses, and Rydberg interaction, which prohibits two neighboring atoms to both occupy Rydberg levels, can be used to generate entanglement between the atoms.

Here we set $|g\rangle \equiv |1\rangle$, $|h\rangle \equiv |2\rangle$, and $|r\rangle \equiv |3\rangle$. Hamiltonian of two Rydberg atoms can then be written as

\[H = \frac{\Omega(t)}{2}\left[e^{i\varphi(t)} \left(\sigma_+^{(1)} + \sigma_+^{(2)}\right) + e^{-i\varphi(t)} \left(\sigma_-^{(1)} + \sigma_-^{(2)}\right)\right] - \Delta(t) (n_1 + n_2) + V n_1n_2,\]

where $\sigma^{(i)}_+ = |3\rangle\langle2|_i$, $\sigma^{(i)}_- = |2\rangle\langle3|_i$, $n_i = |3\rangle\langle3|_i$ and $V$ is the Rydberg interaction strength.

The usual technique to realize a $CZ$ gate (see Levine et al.) consists of shining two Rydberg pulses with a phase jump in between them. The required phase jump, pulse duration and detuning can be analytically calculated from the Rabi frequency of the Rydberg pulse.

In this notebook, we find a single Rydberg pulse profile to do a $CZ$ gate.

using Sisyphus
using QuantumOptics
using Flux, DiffEqFlux
using Plots
using Random
using ProgressMeter

ProgressMeter.ijulia_behavior(:clear)

Since pulses are required to be real-valued, we have to rewrite the Hamiltonian

\[H = Vn_1n_2 - \Delta(t)(n_1 + n_2) + \frac{\Omega(t)\cos\varphi(t)}{2}\left(\sigma^{(1)}_x + \sigma^{(2)}_x\right) + \frac{\Omega(t)\sin\varphi(t)}{2}\left(\sigma^{(1)}_y + \sigma^{(2)}_y\right).\]

V = 2π*10.0 # MHz

bs = NLevelBasis(3)
bsys = bs⊗bs

id = identityoperator(bs)
sp1 = transition(bs, 3, 2)⊗id
sm1 = transition(bs, 2, 3)⊗id
sp2 = id⊗transition(bs, 3, 2)
sm2 = id⊗transition(bs, 2, 3)
n1 = transition(bs, 3, 3)⊗id
n2 = id⊗transition(bs, 3, 3)

H0 = V*n1*n2
H1 = (n1 + n2)
H2 = (sp1 + sp2 + sm1 + sm2)
H3 = 1.0im*(sp1 + sp2 - sm1 - sm2)
states= [nlevelstate(bs, 1)⊗nlevelstate(bs, 1),
         nlevelstate(bs, 1)⊗nlevelstate(bs, 2),
         nlevelstate(bs, 2)⊗nlevelstate(bs, 1),
         nlevelstate(bs, 2)⊗nlevelstate(bs, 2)]

trans = UnitaryTransform(states, [[1.0 0 0 0 ];[0 1.0 0 0 ];[0 0 1.0 0 ]; [0 0 0 -1.0]])
t0, t1 = 0.0, 1.0 # in μs
Ω(p, t) = piecewise_const_interp(p[begin:length(p)÷3], t; t0=t0, t1=t1)
Δ(p, t) = piecewise_const_interp(p[length(p)÷3+1:2*length(p)÷3], t; t0=t0, t1=t1)
φ(p, t) = piecewise_const_interp(p[2*length(p)÷3+1:end], t; t0=t0, t1=t1)

θ = vcat(2π * ones(10), 2π * ones(10), zeros(6), π * ones(6))
println("Number of parameters: ", length(θ))
Number of parameters: 32
green = theme_palette(:auto).colors.colors[3]

ts = t0:t1/1000:t1
plot(ts, [Ω(θ, t) for t in ts], label="Ω/2π")
plot!(ts, [Δ(θ, t) for t in ts], label="Δ/2π")
ylabel!("Waveform (MHz)")
plot!(Plots.twinx(), ts, [φ(θ, t) for t in ts], ylabel="ϕ (Radians)", color=green,
      yguidefont = font(green), legend=false)
xlabel!("Time (µs)", right_margin = 15Plots.mm)

svg

We use (x, y) -> 1.0 - real(x'*y) as a cost function since it takes into account proper phases of final states.

coeffs(p, t) = [-Δ(p, t), Ω(p, t) * cos(φ(p, t)) / 2, Ω(p, t) * sin(φ(p, t)) / 2]
cost = CostFunction((x, y) -> 1.0 - real(x'*y))
H = Hamiltonian(H0, [H1, H2, H3], coeffs)

prob = QOCProblem(H, trans, (t0, t1), cost)
sol = solve(prob, θ, Adam(0.1); maxiter=150, abstol=1e-6, reltol=1e-6)
Progress: 100%|█████████████████████████████████████████| Time: 0:06:18
  distance:     8.055867831952002e-5
  constraints:  0.0
plot(ts, [Ω(sol.params, t) for t in ts], label="Ω/2π")
plot!(ts, [Δ(sol.params, t) for t in ts], label="Δ/2π")
ylabel!("Waveform (MHz)")
plot!(Plots.twinx(), ts, [φ(sol.params, t) for t in ts], ylabel="ϕ (Radians)", color=green,
      yguidefont = font(green), legend=false)
xlabel!("Time (µs)", right_margin = 15Plots.mm)

svg

tout, psit22 = schroedinger_dynamic(ts, nlevelstate(bs, 2)⊗nlevelstate(bs, 2), H, sol.params)
tout, psit21 = schroedinger_dynamic(ts, nlevelstate(bs, 2)⊗nlevelstate(bs, 1), H, sol.params)
tout, psit12 = schroedinger_dynamic(ts, nlevelstate(bs, 1)⊗nlevelstate(bs, 2), H, sol.params)
tout, psit11 = schroedinger_dynamic(ts, nlevelstate(bs, 1)⊗nlevelstate(bs, 1), H, sol.params)

plot(tout, [real((nlevelstate(bs, 2)⊗nlevelstate(bs, 2))'*elm) for elm in psit22], label="|hh⟩")
plot!(tout, [real((nlevelstate(bs, 2)⊗nlevelstate(bs, 1))'*elm) for elm in psit21], label="|hg⟩")
plot!(tout, [real((nlevelstate(bs, 1)⊗nlevelstate(bs, 2))'*elm) for elm in psit12], label="|gh⟩")
plot!(tout, [real((nlevelstate(bs, 1)⊗nlevelstate(bs, 1))'*elm) for elm in psit11], label="|gg⟩")
plot!(xlabel="Time (µs)", ylabel="Overlap (⟨ij|ψ⟩)", legend=:right)

svg