# Control hardware: system identification

**Characterizing quantum hardware with the BOULDER OPAL optimization engine**

The BOULDER OPAL optimization engine provides a large, modular collection of configuration options that allows it to be used for more applications than finding optimal and robust pulses. In this Application note, we show how the optimization engine can be used to determine the parameters that characterize some control hardware. We illustrate this with two examples: first by identifying a linear filter that is present in a single-qubit control line, and then by identifying a constant Hamiltonian to which a qubit is subject. The procedure consists of repeatedly measuring a qubit after applying a set of different pulses or after letting it evolve for different wait times. The results of these measurements are then given to the optimizer, which finds the parameters most likely to have generated those results.

*This notebook contains preview features that are not currently available in BOULDER OPAL. Please contact us if you want a live demonstration for your hardware. Preview features are included through the package qctrlcore, which is not publicly available: cells with qctrlcore will not be executable.*

```
# Essential imports
import numpy as np
from scipy.linalg import expm
import qctrlcore as qc
import tensorflow as tf
# Plotting imports and setup
from qctrlvisualizer import get_qctrl_style
import matplotlib.pyplot as plt
plt.style.use(get_qctrl_style())
# Define standard matrices
sigma_x = np.array([[0, 1], [1, 0]], dtype=np.complex)
sigma_y = np.array([[0, -1j], [1j, 0]], dtype=np.complex)
sigma_z = np.array([[1, 0], [0, -1]], dtype=np.complex)
```

## Characterizing the control lines

In this first section, we will consider the identification of a linear filter that alters the shape of the pulses in the control line. By feeding different control pulses and then measuring the qubit, we are able to extract information about the parameters that characterize the filter. This setup is illustrated in the figure below.

To exemplify this procedure, we will consider a one-qubit system obeying the following Hamiltonian:

$$ H(t) = \frac{1}{2} \mathcal{L} \left[ \alpha(t) \right] \sigma_x, $$where $\alpha(t)$ is a real time-dependent control pulse and $\mathcal{L}$ is a filter applied to the pulse. This linear filter acts on the control pulse $\alpha(t)$ according to the convolution product with some kernel $K(t)$:

$$ \mathcal{L} \left[ \alpha(t) \right] = \int_0^t \mathrm{d} s \; \alpha(s) K(t-s). $$To characterize this filter, we want to determine the parameters that provide the form of the kernel $K(t)$. We will consider a Gaussian kernel characterized by two unknown parameters, a standard deviation $\sigma$ and a mean $\mu$:

$$ K(t) = \frac{1}{\sqrt{2\pi \sigma^2}} \exp \left\{ - \frac{(t-\mu)^2}{2\sigma^2} \right\}. $$The standard deviation $\sigma$ roughly represents the range of frequencies that the filter allows to pass, while the mean $\mu$ represents a time offset between when the pulses were intended to be generated and when they are actually produced. To estimate these two parameters, we will simulate the effect of this filter on a series of pulses and then use the simulated results to find the most likely form of the filter's kernel.

### Simulating the input data

To estimate the parameters $\sigma$ and $\mu$, we subject our system to a series of Gaussian pulses with different standard deviations and timings, and measure their effect on the qubit. Each of these Gaussian wave packets will contain a different range of frequencies, and thus be differently affected by the linear filter. The different timings also result in a different effect of the pulses by the offset of the filter, as later pulses will have a larger portion of their area outside the time interval before the measurement. The values of the final measurements are then used to obtain an estimate of the parameters that characterize the filter.

As the Hamiltonian we are considering is proportional to $\sigma_x$, we will normalize all the Gaussian pulses to produce $\pi$ pulses, prepare the qubit initially in the state $|0 \rangle$, and measure the population in state $| 1 \rangle$. If no filter was present, the initial state would be flipped by the $\pi$ pulse and the population in state $|1 \rangle$ would be 1. Any deviations from this behavior represent the effect of the filter $\mathcal{L}$.

We will assume that the measurement results are given by some physical parameter (given in arbitrary units) linearly proportional to the population in the state $|1 \rangle$, but not necessarily the same, because the populations might not be directly accessible by the measurement apparatus. This results in two additional parameters that connect the measured results $M$ to the population $P_1$ of state $|1 \rangle$:

