Singly-excited resonant open quantum system Tavis-Cummings model with quantum circuit mapping

Tavis-Cummings (TC) cavity quantum electrodynamical effects, describing the interaction of N atoms with an optical resonator, are at the core of atomic, optical and solid state physics. The full numerical simulation of TC dynamics scales exponentially with the number of atoms. By restricting the open quantum system to a single excitation, typical of experimental realizations in quantum optics, we analytically solve the TC model with an arbitrary number of atoms with linear complexity. This solution allows us to devise the Quantum Mapping Algorithm of Resonator Interaction with N Atoms (Q-MARINA), an intuitive TC mapping to a quantum circuit with linear space and time scaling, whose N+1 qubits represent atoms and a lossy cavity, while the dynamics is encoded through 2N entangling gates. Finally, we benchmark the robustness of the algorithm on a quantum simulator and superconducting quantum processors against the quantum master equation solution on a classical computer.


Introduction
The Tavis-Cummings (TC) model [1], which describes interaction of N atoms with an optical cavity has been a cornerstone in the studies of quantum optical systems.The collective interactions in this model give an

√
N increase in the light-matter interaction rate (Fig. 1) and a host of subradiant states with rich phenomenology relevant for the development of quantum networks [2][3][4], all-photonic quantum simulators [5], quantum memories [6,7], quantum transport [8], exciton-polarons in semiconductors [9], superconducting quantum circuits [10], collective interaction of the cavity mode with an ensemble of atoms [11][12][13][14], and entanglement generation [15][16][17][18][19][20].Rapid progress in experimental development in the field of nanophotonics, renders the impracticality and scarceness of theoretical approaches unsatisfactory, especially in the open quantum system setting where the cavity interacts with the environment.Although recent results demonstrate that generalized TC model is integrable and can be solved using a variant of Quantum Inverse Methods (QIM) [21,22], solutions obtained in this way poise difficulties in extracting physical quantities and capturing dynamical correlations in the system.On the other hand, numerical solutions obtained through the quantum master equation [23] are limited by the exponential runtime complexity in Hilbert space size, and have thus far been performed for a single digit number of atoms.Due to the impracticality of analytical approaches based on QIM and exponentially rising cost of numerical solutions of the quantum master equations for such systems, theoretical verifications of experimental results are constrained to low number of atoms.Increasing the size of the Hilbert space has been pursued via approximate methods with polynomial scaling, such as the effective Hamiltonian [24], scattering matrix [4] and quantum trajectories [25] approaches.Furthermore, for applications concerned with the singly-excited regime, exact methods can be derived under linear scaling.
The availability of the Noisy Intermediate Scale Quantum (NISQ) devices has attracted interest for simulating open quantum systems.To date, two prevailing directions have emerged, the first using operator sum representation, where Sz.-Nagy theorem is used to relate Kraus operators with unitary dilatation matrices [26,27] that can then be directly implemented on a quantum circuit.This result has been further generalised and applied to quantum simulate the complex open quantum system, governed by the Fenna-Matthews-Olson Dynamics modelling the quantum theory of electron transfer in biological systems [28].An alternative approach is starting directly from the equations of motion in Lindblad and Gorini-Kossakowski-Sudarshan-Lindblad form and mapping the dynamics to a quantum circuit, which has been applied so far to both Markovian and non-Markovian open quantum systems consisting of 1 or 2 qubits [29].This approach has recently been verified on a canonical model of light matter interaction systems: the Jaynes-Cummings model [29].Cavity quantum electrodynamical models that involve multiple emitters, such as the TC model, have not yet been considered.However, this in particular is the area where classical methods quickly saturate numerical resources and quantum devices may be able to expand the Hilbert size of systems studied in quantum communication, memories and simulators.Moreover, studying a quantum system on purely quantum hardware may provide representations that are intuitive in nature, as both emitters and qubits are two-level systems.
In this work, we study a resonant open Tavis-Cummings model with arbitrary number of atoms and first provide an analytical solution for the singly-excited system with linear complexity.We then design the Quantum Mapping Algorithm of Resonator Interaction with N Atoms (Q-MARINA) which maps the open TC system with N atoms to a gate-based quantum circuit with only N + 1 qubits.We simulate the system on a superconducting quantum computer available through IBM Quantum program [30].

