[1]:

import quanguru as qg
import numpy as np
import matplotlib.pyplot as plt
import platform

18 - Jaynes-Cummings Open System Dynamics#

The Jaynes-Cummings Hamiltonian is written as

\(H_{JC} = \hbar\omega_{c} a^{\dagger}a + \frac{1}{2}\hbar\omega_{q}\sigma_{z} + \hbar g(a^{\dagger}\sigma_{-} + a\sigma_{+})\)

where \(\sigma_{\pm} = (\sigma_{x} \pm i\sigma_{y})/2\) are raising/lowering operators for a two-level system, \(\sigma_{\mu}\) are the Pauli spin operators with \(\mu\in\{x,y,z\}\), \(a^{\dagger}\) and \(a\) are the creation and annihilation operators for the field mode, and \(\omega_{c}\), \(\omega_{q}\), and \(g\) are the cavity-field, qubit, and coupling (angular-) frequencies, respectively.

We will use the below parameters in our simulation

[2]:

# parameters for the Hamiltonian
qubitFreq = 1
cavityFreq = 1
couplingFreq = 0.25
cavityDim = 5

# parameters for the evolution
totalTime = 3*(np.pi/couplingFreq)
timeStep = 0.1

and we can describe the JC-Hamiltonian/system in QuanGuru as

[3]:

# Qubit for the JC model
QubitJC = qg.Qubit(frequency=qubitFreq)

# Cavity for the JC model
CavityJC = qg.Cavity(dimension=cavityDim, frequency=cavityFreq)

# JC model consists of a qubit and cavity
# and this is the 'free evolution' part of the JC-Hamiltonian
JCSystem = QubitJC + CavityJC

The interaction Hamiltonian consists of two terms:

  • \(a^{\dagger}\sigma_{-}\): cavity creation and qubit lowering

  • \(a\sigma_{+}\): cavity annihilation and qubit raising

When creating multi-system terms ensure the order of the operators is the same as the QuantumSystem objects.

e.g.

couplingTerm = totalSystem.createTerm(
    qSystem=[quantumSystem1, quantumSystem2],
    operator=[operator1, operator2])

These terms are combined into a single coupling term with frequency (coupling strength) couplingFreq. If you wish to set individual frequencies for the coupling terms, you can specify them with frequency when creating the terms.

[4]:

# coupling part of the JC-Hamiltonian
# first term of the coupling
JCCouplingT1 = JCSystem.createTerm(
    qSystem=[CavityJC, QubitJC],
    operator=[qg.destroy, qg.sigmap])

# second term of the couplings
# first term of the coupling
JCCouplingT2 = JCSystem.createTerm(
    qSystem=[QubitJC, CavityJC],
    operator=[qg.sigmam, qg.create])

# combine the terms into a single term
JCCoupling = qg.QTerm(qSystem=JCSystem,
                      subSys=[JCCouplingT1, JCCouplingT2],
                      frequency=couplingFreq)

Here, we introduce a dissipatorObj object to include cavity decay in the Lindblad master equation for evolving the JC system.

The dissipator for cavity photon decay can be described as:

$ {D}[a]:nbsphinx-math:rho `= a :nbsphinx-math:rho a^:nbsphinx-math:dagger - :nbsphinx-math:frac{1}{2}` \left`( a^:nbsphinx-math:dagger a :nbsphinx-math:rho + :nbsphinx-math:rho a^:nbsphinx-math:dagger a :nbsphinx-math:right`) $

where \(\rho\) = density matrix of the system, \(a\) = cavity annihilation operator.

The system evolves under both the JC Hamiltonian and the dissipative channel. This captures the irreversible loss of photons from the cavity captures via leakage to the environment.

Here, we set the decay rate \(\kappa\) = 0.1.

[5]:

# add a dissipator term for the cavity decay
cavityDissispator = qg.dissipatorObj(system=CavityJC, jOper=qg.destroy(cavityDim), jRate=.1)

