Simulating a closed two-level quantum system¶
This tutorial provides a step-by-step guide to simulating the dynamics of a simple two-level quantum system using the qruise-toolset.
We'll begin by defining the system Hamiltonian, then demonstrate how to simulate pure states with the Schrödinger equation and mixed states using the von Neumann equation.
1. Defining the Hamiltonian¶
The two-level system (TLS) we are considering is described by the following Hamiltonian:
$$ H(t) = \frac{\omega}{2}\sigma_z + c(t)\sigma_x, \tag{1} $$
where $\omega$ is the TLS resonance frequency, $\sigma_z$ and $\sigma_x$ are Pauli matrices, and $c(t)$ is the drive function. For this we'll use a Gaussian pulse given by:
$$ c(t; \{a, \sigma, \mu \}) = \frac{a}{\sigma\sqrt{2\pi}} \exp\Big[-\frac{(t - \mu)^2}{2\sigma^2}\Big], \tag{2} $$
which is characterised by three parameters:
$a$, a scalar that adjusts the amplitude of the pulse
$\sigma$, the pulse variance
$\mu$, the total pulse duration.
We can start by defining the stationary part of our Hamiltonian, $H_0$, i.e., the first term. For this, we can use sigmaz from the QuTiP library.
Tip: Make sure that your environment is correctly set to handle float64 precision by setting JAX_ENABLE_X64=True or add the following codeblock to your scripts' preamble:
import jax
jax.config.update("jax_enable_x64", True)
import qruise.toolset as qr
import qutip as qt
omega = 2.0 # TLS frequency
# stationary Hamiltonian, H0
H = qr.Hamiltonian(omega / 2 * qt.sigmaz())
Before we add the time-dependent term in our Hamiltonian, we need to define our Gaussian drive pulse.
import jax.numpy as jnp
def drive(t, params):
"""Gaussian pulse"""
a = params["a"] # pulse amplitude
sigma = params["sigma"] # variance (controls pulse width)
factor = a / jnp.sqrt(2.0 * jnp.pi) / sigma
return factor * jnp.exp(-((t - tfinal / 2) ** 2) / (2.0 * sigma**2))
# Note: mu is set to tfinal/2
Note: The signature of the pulse (envelope), in our case the drive function, must strictly satisfy the (t, params) pattern, where t is float and params is Dict (dictionary). It's important to ensure that your functions respect this constraint.
We now need to define the time and pulse parameters for the simulation.
t0 = 0.0 # start time of simulation
tfinal = 1.0 # end time of simulation
N = 200 # number of time points
ts = jnp.linspace(t0, tfinal, N) # time array
# Parameters of the Gaussian pulse
params = {"a": 2.14, "sigma": 0.8 * tfinal}
It’s a good idea to check your drive behaves as desired before proceeding, so let’s quickly plot it.
canvas = qr.PlotUtil(
x_axis_label="t [a.u.]", y_axis_label="Amplitude [a.u.]", notebook=True
)
canvas.plot(ts, drive(ts, params), labels=["Pulse"])
canvas.show_canvas()
Great, our drive looks as expected!
Now, we simply add the drive term to our Hamiltonian.
H.add_term(qt.sigmax(), drive)
2. Simulating pure states using the Schrödinger equation¶
Pure states in a two-level quantum system can be simulated by solving the Schrödinger equation, which describes the unitary evolution of the system's wavefunction over time:
$$ i \frac{\partial}{\partial t} \psi(t) = H(t) \psi(t). \tag{3} $$
To solve the Schrödinger equation, we need to specify the equations and parameters that govern the system. In the qruise-toolset, we combine these using a Problem, which is instantiated with:
- the Hamiltonian (
H) - the initial system state (
y0, here the ground state) - the pulse parameters (
params) - the time range of the simulation (
t0totfinal).
Note: Generally speaking, a Problem is a set of differential or algebraic equations that describe the dynamics of a system. There are well-known equations that describe the dynamics of quantum systems, e.g. the Schrödinger equation for pure states, the von Neumann master equation for mixed states, or the Lindblad equation for open quantum systems with dissipation. The cornerstone of all these equations is the Hamiltonian.
# define initial system state
y0 = jnp.array([1.0, 0.0], jnp.complex128)
# define Problem
prob = qr.Problem(H, y0, params, (t0, tfinal))
We can use Qruise's ODESolver to solve the Schrödinger equation (which is essentially an ordinary differential equation (ODE)), and calculate the wavefunction of the system at each timestamp. For this, we employ the adaptive Runge-Kutta 4(5) method, implemented in JAX through the diffrax library's Tsit5 solver.
from diffrax import Tsit5
# use ODEsolver to solve Schrödinger equation
solver = qr.ODESolver(
Tsit5(), # Runge-Kutta method
dt0=None, # differential of time, set to None as we solve adaptively
saveat={"ts": ts}, # timestamps at which we want to save the solution
atol=1e-8, # Absolute tolerance for the adaptive solver
rtol=1e-6, # Relative tolerance for the adaptive solver
)
solver.set_system(prob.schroedinger()) # specify we want to solve Schrödinger equation
res = solver.evolve(*prob.problem())
We can view the results of our simulation by calculating and plotting the expectation value of the state populations against time.
# calculate the expectation value of each state
pop = qr.get_population(res)
# plot populations against time
canvas.init_canvas(x_axis_label="t [a.u.]", y_axis_label="Population")
canvas.plot(ts, pop, labels=["ground state", "1st excited state"])
canvas.show_canvas()
Nice! You can see we get a (partial) population exchange between the ground and excited states of our TLS.
3. Simulating mixed states using the von Neumann equation¶
Mixed states in a two-level quantum system cannot be simulated by solving the Schrödinger equation. For this we need the von Neumann equation:
$$ \frac{\partial \rho(t)}{\partial t} = -i \left[ H(t), \rho(t) \right] , \tag{4} $$
which describes the time evolution of the system's density matrix, $\rho(t)$.
We start by defining the initial system state, rho0, this time as a density matrix. We then need to redefine our Problem using rho0 instead of y0. (We can leave all other variables as they are.)
# define initial system state (ground state)
rho0 = jnp.asarray(qt.basis(2, 0).full())
rho0 = rho0 @ rho0.T # density matrix
# define Problem
prob = qr.Problem(H, rho0, params, (t0, tfinal))
We can use the solver we defined earlier, which employs the Qruise ODESolver. Then all we need to do is tell the solver to use the von Neumann equation instead of the Schrödinger equation.
solver.set_system(prob.von_neumann()) # specify we want to solve von Neumann equation
rhos = solver.evolve(*prob.problem())
Plotting the state populations will only give us information about the diagonal elements of the density matrix. This means we lose information on the coherence of the system, which is contained in the off-diagonal elements. For mixed states, it's therefore better to plot the expectation values of the Pauli matrices, $\langle \sigma_z \rangle$ and $\langle \sigma_y \rangle$.
# calculate expectation values of sigmaz and sigma y
sigmaz_exp = qr.op_expectation(qt.sigmaz(), rhos)
sigmay_exp = qr.op_expectation(qt.sigmay(), rhos)
# plot expectation values against time
canvas = qr.PlotUtil(
x_axis_label="Time [a.u.]", y_axis_label="Operator expectation value", notebook=True
)
canvas.plot(ts, sigmaz_exp, labels=["<sigma_z>"])
canvas.plot(ts, sigmay_exp, labels=["<sigma_y>"])
canvas.show_canvas()
Great! Now you know how to simulate a closed two-level system using both the Schrödinger equation and the von Neumann equation.