The model
We consider N two-level systems, modeling an ensemble of atoms (or spins), coupling to the environment of discrete bosonic modes.The system and the environment Hamiltonians H S and H E are: while their interaction is described by the Tavis-Cummings Hamiltonian H I : Here, we use the collective system operators To solve the model analytically, we aim to obtain the time dependent Hamiltonian in the interaction picture in the form of: Here, the form of b k (t) = b k e −iω k t is easily derived, however, finding an elegant expression for S ± (t) requires closer consideration.
We first note that, in the Hilbert space of the system, the operator is diagonal in terms we will call x p , 1 ≤ p ≤ 2 N .We find that, x p is a function of the Hamming weight W(p − 1), i.e. the digit sum of the binary representation of the number p − 1, as: Therefore, the term e iH S t = diag({e iω s tx p }) must too be diagonal, which allows us to obtain a closed form solution: thus completing the Eq. ( 3) for the time-dependent interaction Hamiltonian of the TC model with N identical two-level atoms.

Reduced density matrix
The corresponding reduced density matrix ρ S (t) of the TC system with N atoms is given by where c s n are the wavefunction coefficients with the following dependence on the cavity loss and cavity-atom interaction parameters: Real and positive coefficients c s n (0) are subject to the normalization constraint N n=1 ∥c s n (0)∥ 2 = 1, and full derivation of the density matrix is given in section 4.1.One of the key considerations to arrive to an exact solution and study the dynamics of the open TC system with N identical two-level atoms is that the total number of excitations in our system is a constant of motion [H, M] = 0, where H = H I + H S + H T , and is the total number of excitations in the system considered here.It has been shown that if one exploits the permutational symmetry [31] originating from the simplification that we are considering N identical emitters, one can gain further insights into closed systems beyond a single excitation manifold both analytically and numerically [32][33][34][35][36].In the case of the open Tavis-Cummings model studied here, this is seen in the symmetry of our solution for the wavefunction coefficients obtained in Eq. ( 7), where the identical choice of initial conditions would lead to identical behavior of subradiant states.
In the following, we show that this dynamics can be mapped onto a quantum circuit with N + 1 qubit, thus enabling quantum modeling of the Tavis-Cummings open quantum system on a gate-based quantum computer.While the solution derived in this section is a general one for a single-excitation system, for simplicity, we will from now on assume that the first emitter in the system is the one that is initially excited, while others are in the ground state (c S 1 (0) = 1, c S m (0) = 0 for m = 2, . . ., N), and the proposed quantum circuit will reflect that.

Quantum circuit
Here, we devise the Quantum Mapping Algorithm of Resonator Interaction with N Atoms (Q-MARINA), an (N+1)-qubit quantum circuit that evolves an open quantum system of N atoms and a resonant cavity in the single-excitation regime.The quantum circuit consists of N system qubits Q S n and one environment qubit Q E .The initial state is the excited state of one of the atoms, here Q S 1 which is subject to an Xgate.Subsequent application of CU 3 and CNOT gates between Q S 1 and Q E entangles the first atom and the environment, and then N − 1 sequences of CU 3 and CNOT entangling gates are applied to each of the qubits Q S 2 ,...,Q S N paired with Q E , in the opposite direction than for the Q S 1 .The corresponding quantum circuit is shown in Fig. 2. Here, the parameters of the CU 3 gates, CU 3 (2θ n )=CU 3 (2θ n , 0, 0) are selected to implement the Lorentzian density of states of the cavity open to the environment into the circuit:  θ n = arcsin thus resulting in excited state measurement probabilities of the system qubits Q S n equal to ∥c s n (t)∥ 2 .Importantly, this quantum circuit maintains the physical connections typical of the TC model where each atom directly interacts only with the cavity, reflected in entangling gates operating solely on system-environment qubit pairs.

