Setup of a two-qubit chip with \(C^3\)¶
In this example we will set-up a two qubit quantum processor and define a simple gate.
Imports¶
# System imports
import copy
import numpy as np
import time
import itertools
import matplotlib.pyplot as plt
import tensorflow as tf
import tensorflow_probability as tfp
# Main C3 objects
from c3.c3objs import Quantity as Qty
from c3.parametermap import ParameterMap as PMap
from c3.experiment import Experiment as Exp
from c3.model import Model as Mdl
from c3.generator.generator import Generator as Gnr
# Building blocks
import c3.generator.devices as devices
import c3.signal.gates as gates
import c3.libraries.chip as chip
import c3.signal.pulse as pulse
import c3.libraries.tasks as tasks
# Libs and helpers
import c3.libraries.algorithms as algorithms
import c3.libraries.hamiltonians as hamiltonians
import c3.libraries.fidelities as fidelities
import c3.libraries.envelopes as envelopes
import c3.utils.qt_utils as qt_utils
import c3.utils.tf_utils as tf_utils
Model components¶
We first create a qubit. Each parameter is a Quantity (Qty()) object
with bounds and a unit. In \(C^3\), the default multi-level qubit is
a Transmon modelled as a Duffing oscillator with frequency
\(\omega\) and anharmonicity \(\delta\) :
The “name” will be used to identify this qubit (or other component) later and should thus be chosen carefully.
qubit_lvls = 3
freq_q1 = 5e9
anhar_q1 = -210e6
t1_q1 = 27e-6
t2star_q1 = 39e-6
qubit_temp = 50e-3
q1 = chip.Qubit(
name="Q1",
desc="Qubit 1",
freq=Qty(
value=freq_q1,
min_val=4.995e9 ,
max_val=5.005e9 ,
unit='Hz 2pi'
),
anhar=Qty(
value=anhar_q1,
min_val=-380e6 ,
max_val=-120e6 ,
unit='Hz 2pi'
),
hilbert_dim=qubit_lvls,
t1=Qty(
value=t1_q1,
min_val=1e-6,
max_val=90e-6,
unit='s'
),
t2star=Qty(
value=t2star_q1,
min_val=10e-6,
max_val=90e-3,
unit='s'
),
temp=Qty(
value=qubit_temp,
min_val=0.0,
max_val=0.12,
unit='K'
)
)
And the same for a second qubit.
freq_q2 = 5.6e9
anhar_q2 = -240e6
t1_q2 = 23e-6
t2star_q2 = 31e-6
q2 = chip.Qubit(
name="Q2",
desc="Qubit 2",
freq=Qty(
value=freq_q2,
min_val=5.595e9 ,
max_val=5.605e9 ,
unit='Hz 2pi'
),
anhar=Qty(
value=anhar_q2,
min_val=-380e6 ,
max_val=-120e6 ,
unit='Hz 2pi'
),
hilbert_dim=qubit_lvls,
t1=Qty(
value=t1_q2,
min_val=1e-6,
max_val=90e-6,
unit='s'
),
t2star=Qty(
value=t2star_q2,
min_val=10e-6,
max_val=90e-6,
unit='s'
),
temp=Qty(
value=qubit_temp,
min_val=0.0,
max_val=0.12,
unit='K'
)
)
A static coupling between the two is realized in the following way. We
supply the type of coupling by selecting int_XX
\((b_1+b_1^\dagger)(b_2+b_2^\dagger)\) from the hamiltonian library.
The “connected” property contains the list of qubit names to be coupled,
in this case “Q1” and “Q2”.
coupling_strength = 20e6
q1q2 = chip.Coupling(
name="Q1-Q2",
desc="coupling",
comment="Coupling qubit 1 to qubit 2",
connected=["Q1", "Q2"],
strength=Qty(
value=coupling_strength,
min_val=-1 * 1e3 ,
max_val=200e6 ,
unit='Hz 2pi'
),
hamiltonian_func=hamiltonians.int_XX
)
In the same spirit, we specify control Hamiltonians to drive the system. Again “connected” connected tells us which qubit this drive acts on and “name” will later be used to assign the correct control signal to this drive line.
drive = chip.Drive(
name="d1",
desc="Drive 1",
comment="Drive line 1 on qubit 1",
connected=["Q1"],
hamiltonian_func=hamiltonians.x_drive
)
drive2 = chip.Drive(
name="d2",
desc="Drive 2",
comment="Drive line 2 on qubit 2",
connected=["Q2"],
hamiltonian_func=hamiltonians.x_drive
)
SPAM errors¶
In experimental practice, the qubit state can be mis-classified during read-out. We simulate this by constructing a confusion matrix, containing the probabilities for one qubit state being mistaken for another.
m00_q1 = 0.97 # Prop to read qubit 1 state 0 as 0
m01_q1 = 0.04 # Prop to read qubit 1 state 0 as 1
m00_q2 = 0.96 # Prop to read qubit 2 state 0 as 0
m01_q2 = 0.05 # Prop to read qubit 2 state 0 as 1
one_zeros = np.array([0] * qubit_lvls)
zero_ones = np.array([1] * qubit_lvls)
one_zeros[0] = 1
zero_ones[0] = 0
val1 = one_zeros * m00_q1 + zero_ones * m01_q1
val2 = one_zeros * m00_q2 + zero_ones * m01_q2
min_val = one_zeros * 0.8 + zero_ones * 0.0
max_val = one_zeros * 1.0 + zero_ones * 0.2
confusion_row1 = Qty(value=val1, min_val=min_val, max_val=max_val, unit="")
confusion_row2 = Qty(value=val2, min_val=min_val, max_val=max_val, unit="")
conf_matrix = tasks.ConfusionMatrix(Q1=confusion_row1, Q2=confusion_row2)
The following task creates an initial thermal state with given temperature.
init_temp = 50e-3
init_ground = tasks.InitialiseGround(
init_temp=Qty(
value=init_temp,
min_val=-0.001,
max_val=0.22,
unit='K'
)
)
We collect the parts specified above in the Model.
model = Mdl(
[q1, q2], # Individual, self-contained components
[drive, drive2, q1q2], # Interactions between components
[conf_matrix, init_ground] # SPAM processing
)
Further, we can decide between coherent or open-system dynamics using set_lindbladian() and whether to eliminate the static coupling by going to the dressed frame with set_dressed().
model.set_lindbladian(False)
model.set_dressed(True)
Control signals¶
With the system model taken care of, we now specify the control electronics and signal chain. Complex shaped controls are often realized by creating an envelope signal with an arbitrary waveform generator (AWG) with limited bandwith and mixing it with a fast, stable local oscillator (LO).
sim_res = 100e9 # Resolution for numerical simulation
awg_res = 2e9 # Realistic, limited resolution of an AWG
lo = devices.LO(name='lo', resolution=sim_res)
awg = devices.AWG(name='awg', resolution=awg_res)
mixer = devices.Mixer(name='mixer')
Waveform generators exhibit a rise time, the time it takes until the target voltage is set. This has a smoothing effect on the resulting pulse shape.
resp = devices.Response(
name='resp',
rise_time=Qty(
value=0.3e-9,
min_val=0.05e-9,
max_val=0.6e-9,
unit='s'
),
resolution=sim_res
)
In simulation, we translate between AWG resolution and simulation (or “analog”) resolution by including an up-sampling device.
dig_to_an = devices.DigitalToAnalog(
name="dac",
resolution=sim_res
)
Control electronics apply voltages to lines, whereas in a Hamiltonian we usually write the control fields in energy or frequency units. In practice, this conversion can be highly non-trivial if it involves multiple stages of attenuation and for example the conversion of a line voltage in an antenna to a dipole field coupling to the qubit. The following device represents a simple, linear conversion factor.
v2hz = 1e9
v_to_hz = devices.VoltsToHertz(
name='v_to_hz',
V_to_Hz=Qty(
value=v2hz,
min_val=0.9e9,
max_val=1.1e9,
unit='Hz/V'
)
)
The generator combines the parts of the signal generation and assigns a signal chain to each control line.
generator = Gnr(
devices={
"LO": devices.LO(name='lo', resolution=sim_res, outputs=1),
"AWG": devices.AWG(name='awg', resolution=awg_res, outputs=1),
"DigitalToAnalog": devices.DigitalToAnalog(
name="dac",
resolution=sim_res,
inputs=1,
outputs=1
),
"Response": devices.Response(
name='resp',
rise_time=Qty(
value=0.3e-9,
min_val=0.05e-9,
max_val=0.6e-9,
unit='s'
),
resolution=sim_res,
inputs=1,
outputs=1
),
"Mixer": devices.Mixer(name='mixer', inputs=2, outputs=1),
"VoltsToHertz": devices.VoltsToHertz(
name='v_to_hz',
V_to_Hz=Qty(
value=1e9,
min_val=0.9e9,
max_val=1.1e9,
unit='Hz/V'
),
inputs=1,
outputs=1
)
},
chains= {
"d1": {
"LO": [],
"AWG": [],
"DigitalToAnalog": ["AWG"],
"Response": ["DigitalToAnalog"],
"Mixer": ["LO", "Response"],
"VoltsToHertz": ["Mixer"]
},
"d2": {
"LO": [],
"AWG": [],
"DigitalToAnalog": ["AWG"],
"Response": ["DigitalToAnalog"],
"Mixer": ["LO", "Response"],
"VoltsToHertz": ["Mixer"]
}
}
)
Gates-set and Parameter map¶
It remains to write down what kind of operations we want to perform on the device. For a gate based quantum computing chip, we define a gate-set.
We choose a gate time of 7ns and a Gaussian envelope shape with a list of parameters.
t_final = 7e-9 # Time for single qubit gates
sideband = 50e6
gauss_params_single = {
'amp': Qty(
value=0.5,
min_val=0.4,
max_val=0.6,
unit="V"
),
't_final': Qty(
value=t_final,
min_val=0.5 * t_final,
max_val=1.5 * t_final,
unit="s"
),
'sigma': Qty(
value=t_final / 4,
min_val=t_final / 8,
max_val=t_final / 2,
unit="s"
),
'xy_angle': Qty(
value=0.0,
min_val=-0.5 * np.pi,
max_val=2.5 * np.pi,
unit='rad'
),
'freq_offset': Qty(
value=-sideband - 3e6 ,
min_val=-56 * 1e6 ,
max_val=-52 * 1e6 ,
unit='Hz 2pi'
),
'delta': Qty(
value=-1,
min_val=-5,
max_val=3,
unit=""
)
}
Here we take gaussian_nonorm() from the libraries as the function to
define the shape.
gauss_env_single = pulse.Envelope(
name="gauss",
desc="Gaussian comp for single-qubit gates",
params=gauss_params_single,
shape=envelopes.gaussian_nonorm
)
We also define a gate that represents no driving.
nodrive_env = pulse.Envelope(
name="no_drive",
params={
't_final': Qty(
value=t_final,
min_val=0.5 * t_final,
max_val=1.5 * t_final,
unit="s"
)
},
shape=envelopes.no_drive
)
We specify the drive tones with an offset from the qubit frequencies. As is done in experiment, we will later adjust the resonance by modulating the envelope function.
lo_freq_q1 = 5e9 + sideband
carrier_parameters = {
'freq': Qty(
value=lo_freq_q1,
min_val=4.5e9 ,
max_val=6e9 ,
unit='Hz 2pi'
),
'framechange': Qty(
value=0.0,
min_val= -np.pi,
max_val= 3 * np.pi,
unit='rad'
)
}
carr = pulse.Carrier(
name="carrier",
desc="Frequency of the local oscillator",
params=carrier_parameters
)
For the second qubit drive tone, we copy the first one and replace the frequency. The deepcopy is to ensure that we don’t just create a pointer to the first drive.
lo_freq_q2 = 5.6e9 + sideband
carr_2 = copy.deepcopy(carr)
carr_2.params['freq'].set_value(lo_freq_q2)
Instructions¶
We define the gates we want to perform with a “name” that will identify them later and “channels” relating to the control Hamiltonians and drive lines we specified earlier. As a start we write down 90 degree rotations in the positive \(x\)-direction and identity gates for both qubits. Then we add a carrier and envelope to each.
rx90p_q1 = gates.Instruction(
name="rx90p", targets=[0], t_start=0.0, t_end=t_final, channels=["d1", "d2"]
)
rx90p_q2 = gates.Instruction(
name="rx90p", targets=[1], t_start=0.0, t_end=t_final, channels=["d1", "d2"]
)
rx90p_q1.add_component(gauss_env_single, "d1")
rx90p_q1.add_component(carr, "d1")
rx90p_q2.add_component(copy.deepcopy(gauss_env_single), "d2")
rx90p_q2.add_component(carr_2, "d2")
When later compiling gates into sequences, we have to take care of the relative rotating frames of the qubits and local oscillators. We do this by adding a phase after each gate that realigns the frames.
rx90p_q1.add_component(nodrive_env, "d2")
rx90p_q1.add_component(copy.deepcopy(carr_2), "d2")
rx90p_q1.comps["d2"]["carrier"].params["framechange"].set_value(
(-sideband * t_final) * 2 * np.pi % (2 * np.pi)
)
rx90p_q2.add_component(nodrive_env, "d1")
rx90p_q2.add_component(copy.deepcopy(carr), "d1")
rx90p_q2.comps["d1"]["carrier"].params["framechange"].set_value(
(-sideband * t_final) * 2 * np.pi % (2 * np.pi)
)
The remainder of the gates-set can be derived from the RX90p gate by shifting its phase by multiples of \(\pi/2\).
ry90p_q1 = copy.deepcopy(rx90p_q1)
ry90p_q1.name = "ry90p"
rx90m_q1 = copy.deepcopy(rx90p_q1)
rx90m_q1.name = "rx90m"
ry90m_q1 = copy.deepcopy(rx90p_q1)
ry90m_q1.name = "ry90m"
ry90p_q1.comps['d1']['gauss'].params['xy_angle'].set_value(0.5 * np.pi)
rx90m_q1.comps['d1']['gauss'].params['xy_angle'].set_value(np.pi)
ry90m_q1.comps['d1']['gauss'].params['xy_angle'].set_value(1.5 * np.pi)
single_q_gates = [rx90p_q1, ry90p_q1, rx90m_q1, ry90m_q1]
ry90p_q2 = copy.deepcopy(rx90p_q2)
ry90p_q2.name = "ry90p"
rx90m_q2 = copy.deepcopy(rx90p_q2)
rx90m_q2.name = "rx90m"
ry90m_q2 = copy.deepcopy(rx90p_q2)
ry90m_q2.name = "ry90m"
ry90p_q2.comps['d2']['gauss'].params['xy_angle'].set_value(0.5 * np.pi)
rx90m_q2.comps['d2']['gauss'].params['xy_angle'].set_value(np.pi)
ry90m_q2.comps['d2']['gauss'].params['xy_angle'].set_value(1.5 * np.pi)
single_q_gates.extend([rx90p_q2, ry90p_q2, rx90m_q2, ry90m_q2])
With every component defined, we collect them in the parameter map, our object that holds information and methods to manipulate and examine model and control parameters.
parameter_map = PMap(instructions=all_1q_gates_comb, model=model, generator=generator)
The experiment¶
Finally everything is collected in the experiment that provides the functions to interact with the system.
exp = Exp(pmap=parameter_map)
Simulation¶
With our experiment all set-up, we can perform simulations. We first decide which basic gates to simulate, in this case only the 90 degree rotation on one qubit and the identity.
exp.set_opt_gates(['RX90p:Id', 'Id:Id'])
The actual numerical simulation is done by calling exp.compute_propagators().
This is the most resource intensive part as it involves solving the
equations of motion for the system.
unitaries = exp.compute_propagators()
Dynamics¶
To investigate dynamics, we define the ground state as an initial state.
psi_init = [[0] * 9]
psi_init[0][0] = 1
init_state = tf.transpose(tf.constant(psi_init, tf.complex128))
init_state
<tf.Tensor: shape=(9, 1), dtype=complex128, numpy=
array([[1.+0.j],
[0.+0.j],
[0.+0.j],
[0.+0.j],
[0.+0.j],
[0.+0.j],
[0.+0.j],
[0.+0.j],
[0.+0.j]])>
Since we stored the process matrices, we can now relatively inexpensively evaluate sequences. We start with just one gate
barely_a_seq = ['rx90p[0]']
and plot system dynamics.
def plot_dynamics(exp, psi_init, seq, goal=-1):
"""
Plotting code for time-resolved populations.
Parameters
----------
psi_init: tf.Tensor
Initial state or density matrix.
seq: list
List of operations to apply to the initial state.
goal: tf.float64
Value of the goal function, if used.
debug: boolean
If true, return a matplotlib figure instead of saving.
"""
model = exp.pmap.model
dUs = exp.partial_propagators
psi_t = psi_init.numpy()
pop_t = exp.populations(psi_t, model.lindbladian)
for gate in seq:
for du in dUs[gate]:
psi_t = np.matmul(du.numpy(), psi_t)
pops = exp.populations(psi_t, model.lindbladian)
pop_t = np.append(pop_t, pops, axis=1)
fig, axs = plt.subplots(1, 1)
ts = exp.ts
dt = ts[1] - ts[0]
ts = np.linspace(0.0, dt*pop_t.shape[1], pop_t.shape[1])
axs.plot(ts / 1e-9, pop_t.T)
axs.grid(linestyle="--")
axs.tick_params(
direction="in", left=True, right=True, top=True, bottom=True
)
axs.set_xlabel('Time [ns]')
axs.set_ylabel('Population')
plt.legend(model.state_labels)
pass
plot_dynamics(exp, init_state, barely_a_seq)
We can see an ill-defined un-optimized gate. The labels indicate qubit states in the product basis. Next we increase the number of repetitions of the same gate.
barely_a_seq * 10
['rx90p[0]',
'rx90p[0]',
'rx90p[0]',
'rx90p[0]',
'rx90p[0]',
'rx90p[0]',
'rx90p[0]',
'rx90p[0]',
'rx90p[0]',
'rx90p[0]']
plot_dynamics(exp, init_state, barely_a_seq * 5)
plot_dynamics(exp, init_state, barely_a_seq * 10)
Note that at this point, we only multiply already computed matrices. We don’t need to solve the equations of motion again for new sequences.