6. Electric dipole spin resonance and Rabi oscillations
6.1. Requirements
6.1.1. Software components
QTCAD
6.1.2. Mesh file
qtcad/examples/practical_application/meshes/gated_dot.msh2
6.1.3. Python script
qtcad/examples/practical_application/6-EDSR.py
6.1.4. References
6.2. Briefing
In the previous tutorial, we obtained the charge-stability diagram of our gated quantum dot. Using this diagram, we can extract the voltage configurations under which a single electron resides in the quantum dot. Furthermore, we have computed the energy levels of such an electron. Assuming that the second-excited eigenenergy is sufficiently far from the first-excited eigenenergy, our quantum dot’s electron constitutes a qubit. Typically, for a single quantum dot used as a Loss-DiVincenzo spin qubit, the first two levels correspond to an electron in its orbital ground state with Zeeman-split spin states, while the second and third excited states contain an excitation of the electron’s orbital degree of freedom. In this tutorial, we will simulate an electric dipole spin resonance (EDSR) experiment on this qubit, driving transitions between the ground state and first-excited state using a time-varying voltage applied to a gate and an inhomogeneous magnetic field.
6.3. Setting up the device
6.3.1. Header
We start by importing relevant modules.
import numpy as np
import pathlib
from matplotlib import pyplot as plt
import qutip
import pickle
from device import constants as ct
from device.mesh3d import Mesh, SubMesh
from device import io
from device import Device, SubDevice
from device.solver_params import SolverParams
from device.schrodinger import Solver
from device import materials as mt
from qubit.edsr import MicromagnetEDSR
from qubit import dynamics
We have previously seen all of these modules, except from edsr and
dynamics. The former is used to simulate EDSR experiments—specifically
to modulate the voltages applied to the gates of quantum dots and compute the
resulting potential energy matrices. The latter contains various calculators
and visualizers for system dynamics, and harnesses the
QuTiP python library.
6.3.2. Creating a quantum dot with a single electron
The next step is to define our device and set external voltages, as was done
similarly in previous tutorials. Based on the charge-stability diagram obtained
in Calculating the charge-stability diagram of a quantum dot, a plunger gate potential of \(0.4\textrm{ V}\)
results in a single electron residing in the dot. Thus, the applied potentials
are defined to reflect this voltage configuration; notably Vtop_2 and
Vbottom_2 are set to 0.4. Furthermore, the electric potential,
d.phi, is loaded from output/potential_nrg_Vp_0.4.hdf5
(recall that using the
set_potential method of the
Device class ensures that the device’s
confinement potential energy is updated accordingly).
# Load mesh, create device, load confinement potential
script_dir = pathlib.Path(__file__).parent.resolve()
path_mesh = str(script_dir/'meshes'/'gated_dot.msh2')
scaling = 1e-6
mesh = Mesh(scaling,path_mesh)
path_in = str(script_dir/'output'/'electric_potential_0.300V.hdf5')
d = Device(mesh)
d.set_temperature(0.1)
d.set_potential(io.load(path_in))
# Number of orbital states to consider in Schrodinger solver
num_states = 3
# Doping density
doping = 3e18*1e6 # In SI units
# Set up materials in heterostructure stack
d.new_region("cap", mt.GaAs)
d.new_region("barrier", mt.AlGaAs)
d.new_region("dopant_layer", mt.AlGaAs,ndoping=doping)
d.new_region("spacer_layer", mt.AlGaAs)
d.new_region("spacer_dot", mt.AlGaAs)
d.new_region("two_deg", mt.GaAs)
d.new_region("two_deg_dot", mt.GaAs)
d.new_region("substrate", mt.GaAs)
d.new_region("substrate_dot", mt.GaAs)
# Align the bands across hetero-interfaces using GaAs as the reference material
d.align_bands(mt.GaAs)
# Set up boundary conditions
# # Applied potentials
Vtop_1 = -0.5
Vtop_2 = 0.4
Vtop_3 = 0. # Apply 0V to top gate used for EDSR
Vbottom_1 = -0.5
Vbottom_2 = 0.4
Vbottom_3 = -0.5
# # Work function of the metallic gates at midgap of GaAs
barrier = 0.834*ct.e # n-type Schottky barrier height
# Set up regions and boundaries
d.new_schottky_bnd("top_gate_1", Vtop_1, barrier)
d.new_schottky_bnd("top_gate_2", Vtop_2, barrier)
d.new_schottky_bnd("top_gate_3", Vtop_3, barrier)
d.new_schottky_bnd("bottom_gate_1", Vbottom_1, barrier)
d.new_schottky_bnd("bottom_gate_2", Vbottom_2, barrier)
d.new_schottky_bnd("bottom_gate_3", Vbottom_3, barrier)
d.new_ohmic_bnd("ohmic_bnd")
# Create a submesh in which to solve Schrodinger's equation
submesh = SubMesh(mesh, ["spacer_dot","two_deg_dot","substrate_dot"])
subdevice = SubDevice(d,submesh)
6.3.3. Solving Schrödinger’s equation
Having defined the dot subdevice, the Schrödinger equation can be solved
in it, following the same method as in previous tutorials.
# Solve Schrodinger's equation for single electrons
params = SolverParams()
params.num_states = num_states
params.tol = 1e-9 # Decrease tolerance to stabilize Rabi frequency
s = Solver(subdevice,solver_params=params)
s.solve()
# # Plot eigenstates
# subdevice.print_energies()
# from device import analysis as an
# an.plot_slices(submesh, np.absolute(subdevice.eigenfunctions[:,0])**2,
# z=-42e-9)
# an.plot_slices(submesh, np.absolute(subdevice.eigenfunctions[:,1])**2,
# z=-42e-9)
# an.plot_slices(submesh, np.absolute(subdevice.eigenfunctions[:,2])**2,
# z=-42e-9)
The eigenenergies and eigenfunctions may be shown using the commented code.
6.4. The effect of a micromagnet on our quantum dot
6.4.1. Setting up the magnetic fields
The EDSR experiment we wish to simulate uses (1) a homogeneous magnetic field
aligned with the \(z\) axis, and (2) a weaker inhomogeneous
magnetic-field component aligned in the \(x\) direction, whose strength
varies linearly as a function of \(y\). This second magnetic-field
component may be produced by a micromagnet, and is typically called a
slanting field in the EDSR literature. Such a slanting field implements an
effective spin–orbit coupling for the electron confined within our quantum dot
(see Eq. (9.1) in the Quantum control theory
section). The homogeneous magnetic field is defined as homo_B
in the code below, while the slanting magnetic field is defined as
varying_B.
# Set up magnetic field
homo_B = np.array([0,0,0.6]) # Homogeneous B field.
# Position-dependent B - due to micromagnet
b = 0.3/1e-6 # 0.3 T/micron
y = submesh.glob_nodes[:,1]
varying_B = np.array([b*y, 0*y, 0*y]) # Varying B field depends linearly on the y coordinate.
Note that, while homo_B is defined as a 1D NumPy array, varying_B is
defined as a 2D NumPy array; one index runs over the indices of the global
nodes of the mesh while the others runs over \(x\), \(y\), and
\(z\).
6.4.2. The potential energy matrix
Having established the magnetic fields, we define a MicromagnetEDSR object
which models our device and the electron confined in the quantum dot subject to
a micromagnet-induced spin–orbit coupling.
# EDSR calculation
edsr = MicromagnetEDSR(d,subdevice)
The first argument of MicromagnetEDSR, d, represents the full
device, while the second argument, subdevice represents our quantum dot.
Modeling the dynamics of our qubit under EDSR will require solving the Poisson
equation on d and the Schrödinger equation only on subdevice.
We then use the solve method of the edsr object to compute the system
and drive Hamiltonian for our electron under EDSR (see the
Time-dependent Hamiltonian subsection of the Quantum control
theory section).
# Compute the matrix elements of the potential energy operator for
# different values of voltages appled to the 'top_gate_3'
print("Computing matrix elements of the potential")
gate_biases = np.linspace(-0.5, -0.4, num=6)
coupling_params = (homo_B, varying_B)
H0, UU = edsr.solve(gate_biases, 'top_gate_3', *coupling_params)
The edsr.solve method has three arguments. The first, gate_biases,
is an array containing the voltages we wish to apply to a gate. The second,
'top_gate_3', specifies the gate of d on which we apply the voltages
in gate_biases. Finally, *coupling_params is a tuple which specifies
our magnetic field, which will be used to characterize the engineered
spin–orbit coupling. The edsr.solve returns H0, a 2D array containing
the system Hamiltonian in its eigenbasis and shifted so that the ground state
has an energy of \(0\), and UU, a 3D array containing the drive
operator in the same basis as H0 for each element in gate_biases.
We can then plot the matrix elements of UU as a function of
gate_biases; we will do it in the specific case of UU[:, 0, 1]
(the first set of indices, :, runs over the elements of gate_biases).
Each matrix UU[i,:,:] corresponds to the confinement potential
perturbation \(\delta \hat V\) introduced in
Eq. (9.4) of the
Quantum control theory section, projected in the basis of the two
lowest-energy eigenstates of the Zeeman Hamiltonian, and for the
\(i\)-th applied gate bias.
# Check if the potential energy matrix is a linear function of the
# applied voltage
fig = plt.figure(figsize=(7,9))
ax1 = fig.add_subplot(2,1,1)
ax2 = fig.add_subplot(2,1,2)
ax1.plot(gate_biases, np.real(UU[:, 0, 1])/ct.e)
ax2.plot(gate_biases, np.imag(UU[:, 0, 1])/ct.e)
ax1.set_xlabel(r'$\varphi^\mathrm{bias}$', fontsize=16)
ax2.set_xlabel(r'$\varphi^\mathrm{bias}$', fontsize=16)
ax1.set_ylabel(r'$\mathrm{Real}[\delta V_{01}(\varphi^\mathrm{bias})]$ (eV)',
fontsize=16)
ax2.set_ylabel(r'$\mathrm{Imag}[\delta V_{01}(\varphi^\mathrm{bias})]$ (eV)',
fontsize=16)
plt.show(block=True)
This results in the figure below, which clearly shows that this matrix element
is a linear function of the applied voltage. It can similarly be verified that
all other elements of UU are linear functions of the applied voltage.
Fig. 6.4.9 The elements of the confinement potential perturbation matrix \(\delta \hat V\) are linear functions of the applied voltage
This linearity will greatly simplify our simulations of the dynamics of our spin qubit under a time-varying electric field: we can simply use a linear interpolation of the confinement potential perturbation matrix \(\delta \hat V\) to model the potential energy matrices for an arbitrary applied voltage.
In preparation of dynamics simulations, we define delta_V as the element
of UU corresponding to the largest voltage in gate_biases, and shift
its diagonal so as to ensure that delta_V[0,0] is 0. This arbitrary
potential energy shift does not influence the physical solution but will later
eliminate fast oscillating terms proportional to the identity matrix that
would otherwise slow down the solution of the time-dependent Schrödinger
equation.
# Use the fact that the perturbing potential is a linear function of the
# applied voltage.
delta_V = UU[-1]
# Reset the zero of energy (ground state has E!=0).
np.fill_diagonal(delta_V, delta_V.diagonal()-delta_V[0,0])
We can then print the system eigenenergies.
# System Hamiltonian in the eigenbasis is diagonal
print(f'evals (eV):{(np.diagonal(H0) - H0[0,0])/ct.e}')
This results in the following.
evals (eV):[0.00000000e+00 7.03947597e-05 5.78340423e-03 5.85380735e-03 1.15615160e-02 1.16319021e-02]
It can be noticed that the ground state and first-excited state are separated by merely \(70~\mu\textrm{eV}\) while the first-excited state and second-excited state are separated by a quantity two orders of magnitude greater. The spectrum is thus compatible with a well-behaved spin qubit. A good approximation to model the dynamics of our spin qubit is thus a projection onto these two states, ignoring higher-energy states.
6.5. EDSR—Projection on the 2-level subspace
To model the dynamics of such a two-level system, we use the Dynamics class.
# Simple EDSR calculation. The calculation is performed by projecting onto
# the relevant 2 dimensional subspace. Here we are studying transitions
# between the ground state and the first excited state. The optional
# parameter plot is set to true to visualize the results.
dyn = dynamics.Dynamics()
This class has a method, transition_2_levels, which projects the system
Hamiltonian and the drive Hamiltonian onto the appropriate two-level subspace
and compute its dynamics.
The transition_2_levels() method implements the theory described in the
subsection called A simple scenario: a driven two-level system under the rotating-wave approximation in the Quantum control theory section.
result, h0, u, omega0, omega_rabi = dyn.transition_2_levels(1, 0, H0,
delta_V, plot=True, npts=20000)
# Save the projections of H0 and delta_V on the 2-dimensional subspace
# for later use
np.save(str(script_dir/'output'/'h0.npy'),h0)
np.save(str(script_dir/'output'/'u.npy'),u)
The system and drive Hamiltonians are specified as H0 and delta_V,
respectively. The indices describing our two-level subspace are given as the
first two arguments of transition_2_levels: 1 is the state labeled as
\(\ket{\uparrow}\) and 0 is the state labeled as \(\ket{\downarrow}\),
the latter of which is taken as the initial state in the simulation. The
method transition_2_levels returns a number of variables: result is a
QuTiP object containing information about the dynamics, h0 is the system
Hamiltonian projected onto the two-level subspace, u is the drive
Hamiltonian projected onto the two-level subspace, omega0 is the resonance
frequency, and omega_rabi is the Rabi frequency of the transition between
the two states. Furthermore, since we have specified plot=True,
transition_2_levels produces the following figure.
Fig. 6.5.1 Rabi oscillations for our qubit driven by EDSR in the two-level subspace approximation
6.6. EDSR—Beyond the 2-level subspace approximation
We now use QuTiP to go beyond the 2-level subspace and rotating-wave approximations of the previous section. We start by defining the various variables which will be used in the QuTiP simulation.
First, QuTiP assumes units in which \(\hbar=1\). We therefore adjust H0
and delta_V accordingly, converting them to QuTiP Qobj objects.
############################################################################
# More accurate but longer calculation beyond the projection into the 2-level
# subspace.
H0 = qutip.Qobj(H0/ct.hbar) # system Hamiltonian in units of hbar
delta_V = qutip.Qobj(delta_V/ct.hbar) # potential energy matrix amplitude in units of hbar
To define our full Hamiltonian, H, as H0 plus a sinusoidal perturbation with
amplitude delta_V, we define the following.
H = [H0, [delta_V, 'cos(omega*t)']] # Full Hamiltonian - sinusoidal modulation of delta_V
Above, omega remains unspecified. We address this issue with the following.
# Drive the voltage at the resonance frequency of the transition we wish to
# observe
args = {'omega': omega0}
Next, we initialize the system in the ground state as follows.
# Initialize the system in its ground state
# The basis has 2*num_states; num_states orbital states and 2 spin states
psi0 = qutip.basis(2*num_states, 0) # initialize in the ground state
Having specified our Hamiltonian and initial state, we are ready to simulate
dynamics using the QuTiP mesolve function over a given list of times, which
we call t_list.
# Use the mesolve function of QuTiP to solve for the dynamics of the 6-level
# system.
T_Rabi = 2* np.pi / omega_rabi # Rabi cycle in the two-level limit
# array of times for which the solver should store the state vector
tlist = np.linspace(0, 2*T_Rabi, 100000)
result = qutip.mesolve(H, psi0, tlist, [], [], args=args, progress_bar=True)
It remains to save and plot the results, which we do as follows.
# Plot the results with matplotlib and save them
# using the EDSR module to get projectors onto the 6 basis states
Proj = dyn.projectors(2*num_states)
population_data = []
outfile = open(str(script_dir/'output'/'EDSR.pkl'), 'wb')
# Plot projections onto each state.
fig, axes = plt.subplots(2,num_states)
for i in range(2):
population_data_spin = []
for j in range(num_states):
expect = qutip.expect(Proj[num_states*i+j],result.states)
population_data_spin += [expect]
axes[i,j].plot(tlist, expect)
axes[i,j].set_xlabel(r'$t$', fontsize=20)
axes[i,j].set_ylabel(f'$P_{num_states*i+j}$', fontsize=20)
population_data += [population_data_spin]
pickle.dump((tlist,population_data), outfile)
fig.tight_layout()
plt.show()
This results in the following figure.
Fig. 6.6.1 Rabi oscillations for a spin qubit driven by EDSR beyond the two-level subspace approximation
We remark that while most of the population is exchanged between the two qubit states (\(0\) and \(1\)), some significant leakage into higher (non-computational) levels exists for the current drive amplitude and shape. This issue can typically be resolved using pulse shaping, for example so-called DRAG pulses. In addition, going beyond the rotating-wave oscillations, we see that the dynamics contains significant fast oscillatory components in addition to the ideal two-level Rabi oscillations considered above.
6.7. Full code
__copyright__ = "Copyright 2022, Nanoacademic Technologies Inc."
import numpy as np
import pathlib
from matplotlib import pyplot as plt
import qutip
import pickle
from device import constants as ct
from device.mesh3d import Mesh, SubMesh
from device import io
from device import Device, SubDevice
from device.solver_params import SolverParams
from device.schrodinger import Solver
from device import materials as mt
from qubit.edsr import MicromagnetEDSR
from qubit import dynamics
# Load mesh, create device, load confinement potential
script_dir = pathlib.Path(__file__).parent.resolve()
path_mesh = str(script_dir/'meshes'/'gated_dot.msh2')
scaling = 1e-6
mesh = Mesh(scaling,path_mesh)
path_in = str(script_dir/'output'/'electric_potential_0.300V.hdf5')
d = Device(mesh)
d.set_temperature(0.1)
d.set_potential(io.load(path_in))
# Number of orbital states to consider in Schrodinger solver
num_states = 3
# Doping density
doping = 3e18*1e6 # In SI units
# Set up materials in heterostructure stack
d.new_region("cap", mt.GaAs)
d.new_region("barrier", mt.AlGaAs)
d.new_region("dopant_layer", mt.AlGaAs,ndoping=doping)
d.new_region("spacer_layer", mt.AlGaAs)
d.new_region("spacer_dot", mt.AlGaAs)
d.new_region("two_deg", mt.GaAs)
d.new_region("two_deg_dot", mt.GaAs)
d.new_region("substrate", mt.GaAs)
d.new_region("substrate_dot", mt.GaAs)
# Align the bands across hetero-interfaces using GaAs as the reference material
d.align_bands(mt.GaAs)
# Set up boundary conditions
# # Applied potentials
Vtop_1 = -0.5
Vtop_2 = 0.4
Vtop_3 = 0. # Apply 0V to top gate used for EDSR
Vbottom_1 = -0.5
Vbottom_2 = 0.4
Vbottom_3 = -0.5
# # Work function of the metallic gates at midgap of GaAs
barrier = 0.834*ct.e # n-type Schottky barrier height
# Set up regions and boundaries
d.new_schottky_bnd("top_gate_1", Vtop_1, barrier)
d.new_schottky_bnd("top_gate_2", Vtop_2, barrier)
d.new_schottky_bnd("top_gate_3", Vtop_3, barrier)
d.new_schottky_bnd("bottom_gate_1", Vbottom_1, barrier)
d.new_schottky_bnd("bottom_gate_2", Vbottom_2, barrier)
d.new_schottky_bnd("bottom_gate_3", Vbottom_3, barrier)
d.new_ohmic_bnd("ohmic_bnd")
# Create a submesh in which to solve Schrodinger's equation
submesh = SubMesh(mesh, ["spacer_dot","two_deg_dot","substrate_dot"])
subdevice = SubDevice(d,submesh)
# Solve Schrodinger's equation for single electrons
params = SolverParams()
params.num_states = num_states
params.tol = 1e-9 # Decrease tolerance to stabilize Rabi frequency
s = Solver(subdevice,solver_params=params)
s.solve()
# # Plot eigenstates
# subdevice.print_energies()
# from device import analysis as an
# an.plot_slices(submesh, np.absolute(subdevice.eigenfunctions[:,0])**2,
# z=-42e-9)
# an.plot_slices(submesh, np.absolute(subdevice.eigenfunctions[:,1])**2,
# z=-42e-9)
# an.plot_slices(submesh, np.absolute(subdevice.eigenfunctions[:,2])**2,
# z=-42e-9)
# Set up magnetic field
homo_B = np.array([0,0,0.6]) # Homogeneous B field.
# Position-dependent B - due to micromagnet
b = 0.3/1e-6 # 0.3 T/micron
y = submesh.glob_nodes[:,1]
varying_B = np.array([b*y, 0*y, 0*y]) # Varying B field depends linearly on the y coordinate.
# EDSR calculation
edsr = MicromagnetEDSR(d,subdevice)
# Compute the matrix elements of the potential energy operator for
# different values of voltages appled to the 'top_gate_3'
print("Computing matrix elements of the potential")
gate_biases = np.linspace(-0.5, -0.4, num=6)
coupling_params = (homo_B, varying_B)
H0, UU = edsr.solve(gate_biases, 'top_gate_3', *coupling_params)
# Check if the potential energy matrix is a linear function of the
# applied voltage
fig = plt.figure(figsize=(7,9))
ax1 = fig.add_subplot(2,1,1)
ax2 = fig.add_subplot(2,1,2)
ax1.plot(gate_biases, np.real(UU[:, 0, 1])/ct.e)
ax2.plot(gate_biases, np.imag(UU[:, 0, 1])/ct.e)
ax1.set_xlabel(r'$\varphi^\mathrm{bias}$', fontsize=16)
ax2.set_xlabel(r'$\varphi^\mathrm{bias}$', fontsize=16)
ax1.set_ylabel(r'$\mathrm{Real}[\delta V_{01}(\varphi^\mathrm{bias})]$ (eV)',
fontsize=16)
ax2.set_ylabel(r'$\mathrm{Imag}[\delta V_{01}(\varphi^\mathrm{bias})]$ (eV)',
fontsize=16)
plt.show(block=True)
# Use the fact that the perturbing potential is a linear function of the
# applied voltage.
delta_V = UU[-1]
# Reset the zero of energy (ground state has E!=0).
np.fill_diagonal(delta_V, delta_V.diagonal()-delta_V[0,0])
# System Hamiltonian in the eigenbasis is diagonal
print(f'evals (eV):{(np.diagonal(H0) - H0[0,0])/ct.e}')
# Simple EDSR calculation. The calculation is performed by projecting onto
# the relevant 2 dimensional subspace. Here we are studying transitions
# between the ground state and the first excited state. The optional
# parameter plot is set to true to visualize the results.
dyn = dynamics.Dynamics()
result, h0, u, omega0, omega_rabi = dyn.transition_2_levels(1, 0, H0,
delta_V, plot=True, npts=20000)
# Save the projections of H0 and delta_V on the 2-dimensional subspace
# for later use
np.save(str(script_dir/'output'/'h0.npy'),h0)
np.save(str(script_dir/'output'/'u.npy'),u)
############################################################################
# More accurate but longer calculation beyond the projection into the 2-level
# subspace.
H0 = qutip.Qobj(H0/ct.hbar) # system Hamiltonian in units of hbar
delta_V = qutip.Qobj(delta_V/ct.hbar) # potential energy matrix amplitude in units of hbar
H = [H0, [delta_V, 'cos(omega*t)']] # Full Hamiltonian - sinusoidal modulation of delta_V
# Drive the voltage at the resonance frequency of the transition we wish to
# observe
args = {'omega': omega0}
# Initialize the system in its ground state
# The basis has 2*num_states; num_states orbital states and 2 spin states
psi0 = qutip.basis(2*num_states, 0) # initialize in the ground state
# Use the mesolve function of QuTiP to solve for the dynamics of the 6-level
# system.
T_Rabi = 2* np.pi / omega_rabi # Rabi cycle in the two-level limit
# array of times for which the solver should store the state vector
tlist = np.linspace(0, 2*T_Rabi, 100000)
result = qutip.mesolve(H, psi0, tlist, [], [], args=args, progress_bar=True)
# Plot the results with matplotlib and save them
# using the EDSR module to get projectors onto the 6 basis states
Proj = dyn.projectors(2*num_states)
population_data = []
outfile = open(str(script_dir/'output'/'EDSR.pkl'), 'wb')
# Plot projections onto each state.
fig, axes = plt.subplots(2,num_states)
for i in range(2):
population_data_spin = []
for j in range(num_states):
expect = qutip.expect(Proj[num_states*i+j],result.states)
population_data_spin += [expect]
axes[i,j].plot(tlist, expect)
axes[i,j].set_xlabel(r'$t$', fontsize=20)
axes[i,j].set_ylabel(f'$P_{num_states*i+j}$', fontsize=20)
population_data += [population_data_spin]
pickle.dump((tlist,population_data), outfile)
fig.tight_layout()
plt.show()