From integrability to chaos in quantum Liouvillians
FFrom integrability to chaos in quantum Liouvillians ´Alvaro Rubio-Garc´ıa, ∗ Rafael A. Molina, † and Jorge Dukelsky ‡ Instituto de Estructura de la Materia, IEM-CSIC, Serrano 123, Madrid, E-28006, Spain
The dynamics of open quantum systems can be described by a Liouvillian, which in the Markovian approx-imation fulfills the Lindblad master equation. We present a family of integrable many-body Liouvillians basedon Richardson-Gaudin models with a complex structure of the jump operators. Making use of this new re-gion of integrability, we study the transition to chaos in terms of a two-parameter Liouvillian. The transition ischaracterized by the spectral statistics of the complex eigenvalues of the Liouvillian operators using the nearestneighbor spacing distribution and by the ratios between eigenvalue distances.
I. INTRODUCTION
Quantum chaos has gained a renewed interest in the last fewyears due to its relevance to important current topics of re-search like quantum thermalization, non-equilibrium dynam-ics of quantum many-body systems and quantum information[1–5]. In the standard Hermitian quantum mechanics validfor closed systems, the dynamical regime of the system canbe characterized through spectral statistics [6]. Integrable sys-tems present level clustering and the nearest-neighbor spacingdistribution follows the one-dimensional Poisson distribution P ( s ) = e − s [7], while chaotic systems present level repul-sion with the P ( s ) close to the Wigner surmise of RandomMatrix Theory (RMT) depending on their symmetry class [8], P ( s ) ∝ s β , with β = 1 , , for orthogonal, unitary, andsymplectic symmetries respectively. The Bohigas-Gianoni-Schmit (BGS) conjecture [8] is well founded now in the semi-classical theory [9–11] and supported by overwhelming evi-dence in many different quantum systems, experimentally andnumerically [12–14]. The transition between the integrableand chaotic universal limits is non-universal, depending on thespecifics of the particular system under study, although it hasbeen explored in much detail for different systems [15, 16].In the transition between the integrable and chaotic orthogo-nal cases, for example, some systems present fractional levelrepulsion while others present full level repulsion but only fora fraction of levels [17].In open quantum systems, the theory is much less devel-oped, even if the first results came shortly after the proposalof the BGS conjecture [18]. Open quantum systems, can bedescribed by the Liouville equation, which characterizes thetime evolution of the density matrix operator. In the Marko-vian approximation, the Liouvillian is a linear non-Hermitianoperator, and the Liouville equation can be written as a Lind-blad master equation [19]. The Liouvillian, then, has complexeigenvalues instead of the real energies of standard Hermi-tian quantum mechanics. The initial approach to the prob-lem was to study integrable or chaotic Hamiltonians with aweak coupling to the environment. When the Hamiltonianwas integrable, Grobe et al. studied the spectral statistics in ∗ [email protected] † [email protected] ‡ [email protected] the complex plane and found good agreement with the two-dimensional Poisson distribution [18]. In the chaotic limitthere appears universal cubic repulsion P ( s ) ∝ s as in theGinibre ensemble of non-Hermitian random matrices [20],although the details of the P ( s ) distribution depend on thesymmetries of the non-Hermitian matrix [21, 22]. For anopen quantum spin chain, the level spacing distribution inthe transition from integrability to chaos has been fitted by astatic two-dimensional Coulomb gas with harmonic confine-ment where level repulsion is given by the inverse temperatureshowing fractional level repulsion in the transition [23].Very recently, the need for a different approach in the studyof the integrable and chaotic properties of open quantum sys-tems has been triggered as a consequence of the discovery ofnew families of integrable Liouvillians [24–26]. Extendingthe class of exactly solvable and quantum integrable Liouvil-lians is an important important step towards improving ourunderstanding of open quantum systems. The statistical prop-erties of the complex spectra of random chaotic Liouvillianshave been studied in a few recent works [27, 28]. However,the transition between the exactly solvable integrable limit andthe chaotic limit in physical Liouvillians remains mostly un-explored.In this letter we will extend the model of Ref. [25] basedon the SU(2) spin 1 Richardson model to a line of integrabil-ity within the rational Richardson-Gaudin (RG) class of inte-grable models. This new family of integrable Liouvillians hasa rich and complex structure of the jump operators and allowsfor a simple parametrization along the line of integrability. Wethen define in terms of a single parameter a Liouvillian thatinterpolates between integrability and a fully chaotic limit.With these model Liouvillians we characterize these transi-tions through the spectral statistics of the complex eigenval-ues of the Liouvillian operators and by the average propertiesof the complex spacing ratios. II. INTEGRABLE RICHARDSON-GAUDINLIOUVILLIANS
We consider a system described by a Hamiltonian H weakly coupled to a Markovian environment. The evolutionof its density matrix can be expressed in terms of the Lindblad a r X i v : . [ n li n . C D ] F e b master equation [19] ∂ t ρ = L ( ρ ) = − i [ H, ρ ]+ (cid:88) j (cid:18) L j ρL † j − ρL † j L j − L † j L j ρ (cid:19) , (1)where L is the Liouvillian superoperator acting on the spaceof density matrices. The first term on the right hand side is theHamiltonian commutator describing the unitary evolution andthe second term takes into account the interaction with theenvironment through the Lindblad jump operators L j . Thetime evolution of a Liouvillian eigenstate L ( ρ i ) = λ i ρ i isgiven by e L t ρ i = e λ i t ρ i and it can be shown that the realpart of the eigenvalues is always non positive, Re { λ i } ≤ [19, 29, 30]. Furthermore, the Liouvillian spectrum consistsof complex conjugate pairs ( ρ i , λ i ) , ( ρ † i , λ ∗ i ) , implying thata real eigenvalue corresponds to an hermitian density matrix.There is at least one state with zero Liouvillian eigenvalue, the steady state ρ SS , to which any initial state decays in the longtime limit lim t →∞ e L t ρ = ρ SS . As a consequence, the traceof any eigenstate with negative real part must be zero.Similar to Ref. [25], we study a chain of L spins subject toa non uniform magnetic field h i H = L (cid:88) i =1 h i S zi , (2)whose interaction with the environment is mediated by quan-tum jumps of the form L az = (cid:112) L (cid:88) i =1 z ai S zi , L ± = (cid:112) Γ ± L (cid:88) i =1 x i S ± i , (3)where a labels a set of arbitrary different jumps. As wewill see, the equivalence of the gain L + and loss L − op-erator strengths, Γ = Γ + = Γ − , is a necessary conditionfor quantum integrability. In the particular case of a singlecollective spin Hamiltonian, this condition could be relaxed[26]. Note that for this choice the jumps become hermitian, ( L az ) † = L az , ( L + ) † = L − , which immediately leads to atrivial steady state ρ SS = I . The noisy spin model of Ref.[25] is a particular case of the general RG integrable modelscorresponding to the limit x i = 1 , z ai = δ a, (see Eq. (12)below).The vector representation of the Lindblad master equation[30] doubles the Hilbert space H of dimension N by mappingthe density matrix ρ into a vector | ρ (cid:105) in the space H ⊗ H ofdimension N . The mapping transforms every spin operator S acting to the right or to the left of the density matrix into asuperoperator of the form Sρ → S ⊗ I | ρ (cid:105) ˆ= S | ρ (cid:105) ρS → I ⊗ S T | ρ (cid:105) ˆ= J | ρ (cid:105) (4)Thus, the operators acting to the right of the density matrixare mapped to a new set of spin operators J i acting on a dual space with the same dimension N of the original space. Thismapping allows the Lindblad master equation to be written ina matrix-vector form L| ρ (cid:105) with a Liouvillian matrix L of di-mension N whose eigenvalues can be computed by an exactdiagonalization for small size systems or by an appropriatemany-body approximation [31–34].The complete Liouvillian in the Lindblad form is written,after a π rotation around the y axis in the dual space, as L = − i (cid:88) i h i ( S zi + J zi ) − Γ n j (cid:88) a =1 (cid:88) ij z ai z ai ( S zi + J zi ) (cid:0) S zj + J zj (cid:1) − Γ (cid:88) ij x i x j (cid:2)(cid:0) S + i + J + i (cid:1) (cid:0) S − j + J − j (cid:1) + h.c. (cid:3) (5)The fact that the spin S i and the dual spin J i appear alwaysas a sum is a consequence of the equality of the gain and loss Γ + ≡ Γ − jumps operators. We can, therefore, define the totalspin operators K i = S i + J i for each lattice site i . In termsof these new total spin operators the Liouvillian acquires thesimpler expression L = − i (cid:88) i h i K zi − Γ (cid:88) i x i ( K i ) + (cid:88) i (cid:34) Γ x i − Γ (cid:88) a ( z ai ) (cid:35) ( K zi ) − Γ (cid:88) a,i (cid:54) = j z ai z ai K zi K zj − Γ (cid:88) i (cid:54) = j x i x j (cid:0) K + i K − j + K − i K + j (cid:1) . (6)We stress that the magnitude of the spins at each site ofthe chain s i is completely arbitrary, with the only constraintimposed by the mapping that it equals the magnitude of thedual spin s i = j i . Hence, the total spin at site k i can take thefollowing values k i = 0 , , · · · , s i .The Liouvillian (6) presents a set of L +1 weak symmetries[29], given by operators that commute with the total Liouvil-lian [ L , O ] = 0 . These weak symmetries are the total mag-netization operator K z = (cid:80) i K zi and the L SU(2) Casimiroperators at each site K i = ( K zi ) + 12 (cid:0) K + i K − i + K − i K + i (cid:1) = k i ( k i + 1) (7)Since all symmetry operators commute [ K z , K i ] = 0 , theLiouvillian matrix is separated into blocks labelled by the setof L quantum numbers k i and each one, in turn, divided intosub-blocks labelled by the total k z . Due to the symmetry ofthe Liouvillian, the occurrence of a singlet total spin k i =0 at site i effectively removes this site from the chain, sincethe action of any SU(2) generator on this state annihilates it.Therefore, the state with all singlet spins is the steady state ofthe Liouvillian with 0 eigenvalue.For the particular case of our Liouvillian (6) with an in-homogeneous magnetic field along the z axis, there is an-other weak symmetry that generates the complex conjugatepair spectrum, F = C e iπK x . The operator C is complex con-jugation and e iπK x performs a π rotation around the x axisreversing the spins along the y and z directions. This sym-metry commutes with the Liouvillian and with the L Casimiroperators K i but it does not commute with the K z symmetry.In spite of the fact that F cannot be added to the chain of weaksymmetries, it is useful to understand the structure of the spec-trum in each k z block. For every Liouvillian eigenstate | ρ i (cid:105) (invector form) with complex eigenvalue λ i , the eigenstate | F ρ i (cid:105) is also a Liouvillian eigenstate with conjugate eigenvalue L| F ρ i (cid:105) = L F | ρ i (cid:105) = F L| ρ i (cid:105) = F λ i | ρ i (cid:105) = λ ∗ i | F ρ i (cid:105) . (8)Moreover, | F ρ i (cid:105) state has opposite total angular momentum K z to | ρ i (cid:105)(cid:104) F ρ i | K z | F ρ i (cid:105) = (cid:104) ρ i | F † K z F | ρ i (cid:105) = (cid:104) ρ i | ( − K z ) | ρ i (cid:105) . (9)Thus, each Liouvillian block with momentum k z contains thecomplex conjugate eigenvalues of the corresponding blockwith momentum − k z . Blocks with k z = 0 are exceptional,since they contain all complex conjugate pairs. This couldbe explained by the fact that in the subspaces of k z = 0 thesymmetries F and K z do commute.The Liouvillian (6) can be derived from the integrals of mo-tion (IOM) of the Richardson-Gaudin models [35, 36] R i = K zi − G L (cid:88) j ( (cid:54) = i )=1 X ij (cid:0) K + i K − j + K + j K − i (cid:1) + Z ij K zi K zj . (10)The local spins k i are in principle arbitrary. A possible choicefor the matrices X and Z such that the L IOM commuteamong themselves, [ R i , R j ] = 0 , is [36] X ij = (cid:112) (1 − α ) + α η i (cid:113) (1 − α ) + α η j η i − η j Z ij = (1 − α ) + α η i η j η i − η j , (11)with α ∈ [0 , an interpolation variable within the region ofintegrability. Note that the limits α = 0 , assure the equal-ity of the two matrices X and Z and pertain to the rationalor XXX RG family [37]. The IOM with α = 0 lead to theconstant pairing Hamiltonian originally solved by Richard-son [38], while α = 1 leads to the recently studied separablepairing Hamiltonian [39]. For intermediate values of α thematrices X and Z are not equal, corresponding to the XXZor hyperbolic RG family. Moreover, the IOM have a set of L + 1 parameters ( G, η i ) that are arbitrary complex numbers.In what follows, we consider the η (cid:48) s as real parameters and G as pure imaginary, G = ig . While this choice makes theIOM no longer Hermitian, the integrability properties of thesemodels are kept intact. In terms of these IOM we can obtainthe set integrable Liouvillians (6) with the following linear combination L = − i L (cid:88) i =1 η i R i = − i L (cid:88) i η i K zi − g L (cid:88) i (cid:54) = j x i x j (cid:0) K + i K − j + K − i K + j (cid:1) − g L (cid:88) i (cid:54) = j [(1 − α ) + α η i η j ] K zi K zj (12)where x i = (cid:112) (1 − α ) + α η i .This family of integrable Liouvillians can be derived fromam open quantum system with Hamiltonian H = (cid:80) i η i S zi ,gain and loss jumps (3) L ± = (cid:114) g L (cid:88) i =1 (cid:113) (1 − α ) + α η i S ± i (13)and two dephasing jumps L z = (cid:114) (1 − α ) g L (cid:88) i =1 S zi , L z = (cid:114) α g L (cid:88) i =1 η i S zi (14)In order to break integrability we will, in addition, con-sider a set of n j random jumps that lead to a chaotic (non-integrable) Liouvillian. These jumps are given by L a ± = √ Γ L (cid:88) i =1 w ai S ± i , a = 1 , . . . , n j , (15)with ( w a , w a , . . . ) a set of random orthonormal vectors thatare also orthonormal to the integrable gain and loss jump vec-tor ( x , x , . . . ) . The chaotic Liouvillian is then L chaotic = − i (cid:88) i h i K zi + Γ n j (cid:88) a =1 (cid:88) i ( w ai ) (cid:104) ( K zi ) − ( K i ) (cid:105) − Γ n j (cid:88) a =1 (cid:88) i (cid:54) = j w ai w aj (cid:0) K + i K − j + K − i K + j (cid:1) , (16)In general, n j ∼ L/ shows the most chaotic spectralstatistics, while a very high or low n j usually show interme-diate behavior in the spectral statistics.The transition from integrability to chaos will be character-ized by two parameters ( α, β ) , that can interpolate betweenthe two rational integrable liouvillians and the fully chaoticlimit. L ( α, β ) = (1 − β ) L int ( α ) + β L chaotic , (17)with ( α, β ) = (0 , the rational constant model, (1 , theseparable model and ( α, the chaotic limit. -10 -8 -6 -4 -2 0 Re( λ ) -3-2-101 I m ( λ ) a ) (0 , -10 -8 -6 -4 -2 0 Re( λ ) − − − b ) (1 , -10 -8 -6 -4 -2 0 Re( λ ) − − − c ) (0 , s P ( s ) d ) L = 7 L = 9 L = 11 γ = 0 γ = 1 AI † s e ) s f ) FIG. 1. (a-c) Integrable and chaotic Liouvillian spectra for L = 6 , K z = 1 , k i = 1 in the rational constant (a), rational separable (b) andchaotic (c) Liouvillian limits. (d-f) Level spacing distribution P ( s ) of the integrable and chaotic Liouvillians for different sizes L = 7 , , and k i = 1 , k z = 1 with the same ordering as in (a-c). The dashed and dash-dotted lines show the γ = 0 and γ = 1 limits of the interpolationdistribution P ( s ; γ ) (19), respectively, while the dotted line shows the level statistics of the AI † universality class. III. LIOUVILLIANS SPECTRAL STATISTICS IN THETRANSITION FROM INTEGRABILITY TO CHAOS
In this work we study the spectrum of Liouvillians of spin1/2 chains. Because the L + 1 weak symmetries block diago-nalize the Liouvillian space, we will focus on the k i = 1 , k z =1 subspaces, as this represents the block with maximal dimen-sion, adequate for the study of spectral statistics. While thesubspace k z = 0 is bigger, in principle, than k z = 1 , its spec-trum is divided in half by the F symmetry and results into twoblocks of smaller dimension.Figs. 1(a,b) show some typical examples of the spectrumof the rational constant and separable models, respectively,for Liouvillians with L = 6 , while Fig. 1(c) shows the spec-trum of the chaotic Liouvillian limit with the same param-eters and n j = L/ random jumps. The parameters η i are chosen as η i = 1 / i/L + ω i , where the values of ω i are random variables coming from an uniform distribution ω i ∈ [ − /L, /L ] . This is a convenient choice for randomvalues of η i that approximately keeps the relative strength ofthe Hamiltonian and the jumps constant as we increase size.The ordered patterns of the two integrable spectra, as we willsee below, have random local distribution of the eigenvalues,while the chaotic spectrum presents level repulsion.To characterize the integrable to chaos Liouvillian transi- tion we measure the nearest neighbor level statistics. Beforethat, an unfolding procedure is done in order to rescale the lo-cal level density to unity, which results in adimensional quan-tities to be able to compare different Liouvillian systems to theuniversal results of Random Matrix Theory [6]. The unfold-ing procedure can be tricky and if it is not performed carefullymay give rise to misleading results [40]. We use a local un-folding procedure where the level distances are rescaled bythe local level density [18]. This local density is computed bycalculating the area with a number n of nearest neighbors toeach level. The nearest-neighbor spacings are then S i = (cid:114) nπd i,n d i, , (18)where d i,n is the distance in the complex plane between level i and its n -th neighbor. Then, the vector (cid:126)S is rescaled byits mean to obtain the final unfolded spacings s i . There is acertain arbitrariness in choosing the number n , but it has tobe large enough to not depend on small scale fluctuations andsmall enough to capture the local statistics. We have found n = 20 to be a good compromise between these two limits,the results being quite stable and robust in a wide range aroundthis value.For integrable Liouvillians, the level statistics follow aPoisson distribution in the plane P ( s ) = π s exp (cid:0) − π s (cid:1) ,while the level statistics of chaotic Liouvillians belong to the AI † universality class of non-Hermitian symmetric matrices[21, 22]. To observe the transition from quantum integrabilityto chaos we fit the level spacing distribution P ( s ) to a distri-bution governed by a parameter γP ( s ; γ ) = A ( γ ) s γ +1 e − B ( γ ) s , γ ∈ [0 , (19)with A ( γ ) , B ( γ ) such that (cid:82) ∞ P ( s ) ds = (cid:82) ∞ sP ( s ) ds = 1 .Similar to the Brody distribution for the transition betweenintegrability and chaos in standard Hermitian Hamiltonians[15], the limits of this distribution go from the integrable limit γ = 0 , which corresponds to the Poisson distribution in theplane, to the chaotic limit γ = 1 , corresponding to the Wignersurmise of the level spacing distribution of the Ginibre en-semble [20]. Figs. 1(d-f) shows the limit γ = 0 (dashed line)and γ = 1 (dash-dotted). While the distribution in the limit γ = 1 does not exactly correspond to the level statistics ofthe AI † universality class (dotted line in Fig. 1), it shows thesame level repulsion P ( s ) ∼ s and both distributions are suf-ficiently similar so that γ = 1 is a good characterization of thechaotic limit.The unfolded level statistics of the two rational and chaoticlimits are shown in Figs. 1(d-f) for Liouvillians with sizes L = 7 , , and Hilbert space dimension , and , respectively. For each of these sizes L we computea sample of several Liouvillian spectra so that there are in to-tal more than 700.000 eigenvalues to draw statistics from. Thetwo integrable limits show levels statistics that are very closeto the 2d Poisson distribution in the plane, indicating a randomdistribution of the local eigenvalues, while the chaotic limitshows similar level statistics to the AI † universality class.For each Liouvillian limit there is also a convergence to therespective 2d Poisson and AI † distributions from smaller tolarger L ’s.Fig. 2(a) shows the transition from integrable to chaoticlevel statistics characterized by the interpolation parameter γ ,which we fit to the level spacings using a Maximum Likeli-hood Estimation. Each value of the parameters ( α, β ) rep-resents the statistics of a sample of 20 Liouvillians of size L = 10 , with a Hilbert space dimension of eigenvalues,and n j = 5 random jumps. For values of β = 0 the tran-sition parameter γ is equal to 0, showing integrable statisticsfor both rational integrable Liouvillians and for β > thereis a smooth transition from integrable to chaotic statistics upto the limit β = 1 , with the integrability breaking at relativelylow values of β .Apart from the interpolation distribution P ( s ; γ ) we studythe integrable to chaos transition using the complex spacingratios recently introduced by S´a et al. [41] z k = λ NN k − λ k λ NNN k − λ k (20)with λ k the k -th eigenvalue and λ NN k , λ NNN k its first and secondneighbors, respectively. The advantage of the spacing ratios isthat, because they are adimensional quantities, the unfoldingprocedure is not needed, at least if the local density of states does not vary in the scale of the typical minimal distance be-tween levels. The equivalent quantity for closed systems in-troduced in Ref. [42] has become widely used for this reason.The complex spacing ratios have only been computed for afew examples [28, 41, 43] yet so it is particularly importantto compare their results to the more standard spacing distribu-tions in the integrable to chaotic transition in Liouvillians.Figs. 2(b,c) show the means of the absolute and phase an-gle values, respectively, of the complex ratios z k = r k e θ k ,which were determined by S´a et al. to be good indicatorsof the transition. They find in the limit L → ∞ , the val-ues ( (cid:104) r (cid:105) , −(cid:104) cos θ (cid:105) ) = (2 / , for integrable statistics and ≈ (0 . , . for the AI † universality class. Our results for (cid:104) r (cid:105) [Fig. 2(b)] coincide with these limits for both integrable β = 0 and chaotic β = 1 families of Liouvillians, while −(cid:104) cos θ (cid:105) has a slow convergence to the infinite size limit, alsomentioned by S´a et al. and its limits are a bit lower than thoseof L → ∞ . Apart from the finite size effects, both measure-ments of the complex spacing rations show a good agreementwith the P ( s ; γ ) interpolation distribution. IV. CONCLUSIONS
To summarize, we have presented a new continuous fam-ily of integrable many-body Liouvillians based on the rationaland hyperbolic Richardson-Gaudin models. In vector form,these Liouvillians can be written as a fully coupled networkof spins of the XXX or XXZ type, representing open quantumspin systems with gain, loss, and dephasing jumps. This verygeneral family of integrable Liouvillians increases greatly thenumber of systems that can be investigated with exact solu-tions, which can open new avenues of research for the studyof the properties of open quantum many-body systems, a topicthat is still in its infancy.By adding chaotic jumps, we define two-parameter Liou-villians that describe the transition between the different in-tegrable models and chaotic non-integrable Liouvillians withthe same degrees of freedom. We have characterized thesetransitions studying the spectral statistics of the Liouvillians’seigenvalues in the complex plane using two methods, thenearest-neighbor spacing distribution and the average proper-ties of the complex spacing ratios. The complex ratios do notrequire the difficult unfolding procedure, which is an impor-tant advantage as in the case of complex spectra the unfoldingis even more critical than for real spectra.Within the integrable line, we find excellent agreement be-tween the spacing distribution and the distribution of a Pois-son process in the plane, while in the chaotic case we find avery good agreement with the AI † random matrix universalityclass. The agreement improves as we reach larger sizes. Wehave derived a simple interpolation formula between the linearand cubic repulsion cases as a function of one parameter. Theformula captures nicely the transition between integrable andchaotic dynamics in the spacing distribution of the complexLiouvillian eigenvalues, allowing a quantitative description.In the case of the ratios, the behavior of the average value oftheir modulus and phase shows excellent agreement with the .
00 0 .
25 0 .
50 0 .
75 1 . α . . . . . β Rat const Rat sep a ) γ .
00 0 .
25 0 .
50 0 .
75 1 . αb ) h r i .
00 0 .
25 0 .
50 0 .
75 1 . αc ) −h cos θ i FIG. 2. Characterizations of the integrable to chaotic transition using the γ
19 parameter (a) and the (cid:104) r (cid:105) (b) and −(cid:104) cos θ (cid:105) (c) averages of thecomplex spacing ratios 20. Each point represents a sample of 20 Liouvillians of matrix dimension . behavior of the fitted parameter value of the interpolation for-mula as we scan the transition between integrability and chaosin the parameter space of our model. These results give anextra support for the use of the complex ratios for character-izing the onset of chaos in open quantum systems. A similaranalysis in other models could help to elucidate the role ofchaos in information scrambling and thermalization in quan-tum many body open systems. We anticipate that the appli-cations of these tools will be as widespread as the equivalentspectral analysis in closed quantum many body systems. ACKNOWLEDGMENTS
We acknowledge financial support from Projects No.PGC2018-094180-B-I00 (MCIU/AEI/FEDER, EU) and andCAM/FEDER Project No. S2018/TCS-4342 (QUITEMAD-CM). This research has been also supported by CSIC Re-search Platform on Quantum Technologies PTI-001. [1] M. Srednicki, Chaos and quantum thermalization, Phys. Rev. E , 888 (1994).[2] M. Ueda, Quantum equilibration, thermalization and prether-malization in ultracold atoms, Nature Reviews Physics , 669(2020).[3] J. Maldacena, S. H. Shenker, and D. Stanford, A bound onchaos, Journal of High Energy Physics , 106 (2016).[4] R. J. Lewis-Swan, A. Safavi-Naini, J. J. Bollinger, and A. M.Rey, Unifying scrambling, thermalization and entanglementthrough measurement of fidelity out-of-time-order correlatorsin the dicke model, Nature Communications , 1581 (2019).[5] L. M. Sieberer, T. Olsacher, A. Elben, M. Heyl, P. Hauke,F. Haake, and P. Zoller, Digital quantum simulation, trotter er-rors, and quantum chaos of the kicked top, npj Quantum Infor-mation , 78 (2019).[6] F. Haake, Quantum Coherence in Mesoscopic Systems (Springer, 1991) pp. 583–595.[7] M. V. Berry, M. Tabor, and J. M. Ziman, Level clustering in theregular spectrum, Proceedings of the Royal Society of London.A. Mathematical and Physical Sciences , 375 (1977).[8] O. Bohigas, M. J. Giannoni, and C. Schmit, Characterizationof chaotic quantum spectra and universality of level fluctuationlaws, Phys. Rev. Lett. , 1 (1984).[9] M. Sieber and K. Richter, Correlations between periodic orbitsand their role in spectral statistics, Physica Scripta T90 , 128(2001). [10] S. M¨uller, S. Heusler, P. Braun, F. Haake, and A. Altland, Semi-classical foundation of universality in quantum chaos, Phys.Rev. Lett. , 014103 (2004).[11] S. M¨uller, S. Heusler, A. Altland, P. Braun, and F. Haake,Periodic-orbit theory of universal level correlations in quantumchaos, New Journal of Physics , 103025 (2009).[12] H.-J. St¨ockmann, Quantum chaos: an introduction (CambridgeUniversity Press, Cambridge, 1999).[13] M. Gutzwiller, Quantum chaos, Scholarpedia , 3146 (2007).[14] J. G´omez, K. Kar, V. Kota, R. Molina, A. Rela˜no, and J. Re-tamosa, Many-body quantum chaos: Recent developments andapplications to nuclei, Physics Reports , 103 (2011).[15] T. A. Brody, A statistical measure for the repulsion of energylevels, Lettere al Nuovo Cimento (1971-1985) , 482 (1973).[16] M. V. Berry and M. Robnik, Semiclassical level spacings whenregular and chaotic orbits coexist, Journal of Physics A: Math-ematical and General , 2413 (1984).[17] T. Prosen and M. Robnik, Semiclassical energy level statisticsin the transition region between integrability and chaos: tran-sition from brody-like to berry-robnik behaviour, Journal ofPhysics A: Mathematical and General , 8059 (1994).[18] R. Grobe, F. Haake, and H.-J. Sommers, Quantum distinctionof regular and chaotic dissipative motion, Phys. Rev. Lett. ,1899 (1988).[19] H.-P. Breuer and F. Petruccione, The Theory of Open QuantumSystems (Oxford University Press, Oxford, 2007). [20] J. Ginibre, Statistical ensembles of complex, quaternion, andreal matrices, Journal of Mathematical Physics , 440 (1965),https://doi.org/10.1063/1.1704292.[21] A. B. Jaiswal, A. Pandey, and R. Prakash, Universality classesof quantum chaotic dissipative systems, EPL (Europhysics Let-ters) , 30004 (2019).[22] R. Hamazaki, K. Kawabata, N. Kura, and M. Ueda, Universalityclasses of non-hermitian random matrices, Phys. Rev. Research , 023286 (2020).[23] G. Akemann, M. Kieburg, A. Mielke, and T. Prosen, Universalsignature from integrability to chaos in dissipative open quan-tum systems, Phys. Rev. Lett. , 254101 (2019).[24] M. V. Medvedyeva, F. H. L. Essler, and T. Prosen, Exact betheansatz spectrum of a tight-binding chain with dephasing noise,Phys. Rev. Lett. , 137202 (2016).[25] D. A. Rowlands and A. Lamacraft, Noisy spins and therichardson-gaudin model, Phys. Rev. Lett. , 090401 (2018).[26] P. Ribeiro and T. Prosen, Integrable quantum dynamics of opencollective spin models, Phys. Rev. Lett. , 010401 (2019).[27] S. Denisov, T. Laptyeva, W. Tarnowski, D. Chru´sci´nski, andK. ˙Zyczkowski, Universal spectra of random lindblad operators,Phys. Rev. Lett. , 140403 (2019).[28] L. S´a, P. Ribeiro, and T. Prosen, Spectral and steady-state prop-erties of random liouvillians, Journal of Physics A: Mathemati-cal and Theoretical , 305303 (2020).[29] M. van Caspel and V. Gritsev, Symmetry-protected coherent re-laxation of open quantum systems, Phys. Rev. A , 052106(2018).[30] F. Minganti, A. Biella, N. Bartolo, and C. Ciuti, Spectral theoryof liouvillians for dissipative phase transitions, Phys. Rev. A ,042118 (2018).[31] T. E. Lee, C.-K. Chan, and S. F. Yelin, Dissipative phase transi-tions: Independent versus collective decay and spin squeezing,Phys. Rev. A , 052109 (2014). [32] H. S. Dhar, M. Zens, D. O. Krimer, and S. Rotter, Variationalrenormalization group for dissipative spin-cavity systems: Pe-riodic pulses of nonclassical photons from mesoscopic spin en-sembles, Phys. Rev. Lett. , 133601 (2018).[33] B. Buˇca and D. Jaksch, Dissipation induced nonstationarity ina quantum gas, Phys. Rev. Lett. , 260401 (2019).[34] T. Can, V. Oganesyan, D. Orgad, and S. Gopalakrishnan, Spec-tral gaps and midgap states in random quantum master equa-tions, Phys. Rev. Lett. , 234103 (2019).[35] J. Dukelsky, C. Esebbag, and P. Schuck, Class of exactly solv-able pairing models, Phys. Rev. Lett. , 066403 (2001).[36] G. Ortiz, R. Somma, J. Dukelsky, and S. Rombouts, Exactly-solvable models derived from a generalized gaudin algebra, Nu-clear Physics B , 421 (2005).[37] J. Dukelsky, S. Pittel, and G. Sierra, Colloquium: Exactly solv-able richardson-gaudin models for many-body quantum sys-tems, Rev. Mod. Phys. , 643 (2004).[38] R. Richardson and N. Sherman, Exact eigenstates of thepairing-force hamiltonian, Nuclear Physics , 221 (1964).[39] E. Stouten, P. W. Claeys, J.-S. Caux, and V. Gritsev, Integrabil-ity and duality in spin chains, Phys. Rev. B , 075111 (2019).[40] J. M. G. G´omez, R. A. Molina, A. Rela˜no, and J. Retamosa,Misleading signatures of quantum chaos, Phys. Rev. E ,036209 (2002).[41] L. S´a, P. Ribeiro, and T. Prosen, Complex spacing ratios: A sig-nature of dissipative quantum chaos, Phys. Rev. X , 021019(2020).[42] V. Oganesyan and D. A. Huse, Localization of interactingfermions at high temperature, Phys. Rev. B , 155111 (2007).[43] T. Peron, B. M. F. de Resende, F. A. Rodrigues, L. d. F. Costa,and J. A. M´endez-Berm´udez, Spacing ratio characterization ofthe spectra of directed random networks, Phys. Rev. E102