$$ M = a P_1 + b. $$Together, $\sigma$, $\mu$, $a$, $b$ form the four variables that our optimization engine will attempt to determine.

```
# Define parameters to be estimated
cutoff_frequency = 2.0*np.pi * 480.0e3 # Hz
filter_offset = 100e-9 # s
actual_a = 1.7
actual_b = -0.5
```

```
total_duration = 1000e-9 # s
# Define Gaussian pulse tensors, normalized to generate pi-pulses
gauss_pulse = lambda t, mu, width: tf.exp(-0.5*((t-mu)/width)**2.0) * np.sqrt(0.5*np.pi/width**2.0)
max_width = total_duration/6.0
mu_values = np.linspace(max_width, total_duration-max_width, 6)
width_values = np.linspace(max_width/8.0, max_width/2.0, 4)
param_values = np.reshape(np.transpose(np.meshgrid(width_values, mu_values)), (-1, 2))
# Define sampled times
segments = 150
t_values = np.linspace(0.0, total_duration, segments+1)
```

```
def population(width, mean, cutoff, offset):
"""
Calculates population of |1> after a Gaussian of pulse with a given
width and mean is applied to state |0>, while applying a Gaussian filter
with a given cutoff frequency and offset.
"""
gauss_filter = qc.create_gaussian_integral_function(
std=1.0/cutoff,
offset=offset,
)
with tf.Graph().as_default():
alpha = qc.create_pwc_signal(
values=gauss_pulse(t_values, mean, width),
duration=total_duration,
)
alpha_filtered = qc.convolve_pwc(alpha, gauss_filter)
alpha_discretized = qc.discretize_stf(
stf=alpha_filtered,
duration=total_duration,
segments_count=segments)
with tf.compat.v1.Session() as session:
filtered_pulse = session.run(alpha_discretized.values)
simulation_result = qc.calculate_unitary_evolution(
duration=total_duration,
shifts=[{
'operator': qc.HermitianOperator(0.5*sigma_x),
'control': [qc.RealSegment(total_duration/segments, v) for v in filtered_pulse],
}]
)
unitary = simulation_result[-1]['evolution_operator'].operator
final_state = unitary @ np.array([[1.0+0.0j], [0.0+0.0j]])
population_1 = np.abs(final_state[1][0])**2.0
return np.real(population_1)
```

```
populations = [population(width=width, mean=mean, cutoff=cutoff_frequency, offset=filter_offset)
for width, mean in param_values]
```

```
# actual measurement results will have an uncertainty associated to them
# we will estimate the standard deviation of this error as 1% of the population
population_error = 0.01 * np.ones_like(populations)
measurement_results = actual_a*np.random.normal(loc=populations,
scale=population_error) + actual_b
# rescaling error to arbitrary units
measurement_error = np.abs(actual_a)* population_error
```

```
# plot inputs
nb_graphs = len(width_values)
x_values = mu_values*1e6
y_values = np.reshape(measurement_results, [nb_graphs, -1])
y_errorbar = 2.0*np.reshape(measurement_error, [nb_graphs, -1])
for n in range(nb_graphs):
fig = plt.figure()
plt.title("Pulses with width " + str(np.around(width_values[n]*1e6, 3)) + " $\\mu s$")
plt.xlabel("Pulse center ($\\mu s$)")
plt.ylabel("Measured values, M (a.u.)")
plt.errorbar(x_values,
y_values[n],
yerr=y_errorbar[n],
fmt="s")
```

Each of these graphs represents the average measurement results of a series of simulated experiments where the offset of the Gaussian $\pi$ pulses vary. In each graph the width of the Gaussian control pulse is different, as indicated in their titles. The measured values $M$ plotted above are linearly proportional to the populations $P_1$:

$$ M = a P_1 + b. $$The measured values are assumed to be subject to errors with Gaussian probability distributions that have a standard deviation that is 1% of $a$. This means that we are assuming that the error in the value of the of the population is of about 1%. The error bars in the graphs represent two times the standard deviation, having therefore a reliability of 95%.

### Estimating the parameters of the filter

