# How to use QuTiP operators in graphs

**Incorporate QuTiP objects and programming syntax directly into graphs**

It is possible to use the Q-CTRL Python package with QuTiP operators. This is possible in graph execution, where every operator can be replaced by QuTiP's Qobj. In case QuTiP is not available, these operators are NumPy arrays.

## Worked example: QuTiP operators in a graph-based simulation of a single qubit with leakage

In this example we to simulate the dynamics of a single qubit subject to leakage to an additional level. The resulting qutrit is treated as an oscillator (truncated to three levels) with an anharmonicity of $\chi$, described by the Hamiltonian:

\begin{align*} H(t) = & \frac{\chi}{2} (a^\dagger)^2 a^2 + \frac{\Omega(t)}{2} a + \frac{\Omega^*(t)}{2} a^\dagger, \end{align*}where $a = \left|0 \right\rangle \left\langle 1 \right| + \sqrt{2} \left|1 \right\rangle \left\langle 2 \right|$ is the lowering operator and $\Omega(t)$ is a time-dependent Rabi rate.

In the example code below we illustrate how the noise-free infidelity $\mathcal{I}_0$ can be extracted from the coherent simulation results. To do this, we choose as target an X gate between the states $\left| 0 \right\rangle$ and $\left| 1 \right\rangle$. Notice that this target is not unitary in the total Hilbert space, but is still a valid target because it is a partial isometry—in other words, it is unitary in the subspace whose basis is $\{ \left| 0 \right\rangle, \left| 1 \right\rangle \}$.

Note the use of QuTiP to define operators where convenient in defining the relevant graph nodes.

```
import matplotlib.pyplot as plt
import numpy as np
from qctrlopencontrols import new_primitive_control
from qctrlvisualizer import get_qctrl_style
from qctrl import Qctrl
# Starting a session with the API
qctrl = Qctrl()
plt.style.use(get_qctrl_style())
```

```
# Define system parameters
chi = 2 * np.pi * 3 * 1e6 # Hz
omega_max = 2 * np.pi * 1e6 # Hz
total_rotation = np.pi
# Define pulse using pulses from Q-CTRL Open Controls
pulse = new_primitive_control(
rabi_rotation=total_rotation,
azimuthal_angle=0.0,
maximum_rabi_rate=omega_max,
name="primitive",
)
sample_times = np.linspace(0, pulse.duration, 100)
initial_state = np.array([1.0, 0.0, 0.0])
try:
import qutip as qt
# Define matrices for the Hamiltonian operators
a = qt.destroy(3)
ad = qt.create(3)
ad2a2 = ad ** 2 * a ** 2
# Define the target
target_operation = qt.Qobj(
np.array([[0, 1, 0], [1, 0, 0], [0, 0, 0]], dtype=complex)
)
except ModuleNotFoundError:
# Define matrices for the Hamiltonian operators
a = np.diag([1, np.sqrt(2)], 1)
ad = np.diag([1, np.sqrt(2)], -1)
ad2a2 = ad @ ad @ a @ a
# Define the target
target_operation = np.array([[0, 1, 0], [1, 0, 0], [0, 0, 0]], dtype=complex)
pass
# in the graph nodes below we define operators in the Hamiltonian using qutip definitions
graph = qctrl.create_graph()
# Define the anharmonic term
anharmonic_drift = graph.constant_pwc_operator(
duration=pulse.duration, operator=0.5 * chi * ad2a2
)
# Define Rabi drive term
rabi_drive = graph.pwc_operator_hermitian_part(
graph.pwc_operator(
signal=graph.pwc(
durations=pulse.durations,
values=pulse.rabi_rates * np.exp(1j * pulse.azimuthal_angles),
),
operator=a,
)
)
hamiltonian = anharmonic_drift + rabi_drive
# Calculate the time evolution operators
unitaries = graph.time_evolution_operators_pwc(
hamiltonian=hamiltonian, sample_times=sample_times, name="time_evolution_operators"
)
# Calculate the infidelity
infidelity = graph.infidelity_pwc(
hamiltonian=hamiltonian, target=graph.target(target_operation), name="infidelity"
)
evolved_states = unitaries @ initial_state[:, None]
evolved_states.name = "evolved_states"
# Run simulation
graph_result = qctrl.functions.calculate_graph(
graph=graph,
output_node_names=["evolved_states", "infidelity", "time_evolution_operators"],
)
# Extract and print final infidelity
print(
"Noise-free infidelity at end: {}".format(
graph_result.output["infidelity"]["value"]
)
)
# Extract and print final time evolution operator
print("Time evolution operator at end:")
print(graph_result.output["evolved_states"]["value"][-1])
# Extract and plot state populations
state_vectors = graph_result.output["evolved_states"]["value"]
for state in range(3):
plt.plot(
sample_times * 1e6, np.abs(state_vectors[:, state]) ** 2, label=f"$P_{state}$"
)
plt.xlabel("Time (µs)")
plt.ylabel("Probability")
plt.title("State populations")
plt.legend()
plt.show()
```