# How to simulate quantum dynamics subject to noise with graphs

**Simulate the dynamics of closed quantum systems in the presence of Non-Markovian noise**

The Q-CTRL Python package enables you to simulate the dynamics of quantum systems that are affected by various noise processes. Such simulations can provide useful insights into the expected real-world performance of candidate control solutions. In this notebook we show how to use the Q-CTRL Python package to perform simulations of systems affected by different types of noise using graphs.

**Note** This approach is preferable when a simulation becomes part of a larger computation where flexibility is prioritized. For instance, simulation of noisy quantum-system dynamics in the presence of multiple control signals, each distorted by a transmission line.

## Summary workflow: Non-Markovian noise simulation using graphs

The simulation of quantum dynamics subject to Non-Markovian (colored) noise is used to provide insights into the time-domain evolution of a noise process characterized by a power spectral density with strong low-frequency content.

### 1. Create noise power spectrum

An arbitrary frequency-domain representation of a noise power spectrum may simply be created using Python commands.

### 2. Link sampled noise trajectories to a graph node

A simulation graph can be built as explained in the How to simulate quantum dynamics for noiseless systems using graphs user guide. To include noise in the simulation you only need to add a term corresponding to the noise Hamiltonian to this graph.

Given the power-spectral density defined in step 1, you can create random noise trajectories in the time domain with average characteristics reflected by the power spectrum. You can generate batches of noise trajectories using the `graph.random_colored_noise_stf_signal`

graph operation. This node creates a batch of continuous noise trajectories from the power density of the pink noise. These trajectories can then be discretized with a user-defined number of equally spaced points between the start and the end of the operation evolution.

From these batches of noise, you can create a batch of noise Hamiltonians where each Hamiltonian of the batch represents a different realization of the noise variable. The batch of noise Hamiltonians is then added to the (noise-free) system Hamiltonian to create a batch of total Hamiltonians. The Q-CTRL Python package automatically reconciles differing time segmentation while summing over Hamiltonian terms.

Next create a unitary time-evolution operator from the total Hamiltonian via `graph.time_evolution_operators_pwc`

and apply it to evolve the initial state into a final state.

### 3. Create and execute graph

Once the overall graph is created it may be executed as usual via `qctrl.functions.calculate_graph`

.

## Worked example: Non-Markovian-noise simulation of a single qubit subject to dephasing noise

In this example we will show how to simulate the dynamics of a single qubit experiencing $1/f$ dephasing noise, driven by a $\pi/2$ CORPSE pulse. The Hamiltonian of the quantum system is:

\begin{align*} H(t) = & \frac{\Omega(t)}{2} \sigma_- + \frac{\Omega^*(t)}{2} \sigma_+ + \frac{\eta(t)}{2} \sigma_z, \end{align*}where $\Omega(t)$ is a time-dependent Rabi rate, $\eta(t)$ is a stochastic dephasing noise process, $\sigma_\pm = (\sigma_x \mp i \sigma_y)/2$, and $\sigma_k$ are the Pauli matrices.

We assume that $\eta(t)$ represents pink noise ($1/f$ noise), as defined below.

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

```
# Define the non-Markovian noise power spectrum
def pink_spectrum(frequencies, frequency_cutoff):
return 1 / (frequencies + frequency_cutoff)
frequency_step = 2 * 1e3
frequencies = np.arange(0, 2 * 1e6, frequency_step)
power_densities = 4e9 * pink_spectrum(
frequencies=frequencies, frequency_cutoff=0.05 * 1e6
)
fig = plt.figure(figsize=(10, 5))
fig.suptitle("Dephasing noise spectral density")
plt.plot(frequencies / 1e6, power_densities * 1e6, color="#EB64B6")
plt.fill_between(frequencies / 1e6, 0, power_densities * 1e6, color="#F790DB")
plt.xlabel("Frequency (MHz)")
plt.ylabel("Power density (1/MHz)")
plt.show()
```

### Representing the quantum system and noise operator

The drive Hamiltonian is defined as a $\pi/2$ CORPSE pulse, drawn directly from the Q-CTRL Open Controls library.