Implementation of the Quantum Mapping Algorithm on superconducting circuits
As a testbed for quantum simulations of the lossy TC model, we implement the devised Q-MARINA quantum algorithm on the IBM Q Experience hardware.We first demonstrate agreement of the results for open system dynamics obtained through the implementation of the quantum circuit on the IBM QASM simulator provided via Qiskit library [39] .The comparison of the numerical solution of Quantum Master Equation (QME) for N = 7 atoms with the execution of the Q-MARINA quantum circuit in QASM simulator is illustrated in Fig. 3.
We then execute the proposed quantum circuit on the superconducting quantum devices ibmq quito (Falcon r4T processor) and ibm oslo (Falcon r5.11H processor), available through the IBM Quantum program.The quantum circuit requires star connectivity as all system qubits Q S n interact with the environment qubit Q E , therefore we selected devices that can support that layout in a 3-and 4-qubit circuits within the computers' heavy-hexagon topology.The comparison of our quantum device results with the previously obtained benchmarks on the QASM simulator and numerical QME solutions are shown in Fig. 4. The demonstrated close agreement between the solutions of the QME with Q-MARINA executed on QASM simulator and IBM Q quantum devices indicates that NISQ era quantum computers can be used to simulate open quantum system dynamics of highly dimensional models.

Discussion
In this work, we have explored quantum circuit mapping of the dynamics of N two-level atoms in a a lossy optical cavity.By restricting the open quantum system to a single excitation, typical of experimental realizations in quantum optics, we have analytically solved the TC model with an arbitrary number of atoms achieving reduced modeling complexity.This solution enabled us to devise the Quantum Mapping Algorithm of Resonator Interaction with N Atoms (Q-MARINA), an intuitive TC mapping to a quantum circuit with linear space and time scaling.We note here that this work does not aim at quantum advantage, but rather to show that the studied regime of Tavis-Cumming physics in a lossy resonator can be efficiently mapped to N + 1 qubit, as opposed to an infinite number of qubits.
It is interesting to note that the execution of the Q-MARINA quantum circuit illustrated in Fig. 2 on the the superconducting quantum devices ibmq quito and ibm oslo are in good agreement with the numerical solution of the QME (c.f.Fig 4), despite the fact that no error mitigation technique has been considered thus far.These results demonstrate that the open quantum system Tavis-Cummings physics can be simulated on the existing quantum hardware with an intuitive mapping between atoms and qubits and a substantial reduction in complexity implemented through the entangling gates with a single environment qubit.That being said, we acknowledge multiple challenges on the hardware side that need to be resolved before achieving e.g.coherence stability of the quantum devices with the number of qubits comparable to the number of atoms where classical solutions of the master equation become intractable.Therefore, a numerical solution of QME [38], as well as analytical approaches such as mean-field approximation [40], or Keldysh's action formalism [41,42] remain valuable go-to methods for studying the complex dynamics of quantum fluctuations in the TC-like systems.
The devised mapping of the TC system with N identical atoms constitutes a first step toward using superconducting NISQ processors to design new optical quantum devices.The results obtained on existing quantum devices are further limited by the quantum computer size and the corresponding topology which provides the desired star-connectivity to up to 4 qubits.Alternative quantum platforms which provide allto-all connectivity, such as those based on trapped ions [43,44] or atoms [45], may provide options to scale the problem size by at least an order of magnitude [46,47].Once the number of qubits is scaled, the number of entangling gates relative to the qubit coherence time will be the measure of the performance of our algorithm, as the circuit depth scales linearly with the number of atoms.

