Optimal Coherent Control

Quantum computers perform computations by executing circuits consisting of a set of quantum gates, and recording measurement at its output. At this level of abstraction, quantum gates are well-defined mathematical operations on a $n$-qubit Hilbert space. In practice, different hardware realizations engineer a universal set of gates through ia variety of controllable sets of physical interactions between the qubits. Optimal coherent control (OCC) provides is a framework whereby qubit control functions can be optimized to produce a desired target quantum gate.

In this tutorial example, we consider a very simple system made out of only two qubits, described by the Hamiltonian

where $\{\omega_j\}$ are the frequencies of the transmons and $g$ is the qubit exchange coupling.

Trotter simulation of quantum dynamics

Before considering the gate optimization, we first need to simulate the dynamics of the system generate by the Hamiltonian $H$. We begin by importing the relevant packages and setting up the simulation parameters:

using Random
using ITensors
using PastaQ
using Observers
using DataFrames
using Plots

GHz = 2π
MHz = 1e-3 * GHz
plot_args = (
  dpi=1000, size=(600, 300), margin=5Plots.mm, marker=:circle, markersize=2, linewidth=1
)

n = 2            # number of qubits
g = 12 * MHz     # exchange interaction
ω₁ = 5.0 * GHz    # qubit-1 frequency
ω₂ = 5.0 * GHz    # qubit-2 frequency
ω⃗ = [ω₁, ω₂]

q₁, q₂ = 1, 2         # modes ordering
modes = ["q₁", "q₂"]  # modes labels

#generate the Hilbert space
hilbert = qubits(n)

2-element Vector{ITensors.Index{Int64}}:
 (dim=2|id=830|"Qubit,Site,n=1")
 (dim=2|id=574|"Qubit,Site,n=2")

Now that we have created the Hilbert space and set the parameters, we can define the Hamiltonian. For now, we define it as a Vector of Tuple, where each Tuple represent a term in the Hamiltonian. This would be normally defined as an ITensor object (OpSum), but that is not yet fully differentiable.

function hamiltonian(ω⃗::Vector, g::Number)
  H = Tuple[]
  ω₁, ω₂ = ω⃗
  H = vcat(H, [(ω₁, "a† * a", q₁)])
  H = vcat(H, [(ω₂, "a† * a", q₂)])
  H = vcat(H, [(g, "a†b + ab†", (q₁, q₂))])
  return H
end

H = hamiltonian(ω⃗, g)

3-element Vector{Tuple}:
 (31.41592653589793, "a† * a", 1)
 (31.41592653589793, "a† * a", 2)
 (0.07539822368615504, "a†b + ab†", (1, 2))

We would like to simulate the system dynamics and record measurements, such as the average mode occupation. We use the Observers.jl package to keep track of observables. The Observer object is a container of a set of Function, which are called iteratively inside whatever iterative loop we consider. We also need to add a Function that measure the average occupation here.

function population(ψ::MPS, site::Int)
  s = siteinds(ψ)[site]
  orthogonalize!(ψ, site)
  norm2_ψ = norm(ψ)^2
  val = scalar(ψ[site] * op("a† * a", s) * dag(prime(ψ[site], s))) / norm2_ψ
  return real(val)
end;

#define a vector of observables and create the `Observer`.
observables = ["n($α)" => x -> population(x, k)  # actually x -> expect(x, "a† * a"; sites = k)
               for (k, α) in enumerate(modes)]
obs = Observer(observables)