The next step involves performing an optimization in which the four unknown parameters are estimated. This is done by the optimizer by trying to find the form of the evolution that would most closely match the measured points. Given a set of input points $M_i$ and the corresponding values $m_i$ that would be expected from a chosen set of the filter parameters $\sigma, \mu, a, b$, the cost function is:

$$ C = \sum_i \frac{[M_i-m_i(\sigma, \mu, a, b)]^2}{2(\Delta M_i)^2}, $$where $\Delta M_i$ is the standard deviation of each of the measured input points $M_i$. Minimizing this cost function is equivalent to minimizing the negative log likelihood that a certain set of parameters $\sigma, \mu, a, b$ could generate the points $M_i$, for the case where the probability distribution of the errors is a Gaussian. In this case, the probability of a certain curve $m_i(\sigma, \mu, a, b)$ generating the points $M_i$ is the product of all the individual Gaussian probabilities:

$$ P = \prod_i \frac{1}{\sqrt{2\pi (\Delta M_i)^2}} \exp \left\{ - \frac{[M_i - m_i(\sigma, \mu, a, b)]^2}{2 (\Delta M_i)^2} \right\}. $$It is easy to see that the negative logarithm of $P$ is the cost $C$ plus constants.

Minimizing this cost gives us the best choice of parameters that generated the original dynamics of the system, and also allows us to calculate the precision of the estimated parameters. This is done by using the CramÃ©r-Rao bound to identify the Hessian of the cost function with the inverse of the covariance matrix for the variables estimated.

```
def create_optimization_specification(variable_factory):
# parameters to estimate
mean, sigma = tf.unstack(variable_factory.create_unbounded(
2,
initial_lower_bound=0.0,
initial_upper_bound=total_duration,
))
a, b = tf.unstack(variable_factory.create_unbounded(
2,
initial_lower_bound=-1.0,
initial_upper_bound=1.0,
))
params = tf.stack([mean, sigma, a, b])
mean, sigma, a, b = tf.unstack(params)
gauss_filter = qc.create_gaussian_integral_function(
std=sigma,
offset=mean,
)
populations = []
initial_state = np.array([[1.0+0.0j], [0.0+0.0j]])
for width, mu in param_values:
alpha = qc.create_pwc_signal(
values=gauss_pulse(t_values, mu, width),
duration=total_duration,
)
alpha_filtered = qc.convolve_pwc(alpha, gauss_filter)
shift = qc.create_stf_operator(
operator=qc.HermitianOperator(0.5*sigma_x),
signal=alpha_filtered,
)
simulation_result = qc.solve_rk4_function(
sample_times=t_values,
function=qc.map_stf(lambda values: -1.0j * values, shift),
)
unitary = qc.chain_matmul_sequential(simulation_result, reverse=True)
final_state = tf.matmul(unitary, initial_state)
population_1 = - tf.abs(tf.matmul(final_state, initial_state, adjoint_a=True))**2.0 + 1.0
populations += [a*population_1+b]
calculated_points = tf.reshape(tf.stack(populations), [-1])
# calculate cost
cost = tf.reduce_sum((calculated_points - measurement_results)**2.0/(2.0*measurement_error**2.0))
# calculate covariance matrix from the Hessian
hessian_matrix = tf.hessians(cost, params)[0]
covariance_matrix = tf.linalg.inv(hessian_matrix)
return qc.OptimizationSpecification(
cost=cost,
output_callable=lambda session: {"cost": session.run(cost),
"mean": session.run(mean),
"sigma": session.run(sigma),
"a": session.run(a),
"b": session.run(b),
"covariances": session.run(covariance_matrix),
},
)
```

```
result = qc.optimize(create_optimization_specification, number_of_optimizations=10)
cost = result["output"]["cost"]
calculated_mean = result["output"]["mean"]
calculated_sigma = result["output"]["sigma"]
calculated_a = result["output"]["a"]
calculated_b = result["output"]["b"]
# error bars calculated with 95% precision, having 2 sigma
uncertainty_mean = 2.0*np.sqrt(result["output"]["covariances"][0][0])
uncertainty_sigma = 2.0*np.sqrt(result["output"]["covariances"][1][1])
uncertainty_a = 2.0*np.sqrt(result["output"]["covariances"][2][2])
uncertainty_b = 2.0*np.sqrt(result["output"]["covariances"][3][3])
```

