# 8. Charge noise in quantum dots

## 8.1. Requirements

### 8.1.1. Software components

QTCAD

### 8.1.2. Mesh file

`qtcad/examples/practical_application/meshes/gated_dot.msh`

### 8.1.3. Python script

`qtcad/examples/practical_application/8-noise.py`

### 8.1.4. References

## 8.2. Briefing

In the previous tutorial, we have investigated the dynamics of a qubit subject to electric dipole spin resonance. We started out by considering Rabi oscillations within the subspace of the two lowest-energy eigenstates. We then went beyond this approximation by considering leakage into higher-energy states. In this tutorial, we will consider another nonideality to which a qubit may be subject to: charge noise.

Specifically, in Electric dipole spin resonance and Rabi oscillations, the Hamiltonian was given by

where \(\hat H_0\) is the qubit Hamiltonian and \(\delta \hat V\) is the Hamiltonian describing the amplitude of the perturbation produced by the time-dependent electric field (the frequency of which is given by \(\omega\)).

Instead, in this tutorial, we will aim to simulate the dynamics assuming the Hamiltonian

where \(\beta(t)\) is a stationary and ergodic Gaussian random process with vanishing mean—it represents charge noise on the gate driving the qubit. More theoretical background on noisy EDSR is given in the subsection called A more realistic scenario: noisy electric-dipole spin resonance in the Quantum control theory section.

## 8.3. System initialization

We start by importing relevant modules.

```
import numpy as np
import pathlib
from matplotlib import pyplot as plt
from qtcad.device import constants as ct
import qutip
from qtcad.qubit import spectra, dynamics, noise
```

All these modules were previously seen, except from `noise`

, which is used to
model two-level systems subject to noise.

The system and drive Hamiltonians are then loaded from the results of Electric dipole spin resonance and Rabi oscillations and defined in units where \(\hbar=1\) (as required by QuTiP).

```
# Paths to system and drive Hamiltonians
script_dir = pathlib.Path(__file__).parent.resolve()
path_h0 = str(script_dir/'output'/'h0.npy')
path_u = str(script_dir/'output'/'u.npy')
# Load system and drive Hamiltonians and rewrite in units of hbar = 1
H0 = np.load(path_h0)/ct.hbar
delta_V = np.load(path_u)/ct.hbar
```

Several timescales are relevant for the simulation that we will perform. These timescales include the inverse of the qubit frequency (i.e., the difference between the two spin-qubit eigenenergies) and the Rabi period.

```
# Define relevant frequencies and time scales of the problem
omega = H0[0,0] - H0[1,1] # qubit frequency
omega_Rabi = np.absolute(delta_V[0,1]) # Rabi frequency on resonance
T_Rabi = 2*np.pi/omega_Rabi # Rabi period
```

We also need to define a time array on which we will sample the system and noise dynamics.

```
# Define time array
T = 20 * T_Rabi # Total time of the simulations
times = np.linspace(0, T, 1000) # Times at which dynamics are simulated.
```

Finally, we initialize a `Dynamics`

object and a `Noise`

object, which will
be used later.

```
# Initialize objects that will handle dynamics and noise.
dyn = dynamics.Dynamics()
N = noise.Noise()
```

## 8.4. Simulation without noise

To model the dynamics of our system without noise, we use a procedure similar
as in the previous tutorial. To simplify the treatment, we will employ the
`H_RF`

method of the `dyn`

object, which returns a QuTiP Qobj corresponding
to the Hamiltonian \(\hat H(t)=\hat H_0 + \delta \hat V \cos (\omega t)\).

```
# Simple simulation without noise
# The Hamiltonian H = H0 + delta_V cos(omega*t) in a frame rotating
# at frequency omega.
H = dyn.H_RF(H0, delta_V, omega)
```

The system dynamics can then be solved as before.

```
# Solve for the dynamics without any noise.
psi0 = qutip.basis(2, 0) # initial state
result0 = qutip.mesolve(H, psi0, times)
```