# add the dissipation to the evolution of the JC system
cavityDissispator.addToProtocol(JCSystem._freeEvol)

[ ]:

# set the simulation parameters
JCSystem.simStepSize = timeStep
JCSystem.simTotalTime = totalTime

# initial state should be defined as a density matrix due to the altered evolution with dissipation.
# Otherwise, Quanguru will automatically convert it to a density matrix.
JCSystem.initialState = qg.densityMatrix(qg.tensorProd(qg.basis(QubitJC.dimension, 1,), qg.basis(CavityJC.dimension, 1)))

# create a frequency sweep for the cavity
freqSweep = JCSystem.simulation.Sweep.createSweep(
    system=CavityJC,
    sweepKey="frequency",
    sweepMax=qubitFreq+cavityFreq,
    sweepMin=qubitFreq-cavityFreq,
    sweepStep=.025)

We can create a composite operator in the extended Hilbert space so that measurement can be made without needing to project the state in the compute function.

Here, we create the composite qubit Pauli-Z operator \(\langle \sigma_{z} \rangle\) acting on the full qubit \(\otimes\) cavity space. This gives +1 for \(\vert e \rangle\) and -1 for \(\vert g \rangle\).

[7]:

# create sigma Z operator in the extended Hilbert space
compositeSZ = qg.compositeOp(qg.sigmaz(),   dimA=cavityDim)

# calculate the desired results and store
def compute(sim, args):
    stateJC = args[0]
    sim.resultsDict['zJC'].append(qg.expectation(compositeSZ, stateJC))

JCSystem.simCompute = compute

IMPORTANT NOTE FOR WINDOWS USERS : MULTI-PROCESSING (p=True) DOES NOT WORK WITH NOTEBOOK

You can use a python script, but you will need to make sure that the critical parts of the code are under if __name__ == "__main__": We are going to add further tutorials for this later.

[8]:

# do not store the states
JCSystem.simDelStates = True

# run the simulation
# p=True uses multi-processing for the sweep
JCSystem.runSimulation(p=(platform.system() != 'Windows'))

Starting sequential sweep with 81 parameter combinations...
Simulation Start:       2025-11-28 15:56:26
[████████████████████████████████████████] 100.0% | Est. Finish: 2025-11-28 15:57:01
Simulation completed in 00h 00m 35s

[8]:

[]

[9]:

# Plot sigma Z expectation value as a function of time and cavity frequency
Y, X = np.meshgrid(JCSystem.simulation.timeList, freqSweep.sweepList)
plt.pcolormesh(X, Y, JCSystem.simulation.resultsDict['zJC'])
plt.colorbar()
plt.xlabel(r"Frequency $\omega_{c}$")
plt.ylabel("Time")

[9]:

Text(0, 0.5, 'Time')
../../../_images/classes_Tutorials_1_Qubit_18_JaynesCummingsModelOpenSystemDynamics_14_1.png

The plot shows \(\langle \sigma_z \rangle\) as a function of time and cavity frequency for the open Jaynes-Cummings system with cavity decay. Compared to the closed system in Tutorial 17, we observe:

  • Damped oscillations: The Rabi oscillations decay over time due to photon loss from the cavity at rate \(\kappa = 0.1\)

  • Asymptotic behavior: \(\langle \sigma_z \rangle\) approaches -1 (ground state) as the cavity photons are lost to the environment

  • Decoherence: The system loses quantum coherence, with a fading chevron pattern

  • Decay timescale: The oscillations decay on a timescale \(\sim 1/\kappa = 10\), competing with the Rabi period \(\sim \pi/g \approx 12.6\)

The dissipator \(\mathcal{D}[a]\rho = a \rho a^\dagger - \frac{1}{2}(a^\dagger a \rho + \rho a^\dagger a)\) describes the irreversible loss of photons from the cavity, modeling realistic experimental conditions where cavities have finite quality factors.

This Page