```
calculated_curves = []
mu_range = np.linspace(0.1*total_duration, 0.9*total_duration, 30)
for width in width_values:
calculated_curves += [calculated_a*population(width=width,
mean=mu,
cutoff=1.0/calculated_sigma,
offset=calculated_mean)
+calculated_b
for mu in mu_range]
```

```
ideal_curves = []
for width in width_values:
ideal_curves += [actual_a*population(width=width,
mean=mu,
cutoff=cutoff_frequency,
offset=filter_offset)
+actual_b
for mu in mu_range]
```

```
# plot results
nb_graphs = len(width_values)
y_values = np.reshape(measurement_results, [nb_graphs, -1])
y_errorbar = 2.0*np.reshape(measurement_error, [nb_graphs, -1])
y_actual = np.reshape(ideal_curves, [nb_graphs, -1])
y_estimates = np.reshape(calculated_curves, [nb_graphs, -1])
for n in range(nb_graphs):
fig = plt.figure()
plt.title("Pulses with width " + str(np.around(width_values[n]*1e6, 3)) + " $\\mu s$")
plt.xlabel("Pulse center ($\\mu s$)")
plt.ylabel("Measurement values, M (a.u.)")
plt.errorbar(mu_values*1e6,
y_values[n],
yerr=y_errorbar[n],
fmt="s",
color="C0",
label="Measured values")
plt.plot(mu_range*1e6, y_actual[n], color='C0', label="Ideal values", alpha=0.3)
plt.plot(mu_range*1e6, y_estimates[n], '--', label="Estimated values", color='C1')
plt.legend(loc=3)
# print parameter estimates
print ("real sigma =", 1.0/cutoff_frequency/1e-6, "us")
print ("estimated sigma = (", np.around(calculated_sigma/1e-6, 3),
"+-", np.around(uncertainty_sigma/1e-6, 3), ") us")
print ("real mean =", filter_offset/1e-6, "us")
print ("estimated mean = (", np.around(calculated_mean/1e-6, 3),
"+-", np.around(uncertainty_mean/1e-6, 3), ") us")
print ("actual a =", actual_a, "au")
print ("estimated a = (", np.around(calculated_a, 3),
"+-", np.around(uncertainty_a, 3), ") au")
print ("actual b =", actual_b, "au")
print ("estimated b = (", np.around(calculated_b, 3),
"+-", np.around(uncertainty_b, 3), ") au")
```

In the graphs above, we can compare the curves that represent the values that would be measured if the system followed the estimated filter (dashed), and the values that would be measured if perfect measurements were performed on a system subject to the real filter (solid). Each graph corresponds to sets of experiments where pulses of the same width are applied at different points in time, as indicated in their titles and $x$ axis. The close match between the two curves highlights the fact that the estimated parameters are close to the real values.

Above we also provide explicit estimates of the parameters $\sigma$ (`sigma`

) and $\mu$ (`mu`

) of the Gaussian kernel of the linear filter:

The other two parameters `a`

and `b`

represent the relationship between the values actually measured $M$ and the populations $P_1$:

The errors estimated for the variables above have a certainty of 2 times the standard deviation, or 95%.

## Identifying constant terms of a single-qubit Hamiltonian

The flexible optimizer also allows a user to characterize terms of a Hamiltonian to which their qubits are exposed. In this section, we will consider how to estimate a constant Hamiltonian that rotates a single qubit. In this case, preparing different initial states and letting the qubit evolve for different times before measuring it will reveal information about the parameters of the qubit rotation. This setup is illustrated below.