Next, we compute and plot the expectation values of Pauli operators,
\(\braket{\sigma_z(t)}\) and \(\braket{\sigma_y(t)}\) using
the `qutip.expect`

function.

```
# Plot expectation values of sigma_z and sigma_y as a function of time.
fig, axes = plt.subplots(2)
fig.set_size_inches((8,8))
#<sigma_z (t)>
axes[0].plot(result0.times, qutip.expect(qutip.sigmaz(), result0.states))
#<sigma_y (t)>
axes[1].plot(result0.times, qutip.expect(qutip.sigmay(), result0.states))
# Labels
axes[0].set_ylabel(f'$\left< \sigma_z (t) \\right>$', fontsize=20)
axes[1].set_ylabel(f'$\left< \sigma_y (t) \\right>$', fontsize=20)
axes[0].set_xlabel(f'$t \, (s)$', fontsize=20)
axes[1].set_xlabel(f'$t \, (s)$', fontsize=20)
plt.show()
```

This should result in the following figure.

## 8.5. Simulation with noise

### 8.5.1. Modeling the charge noise

Following the procedure described in Electric dipole spin resonance—Noise, to generate a
stationary, Gaussian, and ergodic random process with vanishing mean
\(\beta(t)\), we start by defining a spectral density for the noise.
For this purpose, we take this spectral density to be a Lorentzian centered
at angular frequency `w0 = 0`

with width equal to the Rabi frequency,
which we have previously stored as `omega_Rabi`

. Furthermore, we define
a sampling frequency (or frequency resolution) for the discrete implementation
of this spectral density in such a way that its inverse, `T0`

is much greater
than all other timescales of the problem. Finally, we define the total
power of this noise, `S0`

. With these elements, we can define the Lorentzian
spectral density using the `spectra.lorentz`

method.

```
# Include charge noise which is assumed to affect only the gate which
# generates the drive delta_V.
# Lorentzian spectrum
S0 = 0.01 # Total power
wc = omega_Rabi # Cutoff frequency
w0 = 0 # Central frequency
# T0 should be much larger than the largest time scale of the problem under consideration.
# Used to properly define the Fourier series of the time series related to the noise.
T0 = 50 * T_Rabi # Must take T0 such that 2*pi/T0 << wc to resolve the spectrum
omega_max = 5*wc
wk = np.arange(0, omega_max , 2*np.pi/T0) + w0
spectrum = spectra.lorentz(wk, S0, wc, w0 = w0)
```

This Lorentzian can then be plotted as follows.

```
# Plot spectrum
fig, axes = plt.subplots(1)
fig.set_size_inches((10,6))
axes.plot(wk, spectrum)
axes.set_ylabel(f'$S_\\beta(\omega) \, (s)$', fontsize=20)
axes.set_xlabel('$\omega \, (s^{-1})$', fontsize=20)
axes.set_title('Lorentzian spectral density', fontsize=30)
plt.show()
```

This should result in the following figure.

Having obtained the spectral density, a time series for the charge noise can be
randomly generated using the `N.gen_process`

method.

```
# Plot time series for the charge noise affecting the gate responsible for
# delta_V.
beta = N.gen_process(times, spectrum, wk, T0)
fig, axes = plt.subplots(1)
fig.set_size_inches((8,5))
axes.plot(times, beta)
axes.set_ylabel(f'$\\beta(t)$', fontsize=20)
axes.set_xlabel('$t \, (s)$', fontsize=20)
plt.show()
```

Do to the random nature of this noise generation, simulation results will differ between runs. In our case, we obtained the following noise sample.

### 8.5.2. Dynamics with noise

To simulate the dynamics of our qubit assuming charge noise of the type
described in the previous section, we use the `dynamics`

method
of a QTCAD `Noise`

object

```
# Simulate dynamic including noise.
num_runs1 = 1 # number of runs to average over
sigz_1run = N.dynamics(H0, delta_V, omega, spectrum, T0, psi0, times,
qutip.sigmaz(), vec_omega=wk, num_runs=num_runs1)
```

