Intorduction to CUDA-Q Dynamics (Jaynes-Cummings Model)

Why dynamics simulations vs. circuit simulations?

CUDA-Q Dynamics is a powerful library enabling GPU-accelerated simulations of the time evolution of open and closed quantum systems. In contrast to quantum circuit simulation, which simulates the algorithmic output of the application of discrete logical quantum gates (eg X, Y, Z, CNOT, H, etc), dynamics simulations represent how quantum systems evolve in time, revealing important insights about the non-idealities of a quantum system, including how they interact with the environment and other quantum systems.

To make a classical analogy, the logic of a classical computer can be modeled using binary logic gates applied to transistor states, analogous to simulating the application of quantum logic gates to qubits. Dynamics simulations are an abstraction layer below this. Classically, to design better and higher-performing transistors, it’s useful to fully model the device physics of the transistor itself, including fluctuations in voltage, capacitance, and current. Analogously, quantum dyanmics simulations can fully capture the physics of components beyond just the qubits, including drive pulses, resonators, couplers, and noise.

Figure 1. Classical computing can be a helpful analogy for understanding the difference between algorithm and dynamics simulations. Image credits: Fredrik Brange (top right) and qutip-qip (bottom right)

Functionality

CUDA-Q dynamics provides an easy API, evolve, to solve the Schrodinger and Lindblad Master equation, as well as number of commonly used operators to construct Hamiltonians. The basic order of operations is similar to other quantum dynamics frameworks. 1. Define: - The dimensions of the system - The Hamiltonian - The initial state of the system - Any dissipation terms - The timesteps desired 2. Evolve!

Conveniently, this is the same programming model and API used to submit CUDA-Q jobs to analog quantum computers like QuEra and Pasqal. For simulation workloads, we also provide several commonly used numerical integrators.

In addition to this higher level functionality, we also provide a lower level library called cuDensityMat which provide the building blocks for accelerating custom solvers and other dynamics frameworks.

Performance

CUDA-Q Dynamics shines for large system sizes, as is the case for most GPU workloads. We will see maximal speed-ups over CPU implementations with Hilbert spaces larger than \(O(1000)\) levels. For systems of this size and larger, it is common to see speed-ups of >1000x. A benchmark modeling 100 timesteps of a system consisting of a 32-level transmon, 256-level resonator, and 4-level Purcell filter yielded a 1,140x speed-up in execution from single CPU (dual-socket Intel Xeon 8480CL) to single GPU (H100), or 1.5 days to 2 minutes.

CUDA-Q Dynamics also provides multi-GPU and multi-node capabilities to both increase the performance and scale of simulations. Additional GPUs can be pooled for additional memory to simulate larger systems (weak-scaling). Alternatively, with a fixed system size, additional GPUs can be pooled for increase performance (strong-scaling).

Full details of the benchmarks can be found in the blog here.

Section 1 - Simulating the Jaynes-Cummings Hamiltonian

The Jaynes-Cummings (JC) Hamiltonian is a important theoretical model describing a two-level system interacting with a resonator. It describes how an atom or qubit fundamentally interacts with light. It is a powerful building block for modeling and designing many superconducting qubit systems and can be extended to describe many of the interactions in a superconducting qubit system. In practice, it describes the physics by which qubits can be controlled by external control pulses and the readout process used to determine their state.

Figure 2. A two-level quantum system such as a qubit or an atom interacting with a light field confined by a resonator. The JC model describes this interaction. Image credit: Wikipedia.

The Hamiltonian can be simplified into the following three principal terms, by taking the rotating wave approximation, and dropping constant terms.

$:nbsphinx-math:hat{H}{JC} = :nbsphinx-math:`hat{H}`{cavity} +:nbsphinx-math:hat{H}{qubit} +:nbsphinx-math:`hat{H}`{interaction} $

with

\(\hat{H}_{cavity} = \hbar \omega_c \hat a^\dagger_c \hat a_c\)

\(\hat{H}_{qubit} = \hbar \omega_q \hat \sigma_+ \hat \sigma_-\)

\(\hat{H}_{interaction} = \hbar g_c (\hat \sigma_+ \hat a_c + \hat \sigma_- \hat a^\dagger_c)\)

In this section, we will simulate JC Hamiltonian. You will learn to: - Define the dimensions of the system - Construct the Hamiltonian with CUDA-Q operators - Define the initial state of the system - Execute the simulation with cudaq.evolve()

# Import necessary packages
import cudaq
from cudaq import boson, Schedule, ScipyZvodeIntegrator, spin, RungeKuttaIntegrator
import numpy as np
import cupy as cp
import os
import matplotlib.pyplot as plt