In more detail, if we exclude terms proportional to the identity (which affect all states in the same way and therefore are irrelevant to the qubit's evolution), a constant single-qubit Hamiltonian can be written in terms of the Pauli matrices:

$$ H = \frac{1}{2} \left( \Omega_x \sigma_x + \Omega_y \sigma_y + \Omega_z \sigma_z \right). $$The task of the optimizer in this case will be to determine these three parameters $\Omega_x$, $\Omega_y$, and $\Omega_z$. This can be done by preparing the system in different initial states, and then letting it evolve for different wait times, before measuring it. The sets of measured points are then provided to the optimization engine, which finds the parameters most likely to have generated that series of points.

### Simulating the input data

The experiments used here to determine a constant term in the Hamiltonian consist of simply letting a qubit evolve for different initial states, during different wait periods. We prepare these qubits in the eigenstates of $\sigma_z$, $\sigma_x$, and $\sigma_y$, and then measure the populations in eigenstates of $\sigma_y$, $\sigma_z$, and $\sigma_x$, respectively. Information about the evolution of these three initial states is sufficient to characterize the three parameters that form the Hamiltonian. Moreover, measuring populations that are not in the same eigenbasis as the initial states allows us to get information about the direction of the rotation of the qubit, depending on whether the population initially increases or decreases.

```
# Define parameters to be estimated
actual_Omega_x = (2.0*np.pi) * 0.5e6 # Hz
actual_Omega_y = (2.0*np.pi) * 1.5e6 # Hz
actual_Omega_z = (2.0*np.pi) * 1.8e6 # Hz
actual_Omegas = [actual_Omega_x, actual_Omega_y, actual_Omega_z]
```

```
# range of wait times in the different experiments
wait_times = np.linspace(25e-9, 500e-9, 20)
# list of initial states
initial_states = np.array([[[1.0+0.0j], [0.0+0.0j]],
[[1.0+0.0j], [1.0+0.0j]]/np.sqrt(2.0),
[[1.0+0.0j], [0.0+1.0j]]/np.sqrt(2.0),])
initial_state_names = [r"$\vert 0 \rangle$",
r"$(\vert 0 \rangle + \vert 1 \rangle)/\sqrt{2}$",
r"$(\vert 0 \rangle +i \vert 1 \rangle)/\sqrt{2}$"]
# list of states whose population will be measured at the end
projector_states = np.array([[[1.0+0.0j], [0.0+1.0j]]/np.sqrt(2.0),
[[1.0+0.0j], [0.0+0.0j]],
[[1.0+0.0j], [1.0+0.0j]]/np.sqrt(2.0),])
projector_state_names = [r"$(\vert 0 \rangle +i \vert 1 \rangle)/\sqrt{2}$",
r"$\vert 0 \rangle$",
r"$(\vert 0 \rangle + \vert 1 \rangle)/\sqrt{2}$"]
# creates lists of parameters for all the experiments
all_wait_times = np.tile(wait_times, len(initial_states)).flatten()
all_initial_states = np.repeat(initial_states, len(wait_times), axis=0)
all_projector_states = np.repeat(projector_states, len(wait_times), axis=0)
```

```
# define Hamiltonian
Hamiltonian = lambda Omegas: 0.5*(Omegas[0]*sigma_x+Omegas[1]*sigma_y+Omegas[2]*sigma_z)
# define unitaries for Hamiltonian
Ut = lambda t, Omegas: expm(-1.0j*t*Hamiltonian(Omegas))
# define populations
expectation_value = lambda wait_time, initial_state, projector_state, Omegas: np.abs(projector_state.T.conj()
@ Ut(wait_time, Omegas)
@ initial_state).item()**2.0
```

```
expectation_values = [expectation_value(wait_time=wait_time,
initial_state=initial_state,
projector_state=projector_state,
Omegas=actual_Omegas)
for wait_time, initial_state, projector_state in zip(all_wait_times,
all_initial_states,
all_projector_states,)]
```

```
# actual measurement results will have a certain uncertainty to them
# we will estimate the standard deviation of this error at 0.01
y_error = 0.01 * np.ones_like(expectation_values)
y_values = np.random.normal(loc=expectation_values, scale=y_error)
```

```
# plot inputs
x_values = wait_times*1e6
nb_graphs = len(initial_states)
y_flat_values = np.reshape(y_values, [nb_graphs, -1])
y_flat_errorbar = 2.0*np.reshape(y_error, [nb_graphs, -1])
for n in range(nb_graphs):
fig = plt.figure()
plt.title("Initial state = " + initial_state_names[n])
plt.ylabel("Population in state " + projector_state_names[n])
plt.xlabel("Wait time ($\\mu s$)")
plt.errorbar(x_values,
y_flat_values[n],
yerr=y_flat_errorbar[n],
fmt="s",
color='C0')
```

The graphs above show the values of the measurement results, considering an error whose standard deviation is 1% of the population. Each graph represent a set of experiments with the same initial state (given in the title) and same measured population (given in the $y$ axis). The error bars in the graphs have a length of $2\sigma$, or 95% reliability.

### Estimating the parameters of the Hamiltonian

The parameter estimation for the Hamiltonian terms is then performed by finding the form of the evolution that most closely matches the measured points. Given a set of input points $P_i$ and the corresponding populations $p_i$ that would be expected for a given set of Hamiltonian parameters $\Omega_x$, $\Omega_y$, and $\Omega_z$, the cost function is:

$$ C = \sum_i \frac{[P_i-p_i(\Omega_x, \Omega_y, \Omega_z)]^2}{2(\Delta P_i)^2}, $$where $\Delta P_i$ is the standard deviation of each of the measured points $P_i$. Minimizing this cost function is equivalent to minimizing the negative log likelihood that a certain set of parameters $\Omega_x, \Omega_y, \Omega_z$ could generate the points $P_i$, for the case where the probability distribution of the errors is a Gaussian. In this case, the probability of a certain curve $p_i(\Omega_x, \Omega_y, \Omega_z)$ generating the points $P_i$ is the product of all the individual Gaussian probabilities:

$$ P = \prod_i \frac{1}{\sqrt{2\pi (\Delta P_i)^2}} \exp \left\{ - \frac{[P_i - p_i(\Omega_x, \Omega_y, \Omega_z)]^2}{2 (\Delta P_i)^2} \right\}. $$It is easy to see that the negative logarithm of $P$ is the cost $C$ plus constants.

Minimizing the cost gives us the best choice of parameters that generate the original dynamics of the system, and also allows us to calculate the precision of the estimated parameters. This is done by using the CramÃ©r-Rao bound to identify the Hessian of the cost function with the inverse of the covariance matrix for the variables estimated.

```
def optimization_specification(variable_factory):
# forbids frequencies whose half-periods are shorter than the smaller spacing between points
min_time_resolution = np.min(np.abs(np.diff(np.unique(np.pad(all_wait_times, (1,0))))))
# create the three parameters
Omegas = variable_factory.create_bounded(
3,
lower_bound=-1.0/min_time_resolution,
upper_bound=1.0/min_time_resolution,
)
Omega_x, Omega_y, Omega_z = tf.unstack(Omegas)
# converts input values into constant tensors
wait_times_tensor = tf.constant(all_wait_times)
initial_states_tensor = tf.constant(all_initial_states)
projector_states_tensor = tf.constant(all_projector_states)
# calculates time evolution of the system
sigma_x_tiled = np.tile([sigma_x], [len(all_wait_times), 1, 1])
sigma_y_tiled = np.tile([sigma_y], [len(all_wait_times), 1, 1])
sigma_z_tiled = np.tile([sigma_z], [len(all_wait_times), 1, 1])
hamiltonian = 0.5*(tf.cast(Omega_x, dtype=tf.complex128) * tf.constant(sigma_x_tiled)
+ tf.cast(Omega_y, dtype=tf.complex128) * tf.constant(sigma_y_tiled)
+ tf.cast(Omega_z, dtype=tf.complex128) * tf.constant(sigma_z_tiled))
unitaries = qc.complex_expm_hermitian(
tf.cast(-wait_times_tensor[:, None, None], dtype=tf.complex128)
* hamiltonian
)
final_states_tensor = tf.matmul(unitaries, initial_states_tensor)
# calculate final population in the projector state
inner_products = tf.matmul(projector_states_tensor,
final_states_tensor,
adjoint_a=True)
populations = tf.abs(inner_products)**2.0
calculated_points = tf.reshape(populations, [-1])
# calculate cost
cost = tf.reduce_sum((calculated_points - y_values)**2.0/(2.0*y_error**2.0))
# covariance matrix is the inverse of the Hessian
Hessian_matrix = tf.hessians(cost, Omegas)[0]
covariance_matrix = tf.linalg.inv(Hessian_matrix)
return qc.OptimizationSpecification(
cost=cost,
output_callable=lambda session: {"cost": session.run(cost),
"Omegas": session.run(Omegas),
"covariances": session.run(covariance_matrix),
},
)
```

```
result = qc.optimize(optimization_specification, number_of_optimizations=30)
estimated_Omegas = result["output"]["Omegas"]
# choosing uncertainty as 2 sigma
Omega_x_uncertainty = 2.0*np.sqrt(result["output"]["covariances"][0][0])
Omega_y_uncertainty = 2.0*np.sqrt(result["output"]["covariances"][1][1])
Omega_z_uncertainty = 2.0*np.sqrt(result["output"]["covariances"][2][2])
```

```
# range of wait times in the different experiments
extra_wait_times = np.linspace(1e-9, 500e-9, 100)
# creates lists of parameters for all the points
more_wait_times = np.tile(extra_wait_times, len(initial_states)).flatten()
more_initial_states = np.repeat(initial_states, len(extra_wait_times), axis=0)
more_projector_states = np.repeat(projector_states, len(extra_wait_times), axis=0)
```

```
estimated_values = [expectation_value(wait_time=wait_time,
initial_state=initial_state,
projector_state=projector_state,
Omegas=estimated_Omegas)
for wait_time, initial_state, projector_state
in zip(more_wait_times, more_initial_states, more_projector_states)]
```

```
actual_values = [expectation_value(wait_time=wait_time,
initial_state=initial_state,
projector_state=projector_state,
Omegas=actual_Omegas)
for wait_time, initial_state, projector_state
in zip(more_wait_times, more_initial_states, more_projector_states)]
```

```
# plot results
more_x_values = extra_wait_times/1e-6
y_estimates = np.reshape(estimated_values, [nb_graphs, -1])
y_actual = np.reshape(actual_values, [nb_graphs, -1])
for n in range(nb_graphs):
fig = plt.figure()
plt.title("Initial state = " + initial_state_names[n])
plt.ylabel("Population in state " + projector_state_names[n])
plt.xlabel("Wait time ($\\mu s$)")
plt.errorbar(x_values,
y_flat_values[n],
yerr=y_flat_errorbar[n],
fmt="s",
label="Measured populations",
color='C0')
plt.plot(more_x_values, y_actual[n], color='C0', label="Ideal populations", alpha=0.3)
plt.plot(more_x_values, y_estimates[n], '--', label="Estimated populations", color='C1')
plt.legend(loc=2, bbox_to_anchor=(1, 1))
# print parameter estimates
print ("real Omega_x =", actual_Omegas[0] * 1e-6, "MHz")
print ("estimated Omega_x = (", np.around(estimated_Omegas[0] * 1e-6, 3),
"+-", np.around(Omega_x_uncertainty* 1e-6, 3), ") MHz")
print ("real Omega_y =", actual_Omegas[1] * 1e-6, "MHz")
print ("estimated Omega_y = (", np.around(estimated_Omegas[1] * 1e-6, 3),
"+-", np.around(Omega_y_uncertainty* 1e-6, 3), ") MHz")
print ("real Omega_z =", actual_Omegas[2] * 1e-6, "MHz")
print ("estimated Omega_z = (", np.around(estimated_Omegas[2] * 1e-6, 3),
"+-", np.around(Omega_z_uncertainty* 1e-6, 3), ") MHz")
```

The three parameters estimated above, $\Omega_x$ (`Omega_x`

), $\Omega_y$ (`Omega_y`

), and $\Omega_z$ (`Omega_z`

) correspond to the three coefficients present in the single-qubit Hamiltonian given in the beginning of this section:

The errors in the parameter estimates correspond to 2 times the standard deviation and thus have 95% reliability.

In the graphs above, we can compare the curves that represent the populations that would be measured if the system followed the estimated Hamiltonian (dashed), and the populations that would be measured if perfect measurements were performed on the real Hamiltonian (solid). Each graph corresponds to sets of experiments with a given initial state and final measurement, as indicated in the title and $y$ label of each of them. The close match between the two curves highlights the fact that the estimated parameters are close to the real values.