This method takes several arguments. `H0`

and `delta_V`

specify the system and
drive Hamiltonians (just as for the noiseless simulations). The charge noise
spectrum is specified using the argument `spectrum`

. `T0`

is inversely
proportional to the frequency resolution of the discretized noise process
in our simulation (more precisely, `T0`

sets an artificial period of the
charge noise). `psi0`

is the initial state of the system, and
`times`

is an array containing the instants in time over which we perform
the simulation. The operator whose expectation value we wish to compute is the
Pauli \(\sigma_z\) operator, i.e. `qutip.sigmaz()`

. The angular frequency
components at which the charge noise spectrum is evaluated is specified as
`vec_omega=wk`

. Finally, the number of simulations (or runs) over which the
expectation value is averaged is specified as `num_runs=num_runs1`

. Above, we
perform a single run, but we may also average over a more statistically
significant number of runs, as shown below.

```
num_runs2 = 100 # number of runs to average over
sigz_multiruns = N.dynamics(H0, delta_V, omega, spectrum, T0, psi0, times,
qutip.sigmaz(), vec_omega=wk, num_runs=num_runs2)
```

We may then plot the results, comparing with noiseless simulations.

```
# Plot dynamics
fig, axes = plt.subplots(2)
axes[0].plot(result0.times, qutip.expect(qutip.sigmaz(), result0.states),
label="Without noise")
axes[0].plot(result0.times, sigz_1run, label="With noise")
axes[1].plot(result0.times, qutip.expect(qutip.sigmaz(), result0.states),
label="Without noise")
axes[1].plot(result0.times, sigz_multiruns, label="With noise")
# Format Plots
axes[0].set_ylabel(f'$\left< \sigma_z (t) \\right>$', fontsize=20)
axes[1].set_ylabel(f'$\left< \sigma_z (t) \\right>$', fontsize=20)
axes[0].set_xlabel(f'$t \, (s)$', fontsize=20)
axes[1].set_xlabel(f'$t \, (s)$', fontsize=20)
axes[0].set_title(f'num_runs = {num_runs1}')
axes[1].set_title(f'num_runs = {num_runs2}')
axes[0].legend(loc='center left', fontsize=16)
axes[1].legend(loc='center left', fontsize=16)
fig.tight_layout()
plt.show()
```

We obtain the following results.

While the single-run simulation results may display variability from simulation to simulation, the average over \(100\) runs exhibits a decay of the Rabi oscillations over time.

## 8.6. Full code