We start by setting the backend target to dynamics. Then, we define the degrees of freedom, specifying the number of levels associated with each one. We can then construct the Hamiltonian by applying the correct operators to the relevant degrees of freedom.

# Set the target to our dynamics simulator
cudaq.set_target("dynamics")

# System dimensions: atom (2-level system) and cavity (2-level system)
dimensions = {0: 2, 1: 2}

# Alias for commonly used operators
# Cavity operators
a = boson.annihilate(1)
a_dag = boson.create(1)

# Atom operators
sm = boson.annihilate(0)
sm_dag = boson.create(0)

# Defining the Hamiltonian for the system: self-energy terms and cavity-atom interaction term.
H_cav = 2 * np.pi * boson.number(1)
H_qubit = 2 * np.pi * boson.number(0)
H_int = 2 * np.pi * 0.25 * (sm * a_dag + sm_dag * a)
hamiltonian = H_cav + H_qubit + H_int

Initialize the states of the system in cupy arrays and define the time chedule

# Initial state of the system
# Qubit in ground state
qubit_state = cp.array([[1.0, 0.0], [0.0, 0.0]], dtype=cp.complex128)

# Cavity in a state which has 1 photon initially
cavity_state = cp.zeros((2, 2), dtype=cp.complex128)
cavity_state[1][1] = 1.0

# Construct the Density Matrix
rho0 = cudaq.State.from_data(cp.kron(qubit_state, cavity_state))

# Define our timesteps
steps = np.linspace(0, 10, 201)
schedule = Schedule(steps, ["time"])

Evolve! Various observables can be collected at each time step. For this, we will track the population of both the qubit and the cavity in the absence of noise.

# Evolve the system
evolution_result = cudaq.evolve(
    hamiltonian,
    dimensions,
    schedule,
    rho0,
    observables=[boson.number(1), boson.number(0)],
    collapse_operators=[],
    store_intermediate_results=cudaq.IntermediateResultSave.EXPECTATION_VALUE,
    integrator=ScipyZvodeIntegrator())

# Collect results
get_result = lambda idx, res: [
    exp_vals[idx].expectation() for exp_vals in res.expectation_values()
]

ideal_results = [
    get_result(0, evolution_result),
    get_result(1, evolution_result)
]

# Plot the results
plt.plot(steps, ideal_results[0])
plt.plot(steps, ideal_results[1])
plt.ylabel("Expectation value")
plt.xlabel("Time")
plt.legend(("Cavity Photon Number", "Qubit Excitation Probability"))
plt.title("No decay")

Text(0.5, 1.0, 'No decay')

We see the qubit and cavity coherently exchanging population with a Rabi period of 2. Intuitively, a photon is being absorbed by the qubit and then re-emitted to the cavity.

Recall that we set the coupling rate in the interation Hamiltonian to \(g=\pi/2\). From performing some analysis on the JC Hamiltonian, one can derive that this oscillation rate is \(\Omega_n = 2 g \sqrt{n +1}\), where the effective subspace is spanned by \(\{\ket{e,n},\ket{g,n+1}\}\). Is this consistent with what we observe?

Exercise 1 - Simulating a many-photon Jaynes-Cummings Hamiltonian

Let us now adjust the code to simulate a larger system. Adjust the code below to now simulate 20 photons in the cavity.

# System dimensions: atom (2-level system) and cavity (20 photon system)
dimensions = {0: 2, 1: 21}

# Initial state of the system
# Qubit in ground state
qubit_state = cp.array([[1.0, 0.0], [0.0, 0.0]], dtype=cp.complex128)


# Cavity in a state which has 1 photons initially
cavity_state = cp.zeros((21, 21), dtype=cp.complex128)
cavity_state[20][20] = 1.0

# Construct the Density Matrix
rho0 = cudaq.State.from_data(cp.kron(qubit_state, cavity_state))

# Evolve the system
evolution_result = cudaq.evolve(
    hamiltonian,
    dimensions,
    schedule,
    rho0,
    observables=[boson.number(1), boson.number(0)],
    collapse_operators=[],
    store_intermediate_results=cudaq.IntermediateResultSave.EXPECTATION_VALUE,
    integrator=ScipyZvodeIntegrator())

# Collect results
get_result = lambda idx, res: [
    exp_vals[idx].expectation() for exp_vals in res.expectation_values()
]

ideal_results = [
    get_result(0, evolution_result),
    get_result(1, evolution_result)
]