Observers.Observer with 2 entries:
  "n(q₁)" => NamedTuple{(:f, :results), Tuple{Union{Nothing, Function}, Any}}((…
  "n(q₂)" => NamedTuple{(:f, :results), Tuple{Union{Nothing, Function}, Any}}((…

We are not ready to simulate the system dynamics using a Trotter expansion. The time-evolution propagator up to time $t$ is decomposed as

with $t_g = M\delta t$ being the finla time. Each elementary propagator is then approximated with its Trotter expansion (to order 2 by default):

tg = 30                  # final time (in ns)
trottersteps = 100       # number of Trotter steps
δt = tg / trottersteps   # step size
ts = 0.0:δt:tg           # time list

#build the Trotter circuit
circuit = trottercircuit(H; ts=ts, layered=true)

#set initial state |ψ⟩ = |1,0⟩
ψ₀ = productstate(hilbert, [1, 0])

#perform TEBD simulation and generate output `MPS`
ψ = runcircuit(
  ψ₀, circuit; (observer!)=obs, move_sites_back_before_measurements=true, outputlevel=0
)

ITensors.MPS
[1] ((dim=2|id=830|"Qubit,Site,n=1"), (dim=2|id=819|"Link,n=1"))
[2] ((dim=2|id=574|"Qubit,Site,n=2"), (dim=2|id=819|"Link,n=1"))

The measurements taken during the dynamics (one at each Trotter layer) are store in the Observer and can be retrieved into a DataFrame format (for example). We plot here the average occupation of the two modes as a function of time:

res = DataFrame(results(obs));
p = plot(; xlabel="time (ns)", ylabel="n̂(t)", legend=(0.40, 0.9), plot_args...)
p = plot!(p, ts, res[!, "n(q₁)"]; label="n(q₁)", plot_args...)
p = plot!(p, ts, res[!, "n(q₂)"]; label="n(q₂)", plot_args...)
p

In this simulation we placed the two qubits on resonance ($\omega_1=\omega_2$). By populating one of the qubit with an excitation (qubit 1 above), we observe that the dynamics swaps the excitation between the two qubits at time. In fact, this system implements a perfect iSwap gate.

In practice, in idle mode the two qubits are placed at some detuning, and placed on resonance only when realizing the gate. Here we consider this setting, and we will optimize the modulation of the two qubit frequencies to realize the gate when starting further apart. But first let's just re-run the previous dynamical simulation in this setup:

ω₁ = 5.0 * GHz
ω₂ = 5.3 * GHz
ω⃗ = [ω₁, ω₂]

H = hamiltonian(ω⃗, g)

obs = Observer(observables)

circuit = trottercircuit(H; ts=ts, layered=true)

ψ₀ = productstate(hilbert, [1, 0])

ψ = runcircuit(
  ψ₀, circuit; (observer!)=obs, move_sites_back_before_measurements=true, outputlevel=0
)

res = DataFrame(results(obs));
p = plot(; xlabel="time (ns)", ylabel="n̂(t)", legend=(0.50, 0.9), plot_args...)
p = plot!(p, ts, res[!, "n(q₁)"]; label="n(q₁)", plot_args...)
p = plot!(p, ts, res[!, "n(q₂)"]; label="n(q₂)", plot_args...)
p

Control optimization

After importing relevant packages for the optimization, we now fix the desired gate time to t_g=25`ns. We also define two control functions:

and

using Zygote
using OptimKit
using StatsBase: mean

tg = 25
trottersteps = 100
δt = tg / trottersteps
ts = 0.0:δt:tg

Λ = 20.0 * MHz
fourier_control(ϑ, t) = Λ * tanh(sum([ϑ[i] * sin(π * i * t / tg) for i in 1:length(ϑ)]))

function pulse_control(ϑ, t)
  y₀, ypulse, ton, toff, γ = ϑ
  f = tanh((t - ton) / γ) - tanh((t - toff) / γ)
  return y₀ + 0.5 * (ypulse - y₀) * f
end

pulse_control (generic function with 1 method)

We define the new Hamiltonian as follows. We take qubit 2 and send a pulse to bring its frequency near $\omega_1$. At the same time, we also introduce a frequency modulation $\omega_1(t)$ with small amplitude.

function hamiltonian(θ⃗::Vector, ω⃗::Vector, g::Number, t::Float64)
  ω₁, ω₂ = ω⃗
  ϑ₁, ϑ₂ = θ⃗
  H = Tuple[]
  H = vcat(H, [(ω₁ + fourier_control(ϑ₁, t), "a† * a", q₁)])
  H = vcat(H, [(ω₂ + pulse_control(ϑ₂, t), "a† * a", q₂)])
  H = vcat(H, [(g, "a†b + ab†", (q₁, q₂))])
  return H
end

hamiltonian(θ::Vector, t::Float64) = hamiltonian(θ, ω⃗, g, t);

Let's see how this look like:

Random.seed!(12345)
Ntones = 8
ϑ₁ = rand(Ntones)
ϑ₂ = [0.0, ω₁ - ω₂, 0.1 * tg, 0.9 * tg, 1]
θ⃗₀ = [ϑ₁, ϑ₂]
p = plot(; xlabel="time (ns)", ylabel="ωⱼ(t)", title="", legend=(0.50, 0.9), plot_args...)
p = plot!(
  p, ts, [ω₁ + fourier_control(ϑ₁, t) for t in ts] ./ GHz; label="ω₁(t)", plot_args...
)
p = plot!(
  p, ts, [ω₂ + pulse_control(ϑ₂, t) for t in ts] ./ GHz; label="ω₂(t)", plot_args...
)
p

The cost function to our optimization is computed from the inner products between the time-evolved and the desired wavefunctions:

function loss(Ψ⃗, Φ⃗, θ⃗)
  #build sequence Tuple (OpSum) Hamiltonians at different times
  Ht = [hamiltonian(θ⃗, t) for t in ts]
  #Trotter-Suzuki decomposition
  circuit = trottercircuit(Ht; ts=ts)
  #run the circuit
  UΨ⃗ = [runcircuit(ψ, circuit; cutoff=1e-7) for ψ in Ψ⃗]
  #compute inner product ⟨ϕ|ψ(t)⟩
  ip = [inner(ϕ, ψ) for (ϕ, ψ) in zip(Φ⃗, UΨ⃗)]
  return 1 - abs2(mean(ip))
end;

We not set the ideal gate (Lazy format) and define loss closure.

ideal_gate = [
  [0, 0] => (1, [0, 0]),
  [1, 0] => (-im, [0, 1]),
  [0, 1] => (-im, [1, 0]),
  [1, 1] => (1, [1, 1]),
]

Ψ⃗ = [productstate(hilbert, σ) for σ in first.(ideal_gate)]
Φ⃗ = [ϕ * productstate(hilbert, σ) for (ϕ, σ) in last.(ideal_gate)];
loss(θ) = loss(Ψ⃗, Φ⃗, θ)

loss (generic function with 2 methods)

We initialize the optimizer the run the optimization for $20$ steps:

optimizer = LBFGS(; verbosity=2, maxiter=20)
loss_n_grad(x) = (loss(x), convert(Vector, loss'(x)))
θ⃗, fs, gs, niter, normgradhistory = optimize(loss_n_grad, θ⃗₀, optimizer)
normgradhistory[:, 1]

21-element Vector{Float64}:
 0.4890014517301349
 0.15820242489327863
 0.11266945774372061
 0.04033370367053579
 0.01904319419287115
 0.018339128398784155
 0.018153952398863193
 0.01061965465033321
 0.006624699533632095
 0.00371666644128954
 ⋮
 0.001288202795409421
 0.0011482154687258994
 0.0010982123279993372
 0.0010751277753145994
 0.0010156052662232407
 0.0009813237276582454
 0.0009308468512931878
 0.0009018233244940665
 0.00089522763648342

In markdown sections we can use markdown syntax. For example, we can

Ht = [hamiltonian(θ⃗, t) for t in ts]
circuit = trottercircuit(Ht; ts=ts, layered=true)
ψ₀ = productstate(hilbert, [1, 0])
observables = ["n($α)" => x -> population(x, k) for (k, α) in enumerate(modes)]
obs = Observer(observables)
ψ = runcircuit(
  ψ₀, circuit; (observer!)=obs, move_sites_back_before_measurements=true, outputlevel=0
)
res = DataFrame(results(obs));
p = plot(; xlabel="time (ns)", ylabel="n̂(t)", legend=(0.50, 0.9), plot_args...)
p = plot!(p, ts, res[!, "n(q₁)"]; label="n(q₁)", plot_args...)
p = plot!(p, ts, res[!, "n(q₂)"]; label="n(q₂)", plot_args...)
p

This page was generated using Literate.jl.