Simulating an open two-level quantum system¶
This tutorial provides a step-by-step guide to simulating the dynamics of a simple two-level quantum system subject to dissipation using the qruise-toolset. We will begin by simulating the system on its own, and then explore its behaviour when coupled to a leaky single-mode cavity.
1. Simulating a two-level system¶
The two-level system (TLS) we are considering is described by the following Hamiltonian:
$$ H = 2 \pi \beta \sigma_x , \tag{1} $$
where $\beta$ is the tunnelling rate, describing the coupling strength that drives transitions between the TLS's ground and excited states. $\sigma_x$ is the Pauli-X operator, which we'll import from the QuTip library.
While closed quantum systems are governed by the Schrödinger equation, open systems require the Lindblad master equation, which account for both unitary dynamics and dissipative processes. This is given by
$$ \frac{\partial \rho(t)}{\partial t} = -i \left[ H(t), \rho(t) \right] + \sum_i \gamma_i \left( L_i \rho(t) L_i^\dagger - \frac{1}{2} \left\{ L_i^\dagger L_i, \rho(t) \right\} \right) , \tag{2} $$
where $\rho(t)$ is the density matrix of the system at time $t$, $\gamma_i$ are decay rates (associated with decoherence or dissipation), and $L_i$ are the corresponding collapse operators. In our system, we'll use a single collapse operator to describe dissipation, with $\gamma_i = \gamma = 0.05 \text{ s}^{-1}$ and $L_i = L =\sigma_x$.
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)
We can start by defining the time parameters and decay rates for the simulation.
import jax.numpy as jnp
t0 = 0.0 # start time of simulation
tfinal = 10.0 # end time of simulation
N = 200 # number of time points
ts = jnp.linspace(t0, tfinal, N) # time array
beta = 0.1 # tunnelling rate
gamma = 0.05 # dissipation rate
We then create our Hamiltonian and collapse operator (including the dissipation rate).
import qruise.toolset as qr
import qutip as qt
# define Hamiltonian and collapse operator
H = qr.Hamiltonian(2 * jnp.pi * beta * qt.sigmax())
c_op = jnp.sqrt(gamma) * qt.sigmax()
To solve the Lindblad 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 qubit state (
pho0, here the ground state) - the time range of the simulation (
t0totfinal).
We can then use Qruise's ODESolver with the collapse operator (c_op) to solve the Lindblad master equation and calculate the wavefunction of the system at each timestamp.
from diffrax import Tsit5
# define initial qubit state (in matrix form)
rho0 = jnp.asarray(qt.basis(2, 0).full())
rho0 = rho0 @ rho0.T # density matrix
# define Problem
prob = qr.Problem(H, rho0, {}, (t0, tfinal))
# use ODEsolver to solve Lindblad master equation
solver = qr.ODESolver(Tsit5(), dt0=None, saveat={"ts": ts}, atol=1e-8, rtol=1e-6)
solver.set_system(prob.lindblad([c_op]))
rhos = solver.evolve(*prob.problem())
We can now calculate and plot the expectation values of $\sigma_z$ and $\sigma_y$ ($\langle \sigma_z \rangle$ and $\langle \sigma_y \rangle$, respectively).
# 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()
Nice! You can see that $\langle \sigma_z \rangle$ and $\langle \sigma_y \rangle$ exhibit periodic oscillations driven by tunnelling (governed by $\beta$) between the $|0\rangle$ and $|1\rangle$ states, with their amplitudes decaying over time due to dissipation (governed by $\gamma$).
2. Simulating a two-level system coupled to a cavity¶
Now let's take this one step further by considering our TLS as a two-level atom that we'll couple to a leaky single-mode cavity via a dipole-type interaction. This supports a coherent exchange of energy between the TLS and the cavity.
The Hamiltonian for this system is:
$$ H = \underbrace{2 \pi \sigma^+ \sigma^- \vphantom{\big|}}_{\text{TLS energy}} + \underbrace{2 \pi a^\dagger a \vphantom{\big|}}_{\text{cavity energy}} + \underbrace{2 \pi g (a \sigma^++ a^\dagger \sigma^-) \vphantom{\big|}}_{\text{TLS-cavity interaction}} , \tag{3} $$
where
- $\sigma^+$ and $\sigma^-$ are the raising and lowering operators for the TLS
- $a^\dagger$ and $a$ are the creation and annihilation operators for the single-mode cavity
- $g$ is the TLS-cavity coupling strength.
The Lindblad master equation (Equation 2) stays the same, except we use the new Hamiltonian and a new collapse operator. Here, we have $L_i=L=a$ (instead of $L=\sigma_x$), and $\gamma_i = \kappa$, which represents the leakage rate from the cavity.
We'll use the same time parameters as in the previous simulation, but we need to define the system parameters $g$ and $\kappa$.
g = 0.25 # TLS-cavity coupling strength
kappa = 0.1 # cavity leakage rate
Then we can define our operators and Hamiltonian.
# define TLS and cavity operators
atomic = qt.tensor(qt.destroy(2), qt.qeye(10)) # TLS lowering operator
cavity = qt.tensor(qt.qeye(2), qt.destroy(10)) # cavity annihilation operator
# define Hamiltonian and collapse operator
H = qr.Hamiltonian(
2 * jnp.pi * atomic.dag() * atomic
+ 2 * jnp.pi * cavity.dag() * cavity
+ 2 * jnp.pi * g * (cavity * atomic.dag() + cavity.dag() * atomic)
)
c_op = jnp.sqrt(kappa) * cavity
Note: The Hilbert dimension of the TLS is 2 ($|0\rangle$ and $|1\rangle$) and we'll set it to 10 for the cavity (up to 10 photons).
Now we just need to define our initial system state (rho0), which we'll take as the ground state ($|0\rangle$) for the TLS and as a 5-photon Fock state ($|5\rangle$) for the cavity. We can then define the Problem and solve the Lindblad master equation using Qruise's ODESolver as we did in the previous section.
# define inital system state (in matrix form)
rho0 = jnp.asarray(qt.tensor(qt.fock(2, 0), qt.fock(10, 5)).full())
rho0 = rho0 @ rho0.T
# define Problem
prob = qr.Problem(H, rho0, {}, (t0, tfinal))
# use ODEsolver to solve Lindblad master equation
solver = qr.ODESolver(Tsit5(), dt0=None, saveat={"ts": ts}, atol=1e-8, rtol=1e-6)
solver.set_system(prob.lindblad([c_op]))
rhos = solver.evolve(*prob.problem())
We can now calculate and plot the expectation values of $\sigma^+ \sigma^-$ and $a^\dagger a$ ($\langle \sigma^+ \sigma^- \rangle$ and $\langle a^\dagger a \rangle$, respectively).
# expectation values of sigma^+ sigma^- and a^\dagger a
atomic_dag_atomic_exp = qr.op_expectation(atomic.dag() * atomic, rhos)
cavity_dag_cavity_exp = qr.op_expectation(cavity.dag() * cavity, rhos)
# plotting the expectation values
canvas.init_canvas(
x_axis_label="Time [a.u.]", y_axis_label="Operator expectation value"
)
canvas.plot(ts, cavity_dag_cavity_exp, labels=["cavity photon number"])
canvas.plot(ts, atomic_dag_atomic_exp, labels=["atomic excitation probability"])
canvas.show_canvas()
Great! You can see we still see periodic oscillation, driven by the coupling between the TLS and the cavity ($g$), with decaying amplitudes due to dissipation ($\gamma$).
Congratulations! You now know how to simulate an open two-level quantum system.