# NBVAL_IGNORE_OUTPUT
%load_ext watermark
import os
import qutip
import numpy as np
import scipy
import matplotlib
import matplotlib.pylab as plt
import krotov
from scipy.fftpack import fft
from scipy.interpolate import interp1d
%watermark -v --iversions
krotov 0.3.0+dev matplotlib.pylab 1.16.4 qutip 4.3.1 scipy 1.2.1 numpy 1.16.4 matplotlib 3.1.0 CPython 3.7.3 IPython 7.5.0
$\newcommand{tr}[0]{\operatorname{tr}} \newcommand{diag}[0]{\operatorname{diag}} \newcommand{abs}[0]{\operatorname{abs}} \newcommand{pop}[0]{\operatorname{pop}} \newcommand{aux}[0]{\text{aux}} \newcommand{opt}[0]{\text{opt}} \newcommand{tgt}[0]{\text{tgt}} \newcommand{init}[0]{\text{init}} \newcommand{lab}[0]{\text{lab}} \newcommand{rwa}[0]{\text{rwa}} \newcommand{bra}[1]{\langle#1\vert} \newcommand{ket}[1]{\vert#1\rangle} \newcommand{Bra}[1]{\left\langle#1\right\vert} \newcommand{Ket}[1]{\left\vert#1\right\rangle} \newcommand{Braket}[2]{\left\langle #1\vphantom{#2} \mid #2\vphantom{#1}\right\rangle} \newcommand{op}[1]{\hat{#1}} \newcommand{Op}[1]{\hat{#1}} \newcommand{dd}[0]{\,\text{d}} \newcommand{Liouville}[0]{\mathcal{L}} \newcommand{DynMap}[0]{\mathcal{E}} \newcommand{identity}[0]{\mathbf{1}} \newcommand{Norm}[1]{\lVert#1\rVert} \newcommand{Abs}[1]{\left\vert#1\right\vert} \newcommand{avg}[1]{\langle#1\rangle} \newcommand{Avg}[1]{\left\langle#1\right\rangle} \newcommand{AbsSq}[1]{\left\vert#1\right\vert^2} \newcommand{Re}[0]{\operatorname{Re}} \newcommand{Im}[0]{\operatorname{Im}}$
In the previous examples, we have only optimized for state-to-state transitions, i.e., for a single objective. This example shows the optimization of a simple quantum gate, which requires multiple objectives to be fulfilled simultaneously (one for each state in the logical basis). We consider a superconducting "transmon" qubit and implement a single-qubit Pauli-X gate.
The effective Hamiltonian of a single transmon depends on the capacitive energy $E_C=e^2/2C$ and the Josephson energy $E_J$, an energy due to the Josephson junction working as a nonlinear inductor periodic with the flux $\Phi$. In the so-called transmon limit, the ratio between these two energies lies around $E_J / E_C \approx 45$. The Hamiltonian for the transmon is
$$ \op{H}_{0} = 4 E_C (\hat{n}-n_g)^2 - E_J \cos(\hat{\Phi}) $$where $\hat{n}$ is the number operator, which counts the relative number of Cooper pairs capacitively stored in the junction, and $n_g$ is the effective offset charge measured in Cooper pair charge units. The equation can be written in a truncated charge basis defined by the number operator $\op{n} \ket{n} = n \ket{n}$ such that
$$ \op{H}_{0} = 4 E_C \sum_{j=-N} ^N (j-n_g)^2 \ket{j} \bra{j} - \frac{E_J}{2} \sum_{j=-N}^{N-1} (\ket{j+1} \bra{j} + \ket{j} \bra{j+1}). $$A voltage $V(t)$ applied to the circuit couples to the charge Hamiltonian $\op{q}$, which in the (truncated) charge basis reads
$$ \op{H}_1 = \op{q} = \sum_{j=-N} ^N -2n \ket{n} \bra{n}\,. $$The factor 2 is due to the charge carriers in a superconductor being Cooper pairs. The total Hamiltonian is
$$ \op{H} = \op{H}_{0} + V(t) \cdot \op{H}_{1} $$We use a Gaussian voltage profile as the guess pulse:
tlist = np.linspace(0, 10, 1000)
def eps0(t, args):
T = tlist[-1]
return 4 * np.exp(-40.0 * (t / T - 0.5) ** 2)
def plot_pulse(pulse, tlist, xlimit=None):
fig, ax = plt.subplots()
if callable(pulse):
pulse = np.array([pulse(t, None) for t in tlist])
ax.plot(tlist, pulse)
ax.set_xlabel('time (ns)')
ax.set_ylabel('pulse amplitude')
if xlimit is not None:
ax.set_xlim(xlimit)
plt.show(fig)
plot_pulse(eps0, tlist)
The complete Hamiltonian is instantiated as
def transmon_hamiltonian(Ec=0.386, EjEc=45, nstates=8, ng=0.0, T=10.0):
"""Transmon Hamiltonian
Args:
Ec: capacitive energy
EjEc: ratio `Ej` / `Ec`
nstates: defines the maximum and minimum states for the basis. The
truncated basis will have a total of ``2*nstates + 1`` states
ng: offset charge
T: gate duration
"""
Ej = EjEc * Ec
n = np.arange(-nstates, nstates + 1)
up = np.diag(np.ones(2 * nstates), k=-1)
do = up.T
H0 = qutip.Qobj(np.diag(4 * Ec * (n - ng) ** 2) - Ej * (up + do) / 2.0)
H1 = qutip.Qobj(-2 * np.diag(n))
return [H0, [H1, eps0]]
H = transmon_hamiltonian()
We define the logical basis $\ket{0_l}$ and $\ket{1_l}$ (not to be confused with the charge states $\ket{n=0}$ and $\ket{n=1}$) as the eigenstates of the drift Hamiltonian $\op{H}_0$ with the lowest energy. The optimization goal is to find a potential $V_{opt}(t)$ such that after a given final time $T$ implements an X-gate on this logical basis.
def logical_basis(H):
H0 = H[0]
eigenvals, eigenvecs = scipy.linalg.eig(H0.full())
ndx = np.argsort(eigenvals.real)
E = eigenvals[ndx].real
V = eigenvecs[:, ndx]
psi0 = qutip.Qobj(V[:, 0])
psi1 = qutip.Qobj(V[:, 1])
w01 = E[1] - E[0] # Transition energy between states
print("Energy of qubit transition is %.3f" % w01)
return psi0, psi1
psi0, psi1 = logical_basis(H)
Energy of qubit transition is 6.914
We also introduce the projectors $P_i = \ket{\psi _i}\bra{\psi _i}$ for the logical states $\ket{\psi _i} \in \{\ket{0_l}, \ket{1_l}\}$
proj0 = qutip.ket2dm(psi0)
proj1 = qutip.ket2dm(psi1)
The key insight for the realization of a quantum gate $\Op{O}$ is that (by virtue of linearity)
$$\ket{\Psi(t=0)} \rightarrow \ket{\Psi(t=T)} = \Op{U}(T, \epsilon(t))\ket{\Psi(0)} = \Op{O} \ket{\Psi(0)} $$is fulfilled for an arbitrary state $\Ket{\Psi(t=0)}$ if an only if $\Op{U}(T, \epsilon(t))\ket{k} = \Op{O} \ket{k}$ for every state $\ket{k}$ in logical basis, for the time evolution operator $\Op{U}(T, \epsilon(t))$ from $t=0$ to $t=T$ under the same control $\epsilon(t)$.
The function krotov.gate_objectives
automatically sets up the corresponding
objectives $\forall \ket{k}: \ket{k} \rightarrow \Op{O} \ket{k}$:
objectives = krotov.gate_objectives(
basis_states=[psi0, psi1], gate=qutip.operators.sigmax(), H=H
)
objectives
[Objective[|Ψ₀(17)⟩ to |Ψ₁(17)⟩ via [H₀[17,17], [H₁[17,17], u₁(t)]]], Objective[|Ψ₁(17)⟩ to |Ψ₀(17)⟩ via [H₀[17,17], [H₁[17,17], u₁(t)]]]]
guess_dynamics = [
objectives[x].mesolve(tlist, e_ops=[proj0, proj1]) for x in [0, 1]
]
def plot_population(result):
'''Representation of the expected values for the initial states'''
fig, ax = plt.subplots()
ax.plot(result.times, result.expect[0], label='0')
ax.plot(result.times, result.expect[1], label='1')
ax.legend()
ax.set_xlabel('time')
ax.set_ylabel('population')
plt.show(fig)
plot_population(guess_dynamics[0])
plot_population(guess_dynamics[1])
We define the desired shape of the update and the factor $\lambda_a$, and then start the optimization
def S(t):
"""Scales the Krotov methods update of the pulse value at the time t"""
return krotov.shapes.flattop(
t, t_start=0.0, t_stop=10.0, t_rise=0.5, func='sinsq'
)
pulse_options = {H[1][1]: dict(lambda_a=1, update_shape=S)}
oct_result = krotov.optimize_pulses(
objectives,
pulse_options,
tlist,
propagator=krotov.propagators.expm,
chi_constructor=krotov.functionals.chis_re,
info_hook=krotov.info_hooks.print_table(
J_T=krotov.functionals.J_T_re,
show_g_a_int_per_pulse=True,
unicode=False,
),
check_convergence=krotov.convergence.Or(
krotov.convergence.value_below(1e-3, name='J_T'),
krotov.convergence.delta_below(1e-5),
krotov.convergence.check_monotonic_error,
),
iter_stop=5,
parallel_map=(
qutip.parallel_map,
qutip.parallel_map,
krotov.parallelization.parallel_map_fw_prop_step,
),
)
iter. J_T g_a_int J Delta_J_T Delta J secs 0 1.00e+00 0.00e+00 1.00e+00 n/a n/a 1 1 2.80e-01 3.41e-01 6.22e-01 -7.20e-01 -3.78e-01 5 2 2.12e-01 3.06e-02 2.43e-01 -6.81e-02 -3.75e-02 4 3 1.35e-01 3.28e-02 1.68e-01 -7.72e-02 -4.44e-02 4 4 9.79e-02 1.56e-02 1.13e-01 -3.71e-02 -2.15e-02 4 5 7.13e-02 1.11e-02 8.25e-02 -2.65e-02 -1.54e-02 4
(this takes a while ...)
dumpfile = "./transmonxgate_oct_result.dump"
if os.path.isfile(dumpfile):
oct_result = krotov.result.Result.load(dumpfile, objectives)
else:
oct_result = krotov.optimize_pulses(
objectives,
pulse_options,
tlist,
propagator=krotov.propagators.expm,
chi_constructor=krotov.functionals.chis_re,
info_hook=krotov.info_hooks.print_table(
J_T=krotov.functionals.J_T_re,
show_g_a_int_per_pulse=True,
unicode=False,
),
check_convergence=krotov.convergence.Or(
krotov.convergence.value_below(1e-3, name='J_T'),
krotov.convergence.delta_below(1e-5),
krotov.convergence.check_monotonic_error,
),
iter_stop=1000,
parallel_map=(
qutip.parallel_map,
qutip.parallel_map,
krotov.parallelization.parallel_map_fw_prop_step,
),
continue_from=oct_result
)
oct_result.dump(dumpfile)
oct_result
Krotov Optimization Result -------------------------- - Started at 2019-04-12 17:45:43 - Number of objectives: 2 - Number of iterations: 398 - Reason for termination: Reached convergence: Δ(('info_vals', T[-1]),('info_vals', T[-2])) < 1e-05 - Ended at 2019-04-12 18:20:47
def plot_convergence(result):
fig, ax = plt.subplots()
ax.semilogy(result.iters, np.array(result.info_vals))
ax.set_xlabel('OCT iteration')
ax.set_ylabel('error')
plt.show(fig)
plot_convergence(oct_result)
We obtain the following optimized pulse:
plot_pulse(oct_result.optimized_controls[0], tlist)
The oscillations in the control shape indicate non-negligible spectral broadening:
def plot_spectrum(pulse, tlist, xlim=None):
if callable(pulse):
pulse = np.array([pulse(t, None) for t in tlist])
dt = tlist[1] - tlist[0]
n = len(tlist)
w = np.fft.fftfreq(n, d=dt/(2.0*np.pi))
# the factor 2π in the normalization means that
# the spectrum is in units of angular frequency,
# which is normally what we want
spectrum = np.fft.fft(pulse) / n
# normalizing the spectrum with n means that
# the y-axis is independent of dt
# we assume a real-valued pulse, so we throw away
# the half of the spectrum with negative frequencies
w = w[range(int(n / 2))]
spectrum = np.abs(spectrum[range(int(n / 2))])
fig, ax = plt.subplots()
ax.plot(w, spectrum, '-o')
ax.set_xlabel(r'$\omega$')
ax.set_ylabel('amplitude (arb. units)')
if xlim is not None:
ax.set_xlim(*xlim)
plt.show(fig)
plot_spectrum(oct_result.optimized_controls[0], tlist, xlim=(0, 40))
Lastly, we verify that the pulse produces the desired dynamics $\ket{0_l} \rightarrow \ket{1_l}$ and $\ket{1_l} \rightarrow \ket{0_l}$:
opt_dynamics = [
oct_result.optimized_objectives[x].mesolve(tlist, e_ops=[proj0, proj1])
for x in [0, 1]
]
opt_states = [
oct_result.optimized_objectives[x].mesolve(tlist) for x in [0, 1]
]
plot_population(opt_dynamics[0])
plot_population(opt_dynamics[1])
Looking only at the population dynamics does not verify that the we have a coherent gate, i.e. that the basis transforms with the same global phase.
def relative_phase(result):
overlap_0 = np.angle(result[0].states[-1].overlap(psi1))
overlap_1 = np.angle(result[1].states[-1].overlap(psi0))
rel_phase = (overlap_0 - overlap_1) % (2 * np.pi)
print("Final relative phase / 2π = %.2e" % rel_phase)
relative_phase(opt_states)
Final relative phase / 2π = 5.00e-03
Since the optimized pulse shows some oscillations (cf. the spectrum above),
it is a good idea to check for any discretization error. To this end, we also
propagate the optimization result using the same propagator that
was used in the optimization (instead of qutip.mesolve
). The main difference
between the two propagations is that mesolve
assumes piecewise constant
pulses that switch between two points in tlist
, whereas propagate
assumes
that pulses are constant on the intervals of tlist
, and thus switches on
the points in tlist
.
opt_dynamics2 = [
oct_result.optimized_objectives[x].propagate(
tlist, e_ops=[proj0, proj1], propagator=krotov.propagators.expm
)
for x in [0, 1]
]
The difference between the two propagations gives an indication of the error due to the choice of the piecewise constant time discretization. If this error were unacceptably large, we would need a smaller time step.
print(
"Time discretization error = %.2e" %
abs(opt_dynamics2[0].expect[1][-1] - opt_dynamics[0].expect[1][-1])
)
Time discretization error = 2.86e-05