# Plot the results
plt.plot(steps, ideal_results[0])
plt.plot(steps, ideal_results[1])
plt.ylabel("Expectation value")
plt.xlabel("Time")
plt.legend(("Cavity Photon Number", "Qubit Excitation Probability"))
plt.title("20 Photons")

Text(0.5, 1.0, '20 Photons')

What do you observe about the new Rabi frequency? Is this consistent with what you expect?

The Rabi frequency is proportional to \(2 g \sqrt{n+1}\), so \(\sqrt{20}=4.5\) times faster than the 1 photon case

Section 2 - Simulating open quantum systems with the collapse_operators

In the previous section we simulated an idealized quantum system, isolated from the environment, with no noise. This is equivalent to solving the Schrodinger equation. Let us now add the effects of the environment, by adding noise to the system via the collapse_operators. We will now solve the Lindbladian:

\(\dot \rho = -\frac{i}{\hbar} [H,\rho]+ \sum_i \gamma_i (L_i \rho L^\dagger_i-\frac{1}{2}\{L^\dagger_i L_i,\rho\})\)

where we can assign the dissipation terms \(L_i\) and rates \(\gamma_i\).

In this section we will learn to use the collapse_operators to model noise in the system.

We can add a list of dissipation terms to the collapse_operators argument in cudaq.evolve. Let’s add a photon loss mechanism with rate 0.1 with the a annihilation operator and observe what happens.

gamma=np.sqrt(0.1)

# Then, evolve the system with a collapse operator modeling cavity decay (leaking photons)
evolution_result_decay = cudaq.evolve(
    hamiltonian,
    dimensions,
    schedule,
    rho0,
    observables=[boson.number(1), boson.number(0)],
    collapse_operators=[gamma * a],
    store_intermediate_results=cudaq.IntermediateResultSave.EXPECTATION_VALUE,
    integrator=ScipyZvodeIntegrator())

get_result = lambda idx, res: [
    exp_vals[idx].expectation() for exp_vals in res.expectation_values()
]

decay_results = [
    get_result(0, evolution_result_decay),
    get_result(1, evolution_result_decay)
]

plt.plot(steps, decay_results[0])
plt.plot(steps, decay_results[1])
plt.ylabel("Expectation value")
plt.xlabel("Time")
plt.legend(("Cavity Photon Number", "Qubit Excitation Probability"))
plt.title("With decay")

Text(0.5, 1.0, 'With decay')

Exercise 2 - Adding additional jump operators \(L_i\)

2a. Now instead of considering a photon loss mechanism, consider a photon excitation mechanism with rate 0.1. Initialize the cavity with one photon.

# System dimensions: atom (2-level system) and cavity (20 photon system)
dimensions = {0: 2, 1: 21}

# Initial state of the system
# Qubit in ground state
qubit_state = cp.array([[1.0, 0.0], [0.0, 0.0]], dtype=cp.complex128)

# Cavity in a state which has 1 photon initially
cavity_state = cp.zeros((21, 21), dtype=cp.complex128)
cavity_state[1][1] = 1.0

# Construct the Density Matrix
rho0 = cudaq.State.from_data(cp.kron(qubit_state, cavity_state))

gamma=np.sqrt(0.1)
collapse_operators=[gamma * a_dag]

# Then, evolve the system with a collapse operator modeling cavity decay (leaking photons)
evolution_result_excitation = cudaq.evolve(
    hamiltonian,
    dimensions,
    schedule,
    rho0,
    observables=[boson.number(1), boson.number(0)],
    collapse_operators=collapse_operators,
    store_intermediate_results=cudaq.IntermediateResultSave.EXPECTATION_VALUE,
    integrator=ScipyZvodeIntegrator())

get_result = lambda idx, res: [
    exp_vals[idx].expectation() for exp_vals in res.expectation_values()
]

decay_results = [
    get_result(0, evolution_result_excitation),
    get_result(1, evolution_result_excitation)
]

plt.plot(steps, decay_results[0])
plt.plot(steps, decay_results[1])
plt.ylabel("Expectation value")
plt.xlabel("Time")
plt.legend(("Cavity Photon Number", "Qubit Excitation Probability"))
plt.title("With excitation")

Text(0.5, 1.0, 'With excitation')

2b. Now add a dephasing term to qubit with rate dephasing=0.4.

gamma=np.sqrt(0.1)
dephasing=np.sqrt(0.4)
collapse_operators=[gamma * a_dag, dephasing * spin.z(0)]