```
# Define standard matrices
identity = np.array([[1, 0], [0, 1]], dtype=complex)
sigma_x = np.array([[0, 1], [1, 0]], dtype=complex)
sigma_y = np.array([[0, -1j], [1j, 0]], dtype=complex)
sigma_z = np.array([[1, 0], [0, -1]], dtype=complex)
sigma_m = np.array([[0, 1], [0, 0]], dtype=complex)
# Define control parameters
omega_max = 2 * np.pi * 1e6 # Hz
total_rotation = np.pi / 2
# Import pulse from Q-CTRL Open Controls
pulse = new_corpse_control(
rabi_rotation=total_rotation,
azimuthal_angle=0.0,
maximum_rabi_rate=omega_max,
name="CORPSE",
)
```

In the simulation graph, we define batches of noise trajectories using the `graph.random_colored_noise_stf_signal`

node. This node creates a batch of continuous noise trajectories from the power density, which we then discretize to $50$ equally spaced points in the system evolution using `graph.discretize_stf`

. From these batches of noise, we will create a batch of noise Hamiltonians where each Hamiltonian of the batch represents a different realization of the noise variable.

We will then create the full Hamiltonian by adding the noise and drive terms together. The noise Hamiltonians were sampled on $50$ time points whereas the CORPSE Hamiltonian contains only $3$ time segments. The Q-CTRL Python package automatically reconciles differing time segmentation while summing over Hamiltonian terms.

```
# Run this number of simulations with different noise trajectories
simulation_count = 80
def get_computation_graph():
graph = qctrl.create_graph()
# Signal for Omega extracted from the imported pulse
omega_signal = graph.pwc(
values=pulse.rabi_rates * np.exp(1j * pulse.azimuthal_angles),
durations=pulse.durations,
name="omega",
)
# System Hamiltonian
system_hamiltonian = graph.pwc_operator_hermitian_part(omega_signal * sigma_m)
# Noise Hamiltonian. Generate random noise trajectories.
noise_stf = graph.random_colored_noise_stf_signal(
power_spectral_density=power_densities,
frequency_step=frequency_step,
batch_shape=(simulation_count,),
)
noise_pwc = graph.discretize_stf(
stf=noise_stf, duration=pulse.duration, segment_count=50
)
noise_hamiltonian = noise_pwc * sigma_z
# Combine the system and the noise Hamiltonian to get the total Hamiltonian.
# We can automatically combine a batched and a non-batched term.
hamiltonian = system_hamiltonian + noise_hamiltonian
# Compute the unitary using `time_evolution_operators_pwc`
unitaries = graph.time_evolution_operators_pwc(
hamiltonian=hamiltonian, sample_times=sample_times
)
# Use the unitary to evolve the initial state into the final state.
final_states = unitaries @ initial_state
final_states.name = "final_states"
return graph
```

```
# Define sample times for the output and initial state
sample_times = np.linspace(0, np.sum(pulse.durations), 100)
initial_state = np.array([[1.0], [0.0]])
# Run simulation
simulation_result = qctrl.functions.calculate_graph(
graph=get_computation_graph(), output_node_names=["final_states"]
)
final_states = simulation_result.output["final_states"]["value"]
```

```
def get_density_matrix(final_states):
# Density matrix from simulated states, ρ = |ψ><ψ|
density_matrix = final_states * final_states.conj().swapaxes(-1, -2)
return density_matrix
# Calculate and plot average sigma_x expectation
average_density_matrices = np.average(get_density_matrix(final_states), axis=0)
fig = plt.figure(figsize=(10, 5))
plt.suptitle(r"Average $\sigma_x$ expectation")
plt.plot(
sample_times * 1e6,
[
np.real(np.trace(sigma_x.dot(density_matrix)))
for density_matrix in average_density_matrices
],
)
plt.axhline(y=0.0, c="k", lw=0.5)
plt.xlabel("Time (µs)")
plt.ylabel(r"$\overline{\langle\sigma_x\rangle}$")
plt.show()
```

```
display_bloch_sphere_from_density_matrices(density_matrices=average_density_matrices)
```

```
# Calculate and plot individual sigma_x expectations
all_density_matrices = get_density_matrix(final_states)
fig = plt.figure(figsize=(10, 5))
fig.suptitle(r"$\sigma_x$ expectations")
for density_matrices in all_density_matrices[0:10]:
plt.plot(
sample_times * 1e6,
[
np.real(np.trace(sigma_x.dot(density_matrix)))
for density_matrix in density_matrices
],
)
plt.axhline(y=0, c="k", lw=0.5)
plt.xlabel("Time (µs)")
plt.ylabel(r"$\langle\sigma_x\rangle$")
plt.show()
```