Tutorial
This introductory tutorial illustrates the basic usage of the GRAPE package for an example of a quantum gate on two qubits. It tries to be a self-contained as possible, to highlight some of the core concepts and features of the GRAPE package:
- Using arbitrary data structures for quantum states and operators
- Defining a control problem based on multiple "trajectories", with multi-threading
- Optimizing multiple control fields at the same time (the real and imaginary part of a physical control)
- Using non-trivial "control amplitudes" to ensure smooth switch-on/off
- Using automatic differentiation to minimize arbitrary optimization functionals
- Applying bounds on the amplitude of the control field
A Two-Qubit System
The quantum state of a two-qubit system is described by a complex vector spanned by the possible classical bit configurations, $|00⟩$, $|01⟩$, $|10⟩$, and $|11⟩$ in braket notation. This basis is represented as
using StaticArrays: SVector
basis = [
SVector{4}(ComplexF64[1, 0, 0, 0]), # |00⟩
SVector{4}(ComplexF64[0, 1, 0, 0]), # |01⟩
SVector{4}(ComplexF64[0, 0, 1, 0]), # |10⟩
SVector{4}(ComplexF64[0, 0, 0, 1]), # |11⟩
];An arbitrary vector $|Ψ⟩ = α_{00} |00⟩ + α_{01} |01⟩ + α_{10} |10⟩ + α_{11} |11⟩$ with complex coefficients $α_i = α_{00}$, $α_{01}$, $α_{10}$, $α_{11}$ is understood according to the Born rule to specify the probability as $|α_{i}|^2$ to find the system in the corresponding possible classical state on measurement.
We've used the StaticArrays Julia package to encode the quantum states as an SVector, which provides numerical advantages for very small vectors such as these, but also illustrates an important design choice of the GRAPE package: states can be expressed in any data structure fulfilling a well-specified interface. This allows for custom, problem-specific encodings.
Prior to a measurement, quantum mechanics postulates that a quantum state $|Ψ(t)⟩$ evolves according to the Schrödinger equation, based on a Hamiltonian matrix $\hat{H}$. We consider here a setup inspired by superconducting transmon qubits. Each qubit is driven by an oscillating microwave field with frequency near the resonance for that qubit, and the two qubits are coupled by a shared transmission line. As is common in simulating quantum systems with a fast-oscillating field, the rotating wave approximation allows us to formulate the system in the "rotating frame" of the microwave field center frequency $ω_{mw}$. The Hamiltonian for the two-qubit system is then
\[\hat{H} = \left(δ_1 \hat{n}_1 + \frac{Ω(t)}{2} \hat{a}_1 + \frac{Ω^*(t)}{2} \hat{a}_1^\dagger\right) + \left(δ_2 \hat{n}_2 + \frac{Ω(t)}{2} \hat{a}_2 + \frac{Ω^*(t)}{2} \hat{a}_2^\dagger\right) + J (\hat{a}_1\hat{a}_2^\dagger + \hat{a}_1^\dagger\hat{a}_2) \tag{1}\]
with the usual raising and lowering operators $\hat{a}^\dagger$ and $\hat{a}$ and the number operator $\hat{n} = \hat{a}^\dagger \hat{a}$, and where $δ_{1,2}$ is the detuning of microwave angular central frequency from the transition frequency of the first and second qubit, respectively, $J$ is the static coupling, and $Ω(t)$ is the envelope of the microwave field. In the rotating frame, when the physical microwave can deviate from its central frequency $ω_{mw}$, this is equivalent to a complex amplitude, where the complex phase of $Ω(t)$ is the phase relative to $ω_{mw} t$. Note that Eq. (1) is an oversimplified model for actual transmons, which should include more levels than just the two lowest-lying levels defining the qubit.
Units
Generally, the elements of the Hamiltonian are in units of energy. However, the QuantumPropagators package that GRAPE uses to simulate all dynamics assumes $ħ = 1$, turning all elements of the Hamiltonian into angular frequencies (from the Planck relation $E = ħω$), expressed in units of 2π⋅Hz. It also makes "energy" and "time" directly reciprocal (only the phase $ωt$ is relevant for the dynamics). In numerics generally, it is best for all quantities to have magnitudes between maybe 10⁻³ and 10³ to avoid floating point errors. Here, the numerical quantities are the elements of the operators and the values in the time grid. With superconducting qubits typically having energies in a GHz regime (with detunings in MHz) and a timescale of ns for operations, we can define corresponding "internal" units:
const GHz = 1.0;
const MHz = 0.001GHz;
const ns = 1.0;We can then specify any of the "energy" parameters below with units, e.g.,
using LinearAlgebra: ⋅ # just so that the code looks nicer
δ₁ = 100⋅2π⋅MHz;Control Amplitudes
$Ω(t)$ in Eq. (1) is our control amplitude; the aim of optimal control is to find the particular $Ω(t)$ that steers the system in some particular way. In our case: to implement a quantum gate. However, the implementation of optimal control in GRAPE assumes real-valued controls. Thus, we have to split the complex $Ω(t)$ into an independent real and imaginary part.
We may also want to place some physical constraints onto $Ω(t)$. For example, we may want $Ω(t)$ to smoothly switch on from zero and off to zero at the beginning and end of the time grid. One way of achieving this is to define $Ω(t) = S(t) ϵ(t)$ where $S(t)$ is a static shape that has the smooth switch-on and off, and $ϵ(t)$ is an arbitrary function that we can optimize freely.
The QuantumPropagators framework provides a ShapedAmplitude object to implement this:
using QuantumPropagators.Amplitudes: ShapedAmplitudeFor a fixed time grid ending at
T = 400ns;we can define the function $S(t)$ as
using QuantumPropagators.Shapes: flattop
shape(t) = flattop(t, T = T, t_rise = 15ns);where flattop is a function that switches on smoothly from zero, reaches 1.0 after 15 ns, and remains constant at 1.0 until the last 15 ns before the final time $T =$ 400 ns, where it smoothly switches off to zero.
We can then define an initial guess for $ϵ(t)$ as a constant function and assemble that into the complete $Ω(t)$:
ϵ_re_guess(t; Ω₀ = 35⋅2π⋅MHz) = Ω₀
ϵ_im_guess(t) = 0.0
Ω_re_guess = ShapedAmplitude(ϵ_re_guess; shape);
Ω_im_guess = ShapedAmplitude(ϵ_im_guess; shape);By setting the guess for the imaginary part to zero, we let the system start exactly at the central frequency of the microwave field.
In general, the choice of the initial guess function can have a significant impact on the optimization. What makes a good guess pulse, in terms of which shape to use, or what initial amplitude (35⋅2π MHz, here) is often the result of some physical intuition, or some trial and error.
With the above definition, Ω_re_guess and Ω_im_guess are "control amplitudes" that are the physical control, whereas ϵ_re_guess and ϵ_im_guess are the "controls" from the perspective of GRAPE. The basic task of the GRAPE method is to discretize these controls to the intervals of the time grid ("pulses", and then iteratively update the pulse values at each time interval, based on the gradient of some optimization functional with respect to the pulse value. This distinction between physical control amplitudes and optimization control functions / pulses is a core design aspect of GRAPE, providing great flexibility in the models that can be used for optimization.
We can now define an explicit time grid
tlist = collect(range(0, T, step = 0.1ns));and plot the combined $Ω(t)$:
using Plots
using QuantumPropagators.Controls: discretize
"""
Plot the given complex pulse in two panels: absolute value and phase.
```julia
fig = plot_complex_pulse(tlist, Ω; time_unit=:ns, ampl_unit=:(2π⋅MHz), kwargs...)
```
generates a plot of the complex field `Ω` over the time grid `tlist`.
Arguments:
* `tlist`: A vector of time grid values
* `Ω`: A complex vector of the same length as `tlist` or a function `Ω(t)` returning a complex number
* `time_unit`: A symbol that evaluates to the conversion factor for the time unit
* `ampl_unit`: A symbol that evaluates to the conversion factor of the amplitude unit
All other keyword arguments are forwarded to `Plots.plot`.
```
"""
function plot_complex_pulse(tlist, Ω; time_unit=:ns, ampl_unit=:(2π⋅MHz), kwargs...)
Ω = discretize(Ω, tlist) # make sure Ω is defined on *points* of `tlist`
s_ampl_unit = string(ampl_unit)
if startswith(s_ampl_unit, "(2π) ⋅ ")
s_ampl_unit = "2π " * s_ampl_unit[11:end]
end
ax1 = plot(
tlist ./ eval(time_unit),
abs.(Ω) ./ eval(ampl_unit);
label="|Ω|",
xlabel="time ($time_unit)",
ylabel="amplitude ($s_ampl_unit)",
kwargs...
)
ax2 = plot(
tlist ./ eval(time_unit),
angle.(Ω) ./ π;
label="ϕ(Ω)",
xlabel="time ($time_unit)",
ylabel="phase (π)"
)
plot(ax1, ax2, layout=(2, 1))
end
const 𝕚 = 1im
fig = plot_complex_pulse(tlist, t -> Ω_re_guess(t) + 𝕚 * Ω_im_guess(t))
We see the expected shape with the smooth switch-on/-off and the amplitude of 35⋅2π MHz that we specified, while from the perspective of the GRAPE method, the "guess control" is a constant function.
Hamiltonian
With the separation of $Ω(t)$ into real and imaginary parts, we can now express the Hamiltonian in Eq. (1) numerically:
import QuantumPropagators: hamiltonian
using StaticArrays: SMatrix
"""Construct the time-dependent Hamiltonian for a two-qubit transmon system.
```julia
Ĥ = transmon_hamiltonian(; δ₁, δ₂, J, Ω_re, Ω_im)
```
constructs a [Generator](@extref `QuantumPropagators.Generators.Generator`) for
a two-transmon system truncated to two levels for each qubit, in the rotating
wave approximation (RWA)
Arguments:
* `δ₁`: The detuning of qubit 1 from the RWA angular frequency
* `δ₂`: The detuning of qubit 2 from the RWA angular frequency
* `J`: The static coupling between the two qubits, as angular frequency
* `Ω_re`: The amplitude (angular frequency) of the real part of the microwave drive
* `Ω_im`: The amplitude (angular frequency) of the imaginary part of the microwave drive
Both `Ω_re` and `Ω_im` can be given as a function `Ω_re(t)` and `Ω_im(t)`, or as a
vector of values on the time grid, or on the intervals of the time grid.
"""
function transmon_hamiltonian(; δ₁, δ₂, J, Ω_re, Ω_im)
Ĥ₀ = SMatrix{4,4,ComplexF64}(
[0 0 0 0
0 δ₂ J 0
0 J δ₁ 0
0 0 0 δ₁+δ₂]
)
Ĥ₁_re = (1/2) * SMatrix{4,4,ComplexF64}(
[0 1 1 0
1 0 0 1
1 0 0 1
0 1 1 0]
)
Ĥ₁_im = (𝕚/2) * SMatrix{4,4,ComplexF64}(
[0 -1 -1 0
1 0 0 -1
1 0 0 -1
0 1 1 0]
)
return hamiltonian(Ĥ₀, (Ĥ₁_re, Ω_re), (Ĥ₁_im, Ω_im))
endwhere we have used the QuantumPropagators.Generators.hamiltonian function to construct an object that serves as a time-dependent Hamiltonian (a "Generator", in the terminology of the general QuantumControl framework, ensuring adherence to the required interface). The SMatrix used for the drift and control Hamiltonian operators match the SVector used to encode states.
We use some arbitrary but reasonable values here, with the central frequency of the microwave field centered between the frequencies of the two qubits, resulting in a detuning of ±100⋅2π MHz and a static coupling of 3⋅2π MHz.
Ĥ = transmon_hamiltonian(;
δ₁ = 100⋅2π⋅MHz,
δ₂ = -100⋅2π⋅MHz,
J = 3⋅2π⋅MHz,
Ω_re = Ω_re_guess,
Ω_im = Ω_im_guess,
)Generator with 3 ops and 2 amplitudes
ops::Vector{StaticArraysCore.SMatrix{4, 4, ComplexF64, 16}}:
ComplexF64[0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im -0.6283185307179586 + 0.0im 0.01884955592153876 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.01884955592153876 + 0.0im 0.6283185307179586 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im]
ComplexF64[0.0 + 0.0im 0.5 + 0.0im 0.5 + 0.0im 0.0 + 0.0im; 0.5 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.5 + 0.0im; 0.5 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.5 + 0.0im; 0.0 + 0.0im 0.5 + 0.0im 0.5 + 0.0im 0.0 + 0.0im]
ComplexF64[0.0 + 0.0im -0.0 - 0.5im -0.0 - 0.5im 0.0 + 0.0im; 0.0 + 0.5im 0.0 + 0.0im 0.0 + 0.0im -0.0 - 0.5im; 0.0 + 0.5im 0.0 + 0.0im 0.0 + 0.0im -0.0 - 0.5im; 0.0 + 0.0im 0.0 + 0.5im 0.0 + 0.5im 0.0 + 0.0im]
amplitudes::Vector{QuantumPropagators.Amplitudes.ShapedContinuousAmplitude}:
ShapedAmplitude(::typeof(Main.ϵ_re_guess); shape::typeof(Main.shape))
ShapedAmplitude(::typeof(Main.ϵ_im_guess); shape::typeof(Main.shape))
Optimization Target
We will now use the GRAPE method to find an $Ω(t)$ that implements a CNOT, one of the fundamental building blocks of a quantum computer. In classical terms, a CNOT gate flips the state of the second qubit if and only if the first qubit is "1". This is fundamentally an entangling operation, and the equivalent of an "if" statement in a quantum circuit.
In terms of the basis states, we have corresponding target states:
target_states = [
SVector{4}(ComplexF64[1, 0, 0, 0]), # |00⟩
SVector{4}(ComplexF64[0, 1, 0, 0]), # |01⟩
SVector{4}(ComplexF64[0, 0, 0, 1]), # |10⟩ → |11⟩
SVector{4}(ComplexF64[0, 0, 1, 0]), # |11⟩ → |10⟩
];Beyond a classical interpretation, any quantum state $|Ψ⟩$ should be mapped to the state
\[|Ψ^{(tgt)}⟩ = \hat{O} |Ψ⟩ = \sum_k α_k \hat{O} |k⟩\]
for the four basis states $|k⟩ = |00⟩$, $|01⟩$, $|10⟩$, $|11⟩$, with
\[\hat{O} ≡ \begin{pmatrix} 1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & 0 & 1\\ 0 & 0 & 1 & 0 \end{pmatrix}\,.\]
being the matrix representation of the CNOT gate.
We can use the QuantumPropagators propagators to simulate ("propagate") the $|10⟩$ state exemplarily. QuantumPropagators can use a variety of methods to do this. We use the Cheby method here, which is extremely efficient at solving the Schrödinger equation for piecewise-constant pulses by expanding the time evolution operator into Chebychev polynomials.
using QuantumPropagators: propagate, Cheby
pops = propagate(
basis[3], Ĥ, tlist; method=Cheby, storage=true, observables=(Ψ -> Array(abs2.(Ψ)), )
)
fig = plot(tlist ./ ns, pops'; label=["00" "01" "10" "11"], xlabel="time (ns)", legend=:outertop, legend_column=-1)
From the dynamics, we can see that the $|10⟩$ initial state remains approximately unchanged by the guess controls, instead of evolving into the $|11⟩$ state as it should.
In order to find a field $\Omega(t)$ that implements the desired CNOT gate, the GRAPE method minimizes an appropriate (user-supplied) mathematical functional. We define here
@doc raw"""Square-Modulus functional
```julia
J_T = J_T_sm(Ψ, trajectories)
```
calculates the real-valued scalar
```math
J_{T,sm} = 1 - \frac{1}{N^2} \left\vert\sum_n ⟨Ψ_k(T)|Ψ_k^{(tgt)}⟩\right\vert^2
```
where the state ``|Ψ_k⟩`` is the k'th element of `Ψ` and ``|Ψ_k^{(tgt))⟩`` is
the state stored in the `target_state` attribute of the k'th element of the
`trajectories` list. The ``|Ψ_k⟩`` should generally be the states obtained from
propagating the state stored in the `initial_state` attribute of the k'th
element of `trajectories` forward to some final time ``T``.
Conceptually, this becomes zero if and only if all the forward propagated
states have an overlap of 1 with their respective target state, up to a
global phase.
"""
function J_T_sm(Ψ, trajectories)
N = length(trajectories)
f = zero(ComplexF64)
for (Ψ_k, traj) in zip(Ψ, trajectories)
f += Ψ_k ⋅ traj.target_state
end
return 1.0 - (abs2(f) / N^2)
endto implement
\[J_{T,sm} = 1 - \frac{1}{N^2} \left\vert\sum_n ⟨Ψ_k(T)|Ψ_k^{(tgt)}⟩\right\vert^2\]
appropriate for optimization of a quantum gate if the states $|Ψ_k(T)⟩$ are evolved from the $N = 4$ basis states (for a two-qubit gate) and the target states $|Ψ_k^{(tgt)}⟩$ are chosen as $\hat{O} |Ψ_k(T)⟩$. A generalization of the functional J_T (which is one of the standard functionals of quantum control) is also implemented in QuantumControl.Functionals.J_T_sm. We have re-implemented it here by hand to illustrate the API that the GRAPE package expects for a functional J_T. In particular, GRAPE requires optimization functionals to be defined in terms of trajectories, as the way to implicitly define the final-time states $|Ψ_k(T)⟩$ that enter the functional:
using GRAPE: Trajectory
trajectories = [
Trajectory(Ψ, Ĥ; target_state=Ψ_tgt)
for (Ψ, Ψ_tgt) in zip(basis, target_states)
]4-element Vector{Trajectory{StaticArraysCore.SVector{4, ComplexF64}, QuantumPropagators.Generators.Generator{StaticArraysCore.SMatrix{4, 4, ComplexF64, 16}, QuantumPropagators.Amplitudes.ShapedContinuousAmplitude}}}:
Trajectory with 4-element StaticArraysCore.SVector{4, ComplexF64} with indices SOneTo(4) initial state, Generator with 3 ops and 2 amplitudes, 4-element StaticArraysCore.SVector{4, ComplexF64} with indices SOneTo(4) target state
Trajectory with 4-element StaticArraysCore.SVector{4, ComplexF64} with indices SOneTo(4) initial state, Generator with 3 ops and 2 amplitudes, 4-element StaticArraysCore.SVector{4, ComplexF64} with indices SOneTo(4) target state
Trajectory with 4-element StaticArraysCore.SVector{4, ComplexF64} with indices SOneTo(4) initial state, Generator with 3 ops and 2 amplitudes, 4-element StaticArraysCore.SVector{4, ComplexF64} with indices SOneTo(4) target state
Trajectory with 4-element StaticArraysCore.SVector{4, ComplexF64} with indices SOneTo(4) initial state, Generator with 3 ops and 2 amplitudes, 4-element StaticArraysCore.SVector{4, ComplexF64} with indices SOneTo(4) target stateEach trajectory, at a minimum, defines an initial state Ψ and a dynamical generator Ĥ (i.e., a time-dependent Hamiltonian, in our case). In general, these can be arbitrary problem-specific objects, as long as they implement the required interfaces. Each trajectory can also attach arbitrary additional attributes that may be used by the J_T function, e.g., a target_state in our case.
Defining an optimization in terms of multiple "trajectories" has a number of benefits. For one, simulating the dynamics for the different trajectories can be performed in parallel, resulting in a potential speedup proportional to the number of trajectories. Second, it also enables advanced use cases such as "ensemble optimizations" that consider multiple "copies" of the quantum system, each with a different noisy Hamiltonian, in an effort to find controls that are robust with respect to these variations.
The important point is that all of the trajectories share the same set of controls, two in our case (the real and imaginary part of $Ω(t)$ divided by our shape $S(t)$):
using QuantumPropagators.Controls: get_controls
get_controls(trajectories)(Main.ϵ_re_guess, Main.ϵ_im_guess)The GRAPE method works by numerically evaluating gradients of the given functional J_T with respect to the pulse values at each interval of the time grid. As derived fully in Background, the resulting scheme depends explicitly on J_T via the definition of a set of "adjoint states",
\[|χ_k(T)⟩ \equiv - \frac{\partial J_T}{\partial \bra{Ψ_k(T)}}\,,\]
The calculation of these states must be implemented in a function chi that will be passed to the optimization alongside J_T. For $J_{T,sm}$, the states $|χ_k(T)⟩$ can be calculated analytically, as
\[|χ_k(T)⟩ = \frac{1}{N^2} \left(\sum_j ⟨Ψ_j^{(tgt)}|Ψ_j(T)⟩\right) |Ψ_k^{(tgt)}⟩\,,\]
and is implemented in QuantumControl.Functionals.chi_sm.
However, for more advanced functionals, J_T may not always have an analytic derivative that can be written in closed form, or the derivative is just too cumbersome to write out. The GRAPE package allows constructing chi via automatic differentiation, e.g., via the popular Zygote package. This feature can be enabled as follows:
import Zygote
using GRAPE
GRAPE.set_default_ad_framework(Zygote)[ Info: QuantumControl: Setting Zygote as the default provider for automatic differentiation.Optimization
We can now run the actual optimization by calling GRAPE.optimize. The main positional arguments are the trajectories that enter the functional and the time grid tlist. Beyond that, the required key arguments are J_T to give the functional and prop_method to specify the method to be used for any internal time propagation. We do not specify chi, since we are relying on the automatic differentiation that we activated above.
The optimization will run for as many iterations as given by iter_stop (5000 by default), and we can give a check_convergence callback to stop the optimization earlier, based on some value of the functional being reached. There is also a callback to print some information after each iteration.
By setting use_threads=true, we parallelize the optimization over the different trajectories. This requires that the Julia process was started with the -t <NTHREADS> option, or that the JULIA_NUM_THREADS environment variable was set before starting the Julia process.
Lastly, we make use of the GRAPE package's ability to apply box constraints, i.e., an upper_bound and lower_bound for the values of the control pulse.
result = GRAPE.optimize(
trajectories,
tlist;
prop_method = Cheby,
J_T = J_T_sm,
callback = GRAPE.make_grape_print_iters(),
iter_stop = 200,
check_convergence = (res -> ((res.J_T < 1e-2) && "Gate error < 10⁻²")),
upper_bound = 50⋅2π⋅MHz,
lower_bound = -50⋅2π⋅MHz,
use_threads = true,
)[ Info: make_chi for J_T=J_T_sm: fallback to mode=:automatic
[ Info: make_chi for J_T=J_T_sm: automatic with Zygote
iter. J_T ǁ∇Jǁ ǁΔϵǁ ΔJ FG(F) secs
0 8.32e-01 1.21e+00 n/a n/a 1(0) 24.6
1 8.29e-01 1.21e+00 1.21e+00 -3.43e-03 1(0) 0.8
2 6.18e-01 1.39e+00 2.79e-01 -2.10e-01 5(0) 2.8
3 5.46e-01 1.16e+00 3.47e-01 -7.25e-02 2(0) 1.1
4 4.45e-01 1.15e+00 2.49e-01 -1.01e-01 1(0) 0.5
5 4.29e-01 7.21e-01 4.04e-01 -1.61e-02 1(0) 0.6
6 3.92e-01 6.22e-01 2.31e-01 -3.74e-02 1(0) 0.6
7 3.62e-01 3.70e-01 1.08e-01 -2.97e-02 1(0) 0.5
8 3.18e-01 2.80e-01 2.81e-01 -4.36e-02 2(0) 1.1
9 3.01e-01 4.61e-01 1.39e-01 -1.74e-02 1(0) 0.6
10 2.89e-01 3.18e-01 1.15e-01 -1.23e-02 1(0) 0.6
11 2.86e-01 1.26e-01 5.32e-02 -2.47e-03 1(0) 0.6
12 2.81e-01 1.32e-01 1.62e-01 -5.15e-03 1(0) 0.6
13 2.78e-01 2.11e-01 1.56e-01 -3.14e-03 1(0) 0.6
14 2.74e-01 2.24e-01 4.54e-02 -3.60e-03 1(0) 0.6
15 2.71e-01 1.08e-01 5.31e-02 -2.91e-03 1(0) 0.6
16 2.65e-01 8.78e-02 1.79e-01 -6.39e-03 1(0) 0.6
17 2.56e-01 1.43e-01 4.04e-01 -8.95e-03 1(0) 0.6
18 2.46e-01 2.74e-01 2.89e-01 -9.92e-03 1(0) 0.6
19 2.38e-01 2.98e-01 6.18e-01 -8.12e-03 1(0) 0.6
20 2.33e-01 4.16e-01 1.25e-01 -5.33e-03 1(0) 0.6
21 2.27e-01 2.93e-01 5.17e-02 -6.07e-03 1(0) 0.6
22 2.25e-01 1.32e-01 5.54e-02 -1.67e-03 1(0) 0.6
23 2.22e-01 1.24e-01 1.04e-01 -3.32e-03 1(0) 0.5
24 2.08e-01 1.22e-01 4.43e-01 -1.37e-02 2(0) 1.1
25 1.97e-01 2.14e-01 5.25e-01 -1.10e-02 2(0) 1.1
26 1.87e-01 4.06e-01 4.44e-01 -9.51e-03 1(0) 0.6
27 1.72e-01 5.73e-01 5.91e-01 -1.50e-02 1(0) 0.6
28 1.59e-01 4.11e-01 1.13e-01 -1.35e-02 1(0) 0.6
29 1.52e-01 3.06e-01 7.94e-02 -6.79e-03 1(0) 0.6
30 1.41e-01 2.68e-01 2.60e-01 -1.15e-02 1(0) 0.6
31 1.33e-01 3.56e-01 7.29e-02 -7.60e-03 1(0) 0.6
32 1.21e-01 2.39e-01 2.06e-01 -1.14e-02 1(0) 0.6
33 1.12e-01 1.19e-01 2.26e-01 -9.34e-03 1(0) 0.6
34 1.05e-01 8.67e-02 4.97e-01 -6.96e-03 1(0) 0.6
35 9.85e-02 3.66e-01 9.72e-02 -6.73e-03 1(0) 0.6
36 9.79e-02 9.94e-02 1.61e-02 -6.07e-04 1(0) 0.6
37 9.71e-02 6.93e-02 2.27e-02 -7.73e-04 1(0) 0.6
38 9.47e-02 5.66e-02 1.33e-01 -2.40e-03 1(0) 0.6
39 8.99e-02 8.02e-02 2.72e-01 -4.82e-03 1(0) 0.6
40 8.51e-02 1.13e-01 2.28e-01 -4.73e-03 1(0) 0.6
41 7.92e-02 1.28e-01 3.48e-01 -5.97e-03 1(0) 0.6
42 7.54e-02 1.93e-01 9.98e-02 -3.77e-03 1(0) 0.6
43 7.31e-02 8.78e-02 1.06e-01 -2.32e-03 1(0) 0.6
44 6.90e-02 7.74e-02 2.20e-01 -4.05e-03 1(0) 0.6
45 6.17e-02 7.84e-02 6.19e-01 -7.31e-03 1(0) 0.6
46 5.93e-02 1.40e-01 4.22e-02 -2.45e-03 1(0) 0.6
47 5.73e-02 6.04e-02 8.69e-02 -1.95e-03 1(0) 0.6
48 5.60e-02 7.99e-02 8.99e-02 -1.36e-03 1(0) 0.6
49 5.39e-02 7.42e-02 1.75e-01 -2.01e-03 1(0) 0.6
50 5.32e-02 6.95e-02 4.58e-02 -7.09e-04 1(0) 0.6
51 5.17e-02 6.31e-02 1.39e-01 -1.53e-03 1(0) 0.6
52 5.09e-02 6.98e-02 1.37e-01 -8.36e-04 1(0) 0.6
53 4.99e-02 1.34e-01 4.11e-02 -9.72e-04 1(0) 0.6
54 4.93e-02 6.30e-02 5.58e-02 -6.34e-04 1(0) 0.6
55 4.88e-02 4.08e-02 4.14e-02 -4.52e-04 1(0) 0.6
56 4.71e-02 4.99e-02 1.67e-01 -1.71e-03 1(0) 0.6
57 4.45e-02 7.15e-02 2.72e-01 -2.59e-03 1(0) 0.6
58 4.31e-02 1.17e-01 1.98e-01 -1.46e-03 1(0) 0.6
59 4.22e-02 1.10e-01 2.40e-02 -8.38e-04 1(0) 0.6
60 4.19e-02 4.03e-02 3.54e-02 -3.44e-04 1(0) 0.6
61 4.14e-02 4.52e-02 7.28e-02 -5.25e-04 1(0) 0.6
62 4.08e-02 5.69e-02 1.27e-01 -6.04e-04 1(0) 0.6
63 4.02e-02 1.03e-01 2.88e-02 -5.67e-04 1(0) 0.6
64 3.98e-02 4.77e-02 4.69e-02 -3.67e-04 1(0) 0.6
65 3.95e-02 3.78e-02 4.89e-02 -3.09e-04 1(0) 0.6
66 3.88e-02 4.55e-02 1.12e-01 -7.20e-04 1(0) 0.6
67 3.77e-02 5.35e-02 2.84e-01 -1.05e-03 1(0) 0.6
68 3.66e-02 1.54e-01 2.85e-02 -1.10e-03 1(0) 0.6
69 3.61e-02 6.01e-02 4.89e-02 -5.76e-04 1(0) 0.6
70 3.55e-02 4.17e-02 5.47e-02 -5.50e-04 1(0) 0.6
71 3.35e-02 5.94e-02 2.51e-01 -2.00e-03 1(0) 0.6
72 3.20e-02 8.56e-02 2.35e-01 -1.49e-03 1(0) 0.6
73 3.16e-02 4.81e-02 6.36e-02 -3.97e-04 1(0) 0.6
74 3.13e-02 6.19e-02 3.88e-02 -3.32e-04 1(0) 0.6
75 3.04e-02 5.20e-02 1.46e-01 -9.16e-04 1(0) 0.6
76 2.98e-02 5.37e-02 6.97e-02 -5.32e-04 1(0) 0.6
77 2.88e-02 5.76e-02 1.34e-01 -1.05e-03 1(0) 0.6
78 2.67e-02 6.40e-02 2.61e-01 -2.06e-03 1(0) 0.6
79 2.53e-02 6.51e-02 1.66e-01 -1.43e-03 1(0) 0.6
80 2.39e-02 5.94e-02 1.48e-01 -1.38e-03 1(0) 0.6
81 2.22e-02 3.98e-02 2.25e-01 -1.71e-03 1(0) 0.6
82 2.11e-02 5.78e-02 1.38e-01 -1.10e-03 1(0) 0.6
83 2.01e-02 5.70e-02 1.26e-01 -1.01e-03 1(0) 0.6
84 1.94e-02 3.79e-02 1.10e-01 -7.56e-04 1(0) 0.6
85 1.87e-02 3.60e-02 9.71e-02 -6.68e-04 1(0) 0.6
86 1.73e-02 3.44e-02 2.27e-01 -1.37e-03 1(0) 0.6
87 1.69e-02 4.72e-02 7.00e-02 -4.77e-04 1(0) 0.6
88 1.65e-02 3.31e-02 7.40e-02 -3.90e-04 1(0) 0.6
89 1.60e-02 4.15e-02 8.41e-02 -4.68e-04 1(0) 0.6
90 1.55e-02 3.19e-02 9.27e-02 -5.00e-04 1(0) 0.6
91 1.52e-02 2.90e-02 5.72e-02 -3.18e-04 1(0) 0.6
92 1.43e-02 3.13e-02 1.89e-01 -9.14e-04 1(0) 0.6
93 1.36e-02 4.67e-02 1.28e-01 -6.11e-04 1(0) 0.6
94 1.33e-02 3.23e-02 8.86e-02 -3.70e-04 1(0) 0.5
95 1.31e-02 1.98e-02 5.30e-02 -2.14e-04 1(0) 0.6
96 1.27e-02 1.97e-02 9.88e-02 -3.35e-04 1(0) 0.6
97 1.24e-02 3.15e-02 1.18e-01 -3.63e-04 1(0) 0.6
98 1.19e-02 3.20e-02 2.23e-01 -4.99e-04 1(0) 0.5
99 1.17e-02 5.50e-02 2.97e-02 -1.89e-04 1(0) 0.6
100 1.16e-02 2.00e-02 1.98e-02 -5.73e-05 1(0) 0.6
101 1.15e-02 1.95e-02 6.65e-02 -1.50e-04 1(0) 0.6
102 1.13e-02 2.99e-02 4.70e-02 -1.30e-04 1(0) 0.6
103 1.11e-02 2.45e-02 1.01e-01 -2.55e-04 1(0) 0.6
104 1.10e-02 2.15e-02 2.86e-02 -7.96e-05 1(0) 0.5
105 1.07e-02 2.17e-02 1.45e-01 -3.21e-04 1(0) 0.6
106 1.06e-02 3.23e-02 3.84e-02 -1.34e-04 1(0) 0.5
107 1.04e-02 2.74e-02 4.70e-02 -1.50e-04 1(0) 0.6
108 1.01e-02 2.79e-02 8.75e-02 -2.57e-04 1(0) 0.6
109 9.90e-03 2.41e-02 9.03e-02 -2.45e-04 1(0) 0.6GRAPE Optimization Result
-------------------------
- Started at 2025-10-02T22:25:34.123
- Number of trajectories: 4
- Number of iterations: 109
- Number of pure func evals: 0
- Number of func/grad evals: 118
- Value of functional: 9.90083e-03
- Reason for termination: Gate error < 10⁻²
- Ended at 2025-10-02T22:27:03.622 (1 minute, 29 seconds, 499 milliseconds)
Optimization Result
The call to GRAPE.optimize returns a results-object from which we can extract the optimized controls:
ϵ_opt = result.optimized_controls[1] + 𝕚 * result.optimized_controls[2];While the guess controls ϵ_re_guess and ϵ_im_guess that we defined when setting up the system Hamiltonian could have been either functions (as they were) or vectors of pulse values, the GRAPE method inherently discretizes all control to a time grid (the method is piecewise-constant, by definition). Thus, the optimized_controls are vectors of control values for each point in tlist
Also remember that we used a ShapedAmplitude when setting up the Hamiltonian, with the definition $Ω(t) = S(t)ϵ(t)$. Thus, to get the physical optimized control field, we need to multiply with that same shape S(t), now also discretized to the time grid.
Ω_opt = ϵ_opt .* discretize(Ω_re_guess.shape, tlist)
fig = plot_complex_pulse(tlist, Ω_opt)
We can see how the optimization has tuned both the amplitude and the complex phase of the microwave field, while retaining the overall shape $S(t)$. We can also simulate again the dynamics of the $|10⟩$ state under the optimized fields:
Ĥ_opt = transmon_hamiltonian(;
δ₁ = 100⋅2π⋅MHz,
δ₂ = -100⋅2π⋅MHz,
J = 3⋅2π⋅MHz,
Ω_re = real.(Ω_opt),
Ω_im = imag.(Ω_opt),
)
pops = propagate(
basis[3], Ĥ_opt, tlist; method=Cheby, storage=true, observables=(Ψ -> Array(abs2.(Ψ)), )
)
fig = plot(tlist ./ ns, pops'; label=["00" "01" "10" "11"], xlabel="time (ns)", legend=:outertop, legend_column=-1)
As expected (and demanded by the "conditional NOT"), $|10⟩$ now evolves into $|11⟩$.