```
__copyright__ = "Copyright 2024, Nanoacademic Technologies Inc."
import numpy as np
import pathlib
from matplotlib import pyplot as plt
from qtcad.device import constants as ct
import qutip
from qtcad.qubit import spectra, dynamics, noise
# Paths to system and drive Hamiltonians
script_dir = pathlib.Path(__file__).parent.resolve()
path_h0 = str(script_dir/'output'/'h0.npy')
path_u = str(script_dir/'output'/'u.npy')
# Load system and drive Hamiltonians and rewrite in units of hbar = 1
H0 = np.load(path_h0)/ct.hbar
delta_V = np.load(path_u)/ct.hbar
# Define relevant frequencies and time scales of the problem
omega = H0[0,0] - H0[1,1] # qubit frequency
omega_Rabi = np.absolute(delta_V[0,1]) # Rabi frequency on resonance
T_Rabi = 2*np.pi/omega_Rabi # Rabi period
# Define time array
T = 20 * T_Rabi # Total time of the simulations
times = np.linspace(0, T, 1000) # Times at which dynamics are simulated.
# Initialize objects that will handle dynamics and noise.
dyn = dynamics.Dynamics()
N = noise.Noise()
# Simple simulation without noise
# The Hamiltonian H = H0 + delta_V cos(omega*t) in a frame rotating
# at frequency omega.
H = dyn.H_RF(H0, delta_V, omega)
# Solve for the dynamics without any noise.
psi0 = qutip.basis(2, 0) # initial state
result0 = qutip.mesolve(H, psi0, times)
# Plot expectation values of sigma_z and sigma_y as a function of time.
fig, axes = plt.subplots(2)
fig.set_size_inches((8,8))
#<sigma_z (t)>
axes[0].plot(result0.times, qutip.expect(qutip.sigmaz(), result0.states))
#<sigma_y (t)>
axes[1].plot(result0.times, qutip.expect(qutip.sigmay(), result0.states))
# Labels
axes[0].set_ylabel(f'$\left< \sigma_z (t) \\right>$', fontsize=20)
axes[1].set_ylabel(f'$\left< \sigma_y (t) \\right>$', fontsize=20)
axes[0].set_xlabel(f'$t \, (s)$', fontsize=20)
axes[1].set_xlabel(f'$t \, (s)$', fontsize=20)
plt.show()
# Include charge noise which is assumed to affect only the gate which
# generates the drive delta_V.
# Lorentzian spectrum
S0 = 0.01 # Total power
wc = omega_Rabi # Cutoff frequency
w0 = 0 # Central frequency
# T0 should be much larger than the largest time scale of the problem under consideration.
# Used to properly define the Fourier series of the time series related to the noise.
T0 = 50 * T_Rabi # Must take T0 such that 2*pi/T0 << wc to resolve the spectrum
omega_max = 5*wc
wk = np.arange(0, omega_max , 2*np.pi/T0) + w0
spectrum = spectra.lorentz(wk, S0, wc, w0 = w0)
# Plot spectrum
fig, axes = plt.subplots(1)
fig.set_size_inches((10,6))
axes.plot(wk, spectrum)
axes.set_ylabel(f'$S_\\beta(\omega) \, (s)$', fontsize=20)
axes.set_xlabel('$\omega \, (s^{-1})$', fontsize=20)
axes.set_title('Lorentzian spectral density', fontsize=30)
plt.show()
# Plot time series for the charge noise affecting the gate responsible for
# delta_V.
beta = N.gen_process(times, spectrum, wk, T0)
fig, axes = plt.subplots(1)
fig.set_size_inches((8,5))
axes.plot(times, beta)
axes.set_ylabel(f'$\\beta(t)$', fontsize=20)
axes.set_xlabel('$t \, (s)$', fontsize=20)
plt.show()
# Simulate dynamic including noise.
num_runs1 = 1 # number of runs to average over
sigz_1run = N.dynamics(H0, delta_V, omega, spectrum, T0, psi0, times,
qutip.sigmaz(), vec_omega=wk, num_runs=num_runs1)
num_runs2 = 100 # number of runs to average over
sigz_multiruns = N.dynamics(H0, delta_V, omega, spectrum, T0, psi0, times,
qutip.sigmaz(), vec_omega=wk, num_runs=num_runs2)
# Plot dynamics
fig, axes = plt.subplots(2)
axes[0].plot(result0.times, qutip.expect(qutip.sigmaz(), result0.states),
label="Without noise")
axes[0].plot(result0.times, sigz_1run, label="With noise")
axes[1].plot(result0.times, qutip.expect(qutip.sigmaz(), result0.states),
label="Without noise")
axes[1].plot(result0.times, sigz_multiruns, label="With noise")
# Format Plots
axes[0].set_ylabel(f'$\left< \sigma_z (t) \\right>$', fontsize=20)
axes[1].set_ylabel(f'$\left< \sigma_z (t) \\right>$', fontsize=20)
axes[0].set_xlabel(f'$t \, (s)$', fontsize=20)
axes[1].set_xlabel(f'$t \, (s)$', fontsize=20)
axes[0].set_title(f'num_runs = {num_runs1}')
axes[1].set_title(f'num_runs = {num_runs2}')
axes[0].legend(loc='center left', fontsize=16)
axes[1].legend(loc='center left', fontsize=16)
fig.tight_layout()
plt.show()
```