Reduced density matrix derivation
The wavefunction of an N-atom Tavis-Cummings system in the low-excitation regime is given by the superposition of the vacuum state |g0⟩, single excitations of the n-th atom |e n 0⟩ and the single excitations of the k-th bosonic mode The Schrödinger equation with Hamiltonian given in Eq. (3) yields a system of differential equations: We next note that The system of differential equations transforms to It follows that the k-th cavity mode and the n-th atom amplitude can be expressed as where we approximate the environment coupling terms with Lorentzian density of states modeling the cavity dynamics k ∥g k ∥ 2 = dωJ(ω).The term describes an optical resonator with loss rate κ coupled to an atom at interaction rate g, and represents the channel through which the system interacts with the environment.For a closed system, the cavity would respond to only a singular frequency (κ = 0).Next, similarly to [29], we define and the atomic amplitudes simplify to Taking the Laplace transform of l.h.s. and r.h.s. of the previous equation, we obtain: where cs n (s) and f (s) denote the Laplace transforms of the functions c s n (t) and f (t−t ′ ) defined in Eq. ( 19).Solving the system of coupled equations given in Eq. ( 21) for cs n (s) and performing an inverse Laplace transform gives us wavefunction coefficients which determine the density matrix: where To obtain the reduced density matrix ρ S (t) that describes the state of the system, we remove the environment degrees of freedom through a partial trace: From here, we express the diagonal elements of the (N + 1)-dimensional density matrix as: where the first N diagonal elements correspond to the excited state measurement probabilities of the two-level atoms, represented in Fig. 2 by system qubits Q S n .

Quantum Mapping Algorithm Implementation Details
Here, we give further details on the implementation of the devised Q-MARINA quantum algorithm on the IBM Q Experience hardware.The comparison of the results for open system dynamics illustrated in Fig. 3 is obtained by implementing the quantum circuit on the IBM QASM simulator provided via Qiskit [39] library and contrasting it with the numerical solution of the Quantum Master Equation (QME) modeled in Quantum Toolbox in Python (QuTiP) [37,38] on a classical computer.The combination of the system parameters-loss rate κ and coupling constant g-determine whether the light-matter interaction is considered to be in the weak or or the strong coupling regime.Concretely, in our case with N atoms g √ N < κ/4 corresponds to the weak coupling strength, while for g √ N ≥ κ/4 we reach the strong coupling regime [24], particularly relevant for hybridization of light and matter explored in quantum light generation and extension of coherence in quantum memories.Thus, Fig. 3 compares the QME and the Q-MARINA QASM results for N=7 atoms in the strong coupling regime.
The Q-MARINA implementation on IBM Q hardware shown in Figure 4 studies 3-qubit and 4-qubit circuits on one-to-all connected subgraphs of ibmq quito (Falcon r4T processor) and ibm oslo (Falcon r5.11H processor), respectively, simulates the N = 2 TC system in strong coupling regime and the N = 3 TC system in the borderline regime where an individual atom couples weakly, while the collective coupling is in the strong regime of the cavity QED.The atomic amplitudes follow the exact QME solution closely and leave space for future precision improvement via error mitigation techniques.

Fig. 1 a
Fig. 1 a) An illustration of an open Tavis-Cummings system consisting of an optical cavity of the loss rate κ with N atoms each coupled at the interaction rate g.b) The transmission spectrum of an empty cavity (dashed gray line) featuring a lorentzian profile with linewidth κ and a cavity resonantly coupled to N atoms (solid orange line) featuring two polariton peaks separated by 2g √ N.

Fig. 2 Q
Fig.2Q-MARINA algorithm that maps an open quantum system of N two-level atoms in a lossy cavity to a quantum circuit with N + 1 qubits and 2N entangling gates that encode the interaction of atoms (Q S n ) with the cavity and environment (Q E ).

Fig. 3
Fig.3The evolution of the singly excited open quantum system Tavis-Cummings model of N = 7, g = κ = 5 calculated using a) quantum master equation in QuTiP software[37,38] and b) Q-MARINA algorithm in QASM simulator with 40,000 shots per data point.

Fig. 4 Q
Fig. 4 Q-MARINA simulation of the singly excited open TC system for evolving upon excitation of the Atom 1, executed on a) ibmq quito quantum computer with 10,000 shots per point for N = 2, g = 10, κ = 5, and b) ibm oslo quantum computer with 10,000 shots per point for N = 3, g = 2, κ = 5.The exact QME solution is plotted for comparison.