# Then, evolve the system with a collapse operator modeling cavity decay (leaking photons)
evolution_result_excitation = cudaq.evolve(
    hamiltonian,
    dimensions,
    schedule,
    rho0,
    observables=[boson.number(1), boson.number(0)],
    collapse_operators=collapse_operators,
    store_intermediate_results=cudaq.IntermediateResultSave.EXPECTATION_VALUE,
    integrator=ScipyZvodeIntegrator())

get_result = lambda idx, res: [
    exp_vals[idx].expectation() for exp_vals in res.expectation_values()
]

decay_results = [
    get_result(0, evolution_result_excitation),
    get_result(1, evolution_result_excitation)
]

plt.plot(steps, decay_results[0])
plt.plot(steps, decay_results[1])
plt.ylabel("Expectation value")
plt.xlabel("Time")
plt.legend(("Cavity Photon Number", "Qubit Excitation Probability"))
plt.title("With excitation and qubit dephasing")

Text(0.5, 1.0, 'With excitation and qubit dephasing')

Section 3 - Many qubits coupled to the resonator

Let us now construct a model of several qubits coupled to a resonator. This is known as the Tavis-Cummings model. We define the dimensions of the system, construct the Hamiltonian assuming a distribution of qubit-resonator couplings near \(\pi/2\), qubit energies near \(2\pi\), and photon energies of \(2\pi\). Initialize all qubits in the ground state and with the resonator with 10 photons.

import random

# System dimensions: atom (2-level system) and cavity (2-level system)
no_qubits=6
dimensions = {k: 11 if k == 0 else 2 for k in range(no_qubits+1)}

# Define our timesteps
steps = np.linspace(0, 10, 101)
schedule = Schedule(steps, ["time"])

# Cavity operators
a = boson.annihilate(0)
a_dag = boson.create(0)

# Defining the Hamiltonian for the system: self-energy terms and cavity-atom interaction term.
H_cav = 2 * np.pi * boson.number(0)

# Qubit and interaction Hamiltonians
def H_qubit_i(qubit, energy):
    return energy * boson.number(qubit)
def H_int_i(qubit, coupling):
    return coupling * (boson.annihilate(qubit) * a_dag + boson.create(qubit) * a)

# Qubit and coupling energies
coupling = [np.pi/4 * random.gauss(1,0.1) for i in range(no_qubits)]
qubit_energy = [2 * np.pi * random.gauss(1,0.1) for i in range(no_qubits)]

hamiltonian = H_cav
for qubit in range(1,no_qubits+1):
    hamiltonian+=H_qubit_i(qubit, qubit_energy[qubit-1])+H_int_i(qubit, coupling[qubit-1])

# Initial state of the system

# Cavity in a state which has 10 photons initially
cavity_state = cp.zeros((11, 11), dtype=cp.complex128)
cavity_state[10][10] = 1.0

# Qubits in ground state
qubit_g = cp.array([[1.0, 0.0], [0.0, 0.0]], dtype=cp.complex128)

# Qubit in excited state
qubit_e = cp.array([[0.0, 0.0], [0.0, 1.0]], dtype=cp.complex128)

composite_state = cavity_state
for i in range(1,no_qubits+1):
    composite_state = cp.kron(composite_state,qubit_g)

# Construct the Density Matrix
rho0 = cudaq.State.from_data(composite_state)

# Optional: Add dissipation terms for each term
T1=[abs(random.gauss(100,100)) for i in range(no_qubits+1)]
T1[0]=50

gamma=[1/np.sqrt(T1[i]) for i in range(no_qubits+1)]
collapse_operators=[gamma[i]*boson.annihilate(i) for i in range(no_qubits+1)]

# Define the collapse operators
collapse_operators=[]

# Evolve the system
evolution_result = cudaq.evolve(
    hamiltonian,
    dimensions,
    schedule,
    rho0,
    observables=[boson.number(i) for i in range(no_qubits+1)],
    collapse_operators=collapse_operators,
    store_intermediate_results=cudaq.IntermediateResultSave.EXPECTATION_VALUE,
    integrator=ScipyZvodeIntegrator())

get_result = lambda idx, res: [
    exp_vals[idx].expectation() for exp_vals in res.expectation_values()
]

tc_results = [
    get_result(i, evolution_result) for i in range(no_qubits+1)
]
plt.plot()

for i in range(no_qubits):
    plt.plot(steps, tc_results[i])
plt.ylabel("Expectation value")
plt.xlabel("Time")
plt.legend(("Cavity Photon Number", "Qubit Excitation Probability"))
plt.title("Many qubits coupled to a resonator")

Text(0.5, 1.0, 'Many qubits coupled to a resonator')

We now see the effect of multiple qubits exchanging excitations at different rates with the photons in the resonator