Quantum simulation of open quantum systems in heavy-ion collisions

We present a framework to simulate the dynamics of hard probes such as heavy quarks or jets in a hot, strongly-coupled quark-gluon plasma (QGP) on a quantum computer. Hard probes in the QGP can be treated as open quantum systems governed in the Markovian limit by the Lindblad equation. However, due to large computational costs, most current phenomenological calculations of hard probes evolving in the QGP use semiclassical approximations of the quantum evolution. Quantum computation can mitigate these costs, and offers the potential for a fully quantum treatment with exponential speedup over classical techniques. We report a simplified demonstration of our framework on IBM Q quantum devices, and apply the Random Identity Insertion Method (RIIM) to account for CNOT depolarization noise, in addition to measurement error mitigation. Our work demonstrates the feasibility of simulating open quantum systems on current and near-term quantum devices, which is of broad relevance to applications in nuclear physics, quantum information, and other fields.

Introduction. Considerable advancements in quantum devices, such as qubit coherence times, have recently been achieved [1][2][3][4]. Together with parallel progress in quantum algorithms and executable quantum software, nontrivial quantum computations can be carried out, including hybrid quantum-classical algorithms such as the variational quantum eigensolver [5][6][7][8][9][10] and fully quantum simulations of the unitary time evolution of closed quantum systems [11,12]. In high energy and nuclear physics, a variety of quantum computing applications have emerged . In particular, quantum simulation can be applied to study dynamics of large size systems that are in principle intractable with classical methods. To perform such simulations, quantum circuits compiled into single-and multi-qubit gates can be implemented on digital quantum computers.
Many physical systems of interest are not closed, but consist of a subsystem interacting with an environment. The dynamics of the subsystem can be formulated as an open quantum system. In the Markovian limit (in which the environment correlation time is much smaller than the subsystem relaxation time), the evolution of the subsystem is governed by a generalization of the Schrödinger equation known as the Lindblad equation [38][39][40], where instead of keeping track of all of the environmental degrees of freedom, one only needs to record environment correlators that are relevant for the subsystem evolution.
A key challenge in extending quantum simulation to open quantum systems is that the Lindblad evolution is nonunitary. During the last decade, algorithms have been developed to overcome this issue, most of which couple the subsystem with auxiliary qubits (whose dimension can be significantly smaller than that of the environment) such that the whole system evolves unitarily [41][42][43][44][45][46][47]. More re-cently, simulations of open quantum systems have been carried out on real quantum devices, but without error mitigation [48].
The evolution of hard probes in the QGP can be treated as an open system evolving in a hot medium. A fully field-theoretical description of hard probes in the medium is challenging and typically various approximations are made. Most studies employ semiclassical Boltzmann or Fokker-Planck (equivalent to Langevin) equations [63][64][65][66][67][68][69][70]; semiclassical transport equations are leading order terms in the gradient expansion of the Wigner transformed Lindblad equation [71,72]. Recently, several studies have applied Lindblad equations directly to investigate quarkonia [73][74][75][76][77][78][79][80] and jets [81,82], which are valid if the subsystem and environment are weakly coupled. It is expected that as the size of the subsystem increases (such as the jet radiation phase space, or the number of heavy quarks [83,84] in the subsystem), solving Lindblad equations would challenge the limits of classical computation. Quantum computing offers a possibility to remove the constraint on the subsystem size, and go beyond the approximations made in semiclassical approaches. Moreover, quantum simulation may provide a solution to the notoriously difficult sign problem in classical lattice QCD calculations of real time observables [14,[85][86][87] (the same problem can also appear in open QCD systems).
In this letter, we outline a formulation of the evolution of hard probes in the QGP as a Lindblad equation and explore how simulations on Noisy Intermediate Scale Quantum (NISQ [13]) devices can be used to advance theoretical studies of hard probes in the QGP. Using a quantum algorithm for simulating the Lindblad equation, we study a toy model on IBM Q simulators and quantum devices, and implement error mitigation for measurement and two-qubit gate noise. We demonstrate that quantum algorithms simulating simple Lindblad evolution are tractable on current and near-term devices, in terms of available number of qubits, gate depth, and error rates.
Open quantum system formulation of hard probes in heavy-ion collisions. The Hamiltonian of the full system consisting of the hard probe (subsystem) and the QGP (environment) can be written as Here H S , H E and H I are the Hamiltonians of the subsystem, the environment and their interaction, respectively. A schematic diagram of the setup is shown in Fig. 1. We further split H S into the free H S0 and the interacting part of the subsystem H S1 . In quantum field theories, Hamiltonians are functionals of fields, which require discretization in position space [16]. Here, instead of simulating the dynamics of fields, we focus on simulating the dynamics of particle states, which is valid for hard probes. If we use multi-particle states |p 1 , A 1 ⊗ · · · ⊗ |p n , A n as the basis where p i is the four-momentum, A i represents all discrete quantum numbers, and i = 1, 2, . . . , n, then both H S0 and H S1 are matrices and H S0 is diagonal. Note that H S1 is different from H I : The former is the interaction within the subsystem itself and independent of the environment, while the latter represents the interaction between the subsystem and the environment. For example, for jets in HICs, H S1 can be collinear radiation of collinear particles while H I can describe the Glauber exchange between collinear particles (subsystem) and soft fields from the QGP environment [81]. The total density matrix of the subsystem and the environment evolves under the von Neumann equation. In the interaction picture, this is given by The operators are defined by The levels in S can represent for example: (1) heavy quarkantiquark (QQ) bound states |p, Ai with center-of-mass momentum p and quantum numbers Ai, and (2) unbound QQ pairs |p 1 , p 2 with momenta p 1 , p 2 . For jets the levels of S can represent multi-parton states labeled by momenta |p1, · · · , pn .
The interaction picture used here is special: it is the standard interaction picture for the subsystem but it is the Heisenberg picture for the environment. We will drop the superscript (int) from now on for simplicity but the reader should be reminded that we use the interaction picture throughout. We assume that the initial density matrix factorizes and the environment density matrix is a thermal state 1 where β = 1/T is the inverse of the QGP temperature. After the environment is traced out, the reduced evolution of the subsystem density matrix is generally timeirreversible and non-unitary. If the coupling between the subsystem and the environment is weak, the reduced evolution equation can be cast as a Markovian Lindblad equation [38][39][40]: where H L denotes a thermal correction to H S generated by loop effects of H I , and the L j are called Lindblad operators, whose explicit expressions will be given for a toy 1 The backreaction of the QGP medium to jet energy loss [88][89][90][91][92][93][94][95][96][97], which may further modify jet observables is beyond the scope of our considerations here. For a recent review, see Ref. [98].  [99][100][101][102][103] or the AdS/CFT correspondence [104][105][106][107][108]. For the nonperturbative computation, one needs to formulate the theory such that the relevant correlator is gauge invariant, where effective field theory can be used. A concrete construction of gauge invariant correlators for quarkonium transport can be found in Refs. [72,109]. Quantum algorithm. We will apply a quantum algorithm based on the Stinespring dilation theorem, see for example Refs. [44,110], to simulate the Lindblad equation. The algorithm in terms of the evolution operators J, defined below, and H S , is illustrated in Fig. 2. The algorithm couples the subsystem with auxiliary qubits, which are traced out after each time step ∆t. The dimension of the auxiliary register is m + 1 and the number of qubits needed in practice for the register is ceil(log 2 (m + 1)) ≡ 2 log 2 d . Together with the number of qubits required to record the subsystem state, the total number of qubits needed is 3 log 2 d . We use {|0 a , |1 a · · · , |m a } to label the basis of the auxiliary register, indicated by the subscript a.
We assume the initial state ρ S (0) = |ψ S (0) ψ S (0)| is a pure state 2 . At the beginning of each cycle at time t, the total density matrix of the subsystem and the auxiliary is set to be a (m + 1) × (m + 1) block matrix The J-operator is also a (m + 1) × (m + 1) block matrix where each block is a d × d matrix. One can show that the circuit in Fig. 2 reproduces (9) when ∆t → 0. To simulate the evolution from 0 to t, the size of the time steps is ∆t = t/N cycle where N cycle is the number of cycles, see Fig. 2.
Toy model and simulation on IBM Q. Simulating real jets and heavy quarks on quantum devices requires a large number of fault-tolerant qubits. As a proof of concept, we consider the following toy model that includes qualitative features of hard probes: where we use X, Z to denote the single qubit Pauli gates (Pauli matrices). The subsystem Hamiltonian H S is a two level system with energy difference ∆E. The two levels can correspond to the bound and unbound state of a heavy quark-antiquark pair, exchanging energy with QGP. The environment H E is a 3 + 1D scalar field theory, that together with (8) mimics the thermal QGP.
Here Π is the canonical momentum conjugate to φ. The extension to gauge theories requires a gauge invariant formulation of the environment correlator as mentioned earlier. The environment correlator can be calculated nonperturbatively to all orders in λ. Here for simplicity, we set m = λ = 0. Nonvanishing m and λ lead to different coefficients of the Lindblad operators but do not alter the quantum algorithm. The interaction strength g between the subsystem and the environment is unitless. In the Markovian limit, two Lindblad operators j = 0, 1 are relevant: where Γ 0 = g 2 ∆En B (∆E)/(2π), Γ 1 = g 2 ∆E/(2π) + Γ 0 and n B (∆E) = 1/(exp(β∆E) − 1) is the Bose-Einstein distribution. We will neglect H L in this letter. For our numerical studies, we use a unit system where all quantities are counted in units of T , the temperature of the medium. We initialize the state as ρ S (t = 0) = |0 0| and choose ∆E = 1(T ). The result for this toy model obtained from the IBM Q qiskit simulator [111] is shown in Fig. 3. We measure P 0 (t) ≡ 0|ρ S (t)|0 , which can be interpreted as the timedependent nuclear modification factor. Each time point corresponds to an independent quantum circuit, where the measurement is performed only at the end, as shown in Fig. 2. The results of the quantum algorithm with N cycle = 100 are shown for different values of the coupling g. They are consistent with the results obtained with a 4th order Runge-Kutta method that solves Eq. (9) classically. This agreement demonstrates that the circuit successfully solves the Lindblad equation. As expected, the strength of the coupling g controls the rate of approaching thermalization.
In order to run the circuit on a quantum device, we select N cycle = 1 in order to achieve a sufficiently small circuit depth. Modern quantum software packages are available to compile quantum circuits that approximate general unitary operators with minimal error and optimal depth [10,[112][113][114]. We synthesize a circuit for the e −iJ √ ∆t operator in terms of single qubit and cnot gates using the qsearch compiler [114]. The compiler yields circuits with 70 gates on average, including approximately 10 cnots per cycle; an example circuit for one cycle is shown in the supplemental material.
The results obtained from IBM Q Vigo device [115] are shown in Fig. 4. In addition to the uncorrected result, the results with readout and cnot error mitigation are also shown. We correct the readout error using the constrained matrix inversion approach in IBM's qiskit-ignis package. The response matrix can be found in the supplemental material. We also correct for cnot noise using a leading order zero-noise extrapolation based on the recently developed resource efficient Random Identity Insertion Method (RIIM) [116]. This procedure corrects for depolarization noise using a set of additional (cnot) 2 identity insertions, at the expense of amplifying statistical noise. Each data point corresponds to 5 evenly spaced time points that are averaged together. Each time point is calculated from the average of 49152 shots (runs). We observe that the error mitigation is more important at small values of t. Similar results were reproduced on the IBM Q Valencia and Santiago devices [117,118]. Overall, we observe good agreement of the results from the quantum device with the results from the simulator for N cycle = 1 after the error mitigation is applied. The choice of N cycle = 1 is seen to be a reasonable approximation for sufficiently small t. Moreover, a modest increase to N cycle = 3, as shown by the simulator in Fig. 4, yields considerably improved convergence, which is promising for near-term applications. These results demonstrate that the simulation of open quantum system dynamics relevant for HICs should be feasible on current and nearterm quantum devices.
Conclusions and Outlook. We performed simulations of open quantum systems using quantum devices from IBM Q. In particular, we focused on simulating the nonunitary evolution of a subsystem governed by the Lindblad equation. We demonstrated that digital quantum simulations with a few qubits and a circuit depth of ∼ 70 gate operations with ∼ 10 cnot gates are feasible on current quantum devices. We used the qsearch compiler to construct the quantum circuit, and implemented twoqubit gate error mitigation using zero noise extrapolation with the Random Identity Insertion Method (RIIM), in addition to readout error mitigation. Simulating open quantum systems is of great importance for theoretical studies of hard probes in heavy-ion collisions. The open quantum system formulation allows one to go beyond semiclassical transport calculations currently used in most phenomenological studies. Future calculations, using a time dependent environment density matrix may allow one to explore a broad range of physical models by varying medium properties such as the initial temperature, microscopic structure, or the probe-medium coupling. Open quantum systems are also relevant for various other systems in nuclear and high-energy physics such as studies of Cold Nuclear Matter effects at the future Electron-Ion Collider [119], the resummation of large logarithms relevant for jet physics [120][121][122][123] and studies of the Color Glass Condensate [124,125].  [115] which is used for the readout error mitigation in Fig. 4. The 2 3 states are prepared by applying X gates and then corresponding measurements are performed. The error mitigation is implemented using the constrained matrix inversion approach which is implemented in IBM's qiskit-ignis package [111].