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

\[\hat H(t)=\hat H_0 + \delta \hat V \cos (\omega t)\]

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

\[\hat H(t)=\hat H_0 + \delta \hat V[1+\beta(t)]\cos (\omega t)\]

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.

Expectation values of Pauli operators ignoring charge noise

Fig. 8.4.4 Expectation values of Pauli operators ignoring charge noise

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.

Lorentzian spectral density for the charge noise

Fig. 8.5.1 Lorentzian spectral density for the charge noise

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.

A sample charge noise

Fig. 8.5.2 A sample charge noise

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.

Dynamics with and without charge noise

Fig. 8.5.3 Dynamics with and without charge noise

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()