Random Lindblad Dynamics

TL;DR

This paper investigates the mixing times of random Lindblad generators, revealing how the number of jump operators influences whether the system relaxes algebraically or exponentially, with implications for quantum dissipative dynamics.

Contribution

It introduces a novel diagrammatic expansion and exact form factor calculations for random Lindblad generators, elucidating the spectral gap's role in mixing behavior.

Findings

Single jump operator generators relax algebraically, closing the spectral gap.

Multiple jump operators or Hamiltonian terms induce a finite spectral gap.

Numerical results support the analytical predictions.

Abstract

We study the mixing behavior of random Lindblad generators with no symmetries, using the dynamical map or propagator of the dissipative evolution. In particular, we determine the long-time behavior of a dissipative form factor, which is the trace of the propagator, and use this as a diagnostic for the existence or absence of a spectral gap in the distribution of eigenvalues of the Lindblad generator. We find that simple generators with a single jump operator are slowly mixing, and relax algebraically in time, due to the closing of the spectral gap in the thermodynamic limit. Introducing additional jump operators or a Hamiltonian opens up a spectral gap which remains finite in the thermodynamic limit, leading to exponential relaxation and thus rapid mixing. We use the method of moments and introduce a novel diagrammatic expansion to determine exactly the form factor to leading order in Hilbert space dimension $N$. We also present numerical support for our main results.

Figures

Equations

$$\dot{\rho }=\mathcal{L}[\rho ],$$

$${\rho }^{\prime }=\sum _{k=1}^{N^2}{A}_k\rho A_k^{†},\sum _k{A}_k^{†}{A}_k=\text{mathbbm}1,$$

$$(1)K_{t+s}=K_t{K}_s,(2)K_0=\text{mathbbm}1,t,s\ge 0.$$

$$\mathcal{L}$$

$${\mathcal{C}}_H[\rho ]$$

$${\mathcal{C}}_H$$

$${\mathcal{D}}_L[\rho ]$$

$${\mathcal{D}}_L$$

$$K_t=e^{t\mathcal{L}},t\ge 0,$$

$$F(t)=\frac{1}{N^2}⟨\mathrm{tr}{K}_t⟩,$$

$$M_n=\frac{1}{N^2}⟨\mathrm{tr}{\mathcal{L}}^n⟩.$$

$$f(t)=\text{mbox}tr\left(\rho (t)\rho (0)\right),$$

$$⟨f(t)⟩$$

$$C(t)=\frac{N^{2}F(t)-1}{N^2-1}C(0),$$

$$\Delta =-\max \mathrm{Re}(z_i),z_i\in \sigma (\mathcal{L})/\{0\},$$

$$\Delta =-\underset{t\to \infty }{\lim }\frac{1}{t}\log \left[\frac{1}{N^2}\left(\text{mbox}tr{K}_t-1\right)\right].$$

$$⟨\Delta ⟩=-\underset{t\to \infty }{\lim }\frac{1}{t}\log F(t)+O(N^{-1}).$$

$${\tilde{\mathcal{D}}}_L=-L^{†}L{\otimes }_t\text{mathbbm}1-\text{mathbbm}1{\otimes }_t{L}^{†}L,$$

$$M_n=\underset{N\to \infty }{\lim }\frac{1}{N^2}⟨\mathrm{tr}{\tilde{\mathcal{L}}}^n⟩+O(N^{-1}).$$

$$F(t)=\frac{1}{N^2}⟨\mathrm{tr}{e}^{{\mathcal{D}}_{L}t}⟩=e^{-4vt}{\left(I_0(2vt)+I_1(2vt)\right)}^2+O(N^{-1}),$$

$$F(t)=\frac{1}{N}+\frac{1}{t^2}{\left(1-{\left(1+\frac{t}{N}\right)}^{-N}\right)}^2-\frac{1}{t^2}\left(1-\frac{(1+2t/N{)}^N}{(1+t/N{)}^{2N}}\right).$$

$$F(t)=C{t}^{-3}{e}^{-\Delta vt}+O(N^{-1}),$$

$$\Delta =2(1-\sqrt{m}{)}^2.$$

$$F(t)\le \frac{1}{4{\pi }^2}{t}^{-4}{e}^{-\Delta t}+O(N^{-1}),t\to \infty ,$$

$$\mathcal{D}=2L{\otimes }_t{L}^{†}-L^{†}L{\otimes }_t\text{mathbbm}1-\text{mathbbm}1{\otimes }_t{L}^{†}L,$$

$$⟨{\bar{L}}_{ij}{L}_{kl}⟩=\frac{v}{N}{\delta }_{ik}{\delta }_{jl},L_{ij}\in \text{mathbbm}C.$$

$$\sum _{i=1}^N{\mathcal{D}}_{iikl}=0.$$

$${\mathcal{D}}_{ijkl}=2L_{ik}{L}_{lj}^{†}-(L^{†}L{)}_{ik}{\delta }_{jl}-(L^{†}L{)}_{lj}{\delta }_{ik}.$$

$$\sum _i{\left({\mathcal{D}}^n\right)}_{iikl}=0,$$

$${\left({\mathcal{D}}^2\right)}_{ijkl}={\mathcal{D}}_{ijab}{\mathcal{D}}_{abkl},$$

Full text

1. 1 Introduction

2. 1.1 Basic Definitions and Main Results

3. 1.2 Discussion and Background

4. 2 Simple Dissipator with a Random Jump Operator

5. 2.1 Complex Ginibre Ensemble

6. 2.2 Complex Random Normal Matrices

7. 2.3 Gaussian Unitary Ensemble

8. 3 Dissipators with Multiple Jump Operators

9. 4 Simple Lindbladian

10. 5 Conclusion

11. A Asymptotic Decay Rate and Spectral Gap

12. B Complete Diagrams for Second Moment of the Simple Dissipator

13. C Topological Expansion and Proof of Lemma (1)

14. D Random Normal Matrices: Correlation Functions and Exact Form Factor

15. E Random Liouville-von Neumann Generator

16. F Random Non-Hermitian Hamiltonian: Resolvent and Eigenvalue Density

17. G Non-Hermitian Hamiltonian with Multiple Jump Operators

1 Introduction

A foundational problem in statistical mechanics concerns the existence, or lack thereof, of thermalization. Indeed, Boltzmann was among the first to address the question of whether many-body systems reach equilibrium Boltzmann1872 . Closely related to this issue is the question of how non-equilibrium systems approach equilibrium, in particular the dynamics of relaxation.

Since Boltzmann, the conditions under which thermalization is achieved have been clarified considerably, and there now seems to be a clear picture of the roles of chaos, ergodicity, and integrability in determining whether thermalization occurs in classical and quantum systems. However, there seems to be much left to say about the story of non-equilibrium relaxation, and characterizing for instance the typical time scales involved in the approach to equilibrium Goldstein2013 , and the effects of metastability Lesanovsky2013 ; Macieszczak2016 .

Thermalization in classical systems is intimately tied to non-integrability and chaoticity of the microscopic dynamics. The quantum manifestation of chaos is the appearance of a universal random matrix theory (RMT) description of late-time dynamics and spectral statistics Berry1985 ; Bohigas1984 . The eigenstate thermalization hypothesis (ETH) provides the bridge between the universal RMT behavior and thermalization in quantum statistical mechanics Srednicki1999 ; Rigol2008 ; Dymarsky2018a ; DAlessio2016 ; Schiulaz2018 .

While it is valuable to clarify this story for isolated systems, they are still rather the exception in nature. It is natural to wonder if driven and dissipative systems also exhibit universal dynamical or steady state properties. Scattering theory has proven to be a fruitful approach to addressing the question of universality in open quantum systems. A conjecture for open chaotic scattering involving the distribution of phase shifts was proposed in Blumel1990 , as the open analog of the Bohigas-Giannoni-Schmit conjecture for level statistics in chaotic billiards Bohigas1984 . In the scattering problem, the object of interest is the scattering $S$ matrix whose poles are the scattering resonances. For a chaotic scattering region, the $S$ matrix is described by Dyson’s circular ensembles of random unitary matrices Dyson1962a ; Schomerus2016 . Furthermore, the resonance width distribution in the limit of weak coupling to the continuum is described by the Porter-Thomas distribution, which is also reproduced by random matrix theory Fyodorov1997a (for reviews, see e.g. Schomerus2016 ; Fyodorov2010 ). Many of these results have also received experimental support Alt1995 ; Kuhl2005 ; Kuhl2008 ; Kuhl2005a .

The dynamical signatures of quantum chaos in the presence of dissipation and decoherence have been taken up in, e.g. Zurek1994 ; Miller1998a ; Casati1997 ; Waltner2018 . In the context of open chaotic scattering, the time evolution of decay functions or survival probabilities have been studied using random matrix theory Dittes1992 ; Harney1992 ; Dittes2000 .

Motivated by the success of RMT in describing the dynamical signatures of thermalization in closed quantum systems, we seek to expand this to open quantum systems by studying a random matrix theory of quantum master equations. In particular, we are interested in the approach to steady state, and the typical relaxing behavior of open quantum systems. Our results appear to be consonant with the literature on classical open chaotic scattering, suggesting that these domains can be bridged, and that a RMT description of open quantum master equations can capture universal relaxation dynamics. This discussion mainly serves to provide a larger context for our work, and we leave the precise connection to quantum chaos as a direction for future research.

The spirit of the present paper is more in line with the original approach in nuclear physics to develop a statistical description of complex quantum phenomena. A statistical description of deterministic processes can only ever have a chance of being meaningful when there are many degrees of freedom involved, allowing us to retain only the essential structure (e.g. hermitian vs non-hermitian) and symmetries (e.g. time-reversal, parity, etc. Dyson1962 ; Altland1997 ), while eschewing complex microscopic details Dyson1962a ; Wigner1955 ; Dyson1962 . These are then maximum entropy ensembles subject to certain constraints imposed by the physical problem of interest Balian1968 . This rule of thumb is borne out by the success of RMT in describing chaotic systems, whose motion in phase space is ergodic and highly mixing, as well as heavy nuclei which are highly complex interacting systems with many degrees of freedom. The universality of the statistical description is observed in properties which depend only on the general structure and symmetries. In developing a statistical description of non-equilibrium dissipative quantum dynamics, the ensembles we consider in this paper are the maximum entropy ensembles of generators of the quantum dynamical semigroup, which govern the quantum Markov master equation, more familiarly known as the Lindblad equation. The general representation of Lindblad generators splits the evolution into a unitary part, governed by a Hamiltonian, and a dissipative part involving jump operators. Selecting the Hamiltonian and jump operators from a random matrix ensemble then naturally leads to a random Lindblad generator. Furthermore, one may consider discrete symmetries as in Buca2012 ; albert2014symmetries , and their effect on the statistical properties of the master equation; we save this for future work, and consider here only ensembles without any additional symmetries.

We now proceed to review the relevant results from the theory of open quantum systems which we will use to define our random matrix ensembles and introduce our main object of study. We refer to Petruccione for more details on this standard textbook material.

Open quantum systems

Quantum master equations provide a framework to study driven and dissipative quantum systems through the continuous time evolution of the density matrix

[TABLE]

where the Liouvillian $\mathcal{L}$ is in general an integral operator, non-local in time, and possibly non-linear in $\rho$. Formally integrating the master equation one can define the dynamical map $K_{t}:\mathcal{S}(\mathcal{H})\to\mathcal{S}(\mathcal{H})$ as a map acting on the space of density matrices $\mathcal{S}(\mathcal{H})$ over the Hilbert space $\mathcal{H}$. We consider finite but large $N$-dimensional Hilbert spaces. The Liouvillian is a superoperator acting as a linear map on the space of bounded operators $\mathcal{B}(\mathcal{H})=\mathcal{H}\otimes\mathcal{H}^{*}$. It also has a non-unique $N^{2}\times N^{2}$ matrix representation $\mathcal{L}_{ab}={\rm tr}(e_{a}^{\dagger}\mathcal{L}[e_{b}])$ whose specific form depends on the basis vectors $e_{a}\in\mathcal{B}(\mathcal{H})$. We use the standard basis below.

There exists a useful operator sum representation of linear maps on the space of density matrices due to Kraus Kraus1971 (and independently Choi Choi1975 ). Such maps have come to be known as quantum operations or channels. The space of density matrices consists of complex-valued positive-definite hermitian matrices with unit trace. Any physically sensible dynamical evolution of the density matrix must preserve these properties. Kraus’ theorem states that any completely positive and trace preserving (CPTP) quantum channel which takes $\rho\to\rho^{\prime}$ can be represented by

[TABLE]

for some Kraus operators $A_{k}$ subject to the stated constraint. The crucial assumption which allows the explicit representation of the Kraus map is complete positivity, which ensures that every trivial extension of the map to an enlarged Hilbert space, which acts as the identity on this additional space, will preserve the positivity of the density matrix.

A subset of quantum channels are those which are parametrized by a continuous time variable and form a semigroup. A dynamical map forms a quantum dynamical semigroup (QDS) if it satisfies the following axioms

[TABLE]

Axiom (1) is equivalent to the Markov property, and implies that a QDS necessarily describes Markovian dynamics. The Markovian assumption allows us to utilize a fundamental theorem independently discovered by Lindblad Lindblad1976 and Gorini, Kossakowski, & Sudarshan (GKS) Gorini1976 , concerning the classification of the generators of a QDS. The generator $\mathcal{L}$ of the QDS can be written (in the “diagonal" Lindblad basis)

[TABLE]

We refer to the Liouvillian in this form as the Lindblad generator or superoperator, or simply the Lindbladian 111Also known as the GKS-Lindblad or GKLS generator. We use generator and superoperator interchangeably when referring to the Lindbladian.. We have split it up into two pieces: the first describes unitary Liouville-von Neumann evolution under a (Hermitian) Hamiltonian $H$

[TABLE]

which we refer to as the Liouville-von Neumann (LvN) generator for $H$. We work mainly with the following matrix representation

[TABLE]

where $\mathbbm{1}$ is the identity matrix of dimension $N$. We have introduced the notation $\otimes_{t}$ for clarity of presentation, which is defined $A\otimes_{t}B\equiv A\otimes B^{T}$, with $\otimes$ the conventional tensor (Kronecker) product of matrices, and $B^{T}$ denoting the matrix transpose of $B$. In this paper, we work almost exclusively with the matrix representation (6), and distinguish the superoperator as a map by explicitly including an argument, as in (5). The dissipative part of the Lindbladian is described by a sum of dissipators $\mathcal{D}_{L}$ for the operator $L$ defined by the superoperator

[TABLE]

which has the matrix representation

[TABLE]

The theorem imposes no constraints on the “jump operators" $L_{a}$. The dissipative couplings $\gamma_{a}$ (also called dephasing rates), are required to be positive to ensure the convergence of the flow to a stationary state. The most general CPTP map of the Kraus form has $N^{4}-N^{2}$ independent parameters, which translates to a general dissipator consisting of $m=N^{2}-1$ jump operators, or dissipation channels. However, the number of the jump operators, and the shape they take, is often informed by the physical problem at hand. We only consider Lindbladians with a number of dissipation channels that remains finite in the thermodynamic limit.

Sometimes, the entire dissipative contribution to the Lindblad generator is collectively referred to as the dissipator. To avoid confusion, we refer to (8) as the simple dissipator (also called the simple generator Baumgartner2008 ). Following the literature, we refer to the first term in (8) as the “recycling term".

Given the GKS-Lindblad representation of the generator for the quantum Markov master equation, we can define an ensemble of generators by taking the Hamiltonian and jump operators to be random matrices. In this paper, we consider matrices without any additional discrete symmetries, and for simplicity assume Gaussian distributed matrix elements. Thus, we take the Hamiltonian from the Gaussian unitary ensemble (GUE) of complex Hermitian matrices. We take the jump operators from the complex Ginibre ensemble of complex matrices with no symmetries. We also present results for the simple dissipator with jump operators drawn from the GUE and normal random matrix ensembles.

In this work, we study relaxation times under random Lindblad evolution by focusing on the dynamical map or propagator $K_{t}$. Our approach is sensitive to the important question of whether relaxation occurs exponentially fast, and is thus rapidly mixing, or algebraically in time. From this relaxation behavior we can infer the mean spectral gap in the large $N$ limit. We develop a diagrammatic approach to studying random Lindblad generators using the method of moments. Our main results are summarized in the next section.

1.1 Basic Definitions and Main Results

The primary object of study in this paper is the dynamical map or propagator

[TABLE]

which forms a single parameter representation of a quantum dynamical semigroup. We take the Liouvillian to be of Lindblad form. Since we consider a finite $N$-dimensional Hilbert space, we work with the $N^{2}\times N^{2}$ matrix representation of the Lindblad superoperator, in which the propagator is just the exponential of the Lindbladian. We access asymptotic relaxation by studying the dissipative form factor (DFF) for the Lindbladian:

Definition 1** (Dissipative Form Factor).**

The dissipative form factor is defined as the ensemble averaged trace of the propagator (9)

[TABLE]

normalized such that $F(0)=1$. Angular brackets indicate averaging over the relevant ensemble of Hamiltonian and/or jump operator(s).

The DFF reduces to the spectral form factor (normalized by $N^{2}$) for unitary Liouville-von Neumann evolution without a dissipator. The form factor is nothing but the moment generating function for the Lindblad generator, where the trace moments are defined

[TABLE]

Here, $\mathcal{L}^{n}$ is the $n^{th}$ matrix power of the Lindblad generator, which is equivalent to the matrix representation of the $n$-fold composition of the superoperator $\mathcal{L}:\mathcal{B}(\mathcal{H})\to\mathcal{B}(\mathcal{H})$.

The form factor (10) is closely related to Uhlmann’s “transition probability"Uhlmann1976 , or mixed state fidelity Jozsa1994

[TABLE]

which provides a measure of decoherence due to interaction with an environment 222Recently dubbed the Uhlmann fidelity Tonielli2018 , but also referred to as the mixed-state memory fidelity Lidar1998; Bacon1999, relative purity DelCampo2013 , linear fidelity, or simply the overlap Audenaert2012.. In the thermodynamic limit, where the system’s Hilbert space dimension $N\to\infty$, the ensemble-averaged fidelity for an initial pure state is given by

[TABLE]

The asymptotic behavior of Uhlmann’s fidelity was considered in Tonielli2018 and shown to exhibit non-trivial scaling behavior tied to the non-orthogonality of eigenvectors of the Lindblad generator. As a result of the ensemble averaging, the dissipative form factor, and thus the fidelity, does not carry any direct information about eigenvectors. Our results for $F(t)$ can be imported into (13) to recover the long-time behavior of the fidelity. In Refs. Miller1998a ; Tameshtit1993 , (13) was referred to as the survival probability function DelCampo2013 . The adjoint Lindblad generator $\mathcal{L}^{\dagger}$ (which is the conjugate transpose in the matrix representation) governs the time evolution of observables $A(t)=e^{t\mathcal{L}^{\dagger}}A(0)$, and can be used to compute ensemble-averaged autocorrelation functions $C(t)=\langle A(t)A(0)\rangle$ (assuming $A(0)$ is traceless for simplicity)

[TABLE]

once again illustrating the central role of the DFF in the ensemble-averaged dynamics. Most significantly, such autocorrelation functions can be measured in NMR experiments on echo dynamics Goussev2012 .

It is also useful to introduce here the spectral gap of the Lindblad superoperator:

Definition 2** (Spectral Gap).**

For finite $N$, the spectral gap is defined as

[TABLE]

where $\sigma(\mathcal{L})/\{0\}$ denotes the eigenvalue spectrum of the Lindbladian excluding the zero modes.

Since ${\rm Re}(z_{i})\leq 0$ for all eigenvalues of $\mathcal{L}$, the spectral gap is a positive real number. At finite $N$ and for a quenched random Lindbladian, the spectral gap follows directly from the limit 333This definition is sufficient for our purposes since the steady state is unique and purely imaginary eigenvalues require additional structure and symmetry, and thus have a vanishingly small probability to occur in our ensembles. ,

[TABLE]

In Appendix (A), we argue that in the large $N$ limit, the ensemble averaged spectral gap will be given by

[TABLE]

This can be viewed as a reflection of the self-averaging property of the trace of the propagator. Consequently, we identify the asymptotic decay rate, defined by the right-hand-side of Eq.(17) with the expected value of the spectral gap. Exponential relaxation of the DFF then implies a nonvanishing average spectral gap, whereas algebraic relaxation implies a vanishing gap.

We state our main results as propositions below. To get there, we must first introduce the concept of a non-crossing truncation of the Lindblad superoperator.

Definition 3** (non-crossing truncation).**

The non-crossing truncation of the simple dissipator is defined by

[TABLE]

which consists of removing the recycling term $2L\otimes_{t}L^{\dagger}$. The non-crossing truncation of the Lindbladian $\tilde{\mathcal{L}}$ consists of replacing all simple dissipators with their non-crossing truncation.

The terminology “non-crossing" comes from the diagrammatic analysis we develop in Sec.(2.1), where its meaning is made clear. Briefly, the density matrix can be viewed as an element of $\mathcal{H}_{L}\otimes\mathcal{H}_{R}$, where $\mathcal{H}_{L}$ and $\mathcal{H}_{R}$ are left (bra) and right (ket) identical copies of the Hilbert space. Some operators appearing in the Lindbladian act as the identity on one of these spaces, such as the term $L^{\dagger}L\otimes_{t}\mathbbm{1}$ which acts trivially on the right Hilbert space. On the other hand, the recycling term $L\otimes_{t}L^{\dagger}$ acts nontrivially on both the left and right Hilbert space. The non-crossing truncation involves just removing the terms in the Lindbladian which act nontrivially on both left and right Hilbert spaces; these are precisely the recycling terms.

Our next result is crucial for establishing most of the following propositions (unless otherwise noted):

Lemma 1**.**

Let the jump operators be complex Ginibre random matrices Ginibre1965 . To leading order in $N$, the trace moments of the Lindblad superoperator are determined by the non-crossing truncation

[TABLE]

Our proof of this lemma relies on the graphical calculus developed in Sec.(2.1). The proof is completed in two stages: first for a single jump operator in Sec.(2.1), and next for multiple jump operators in Sec.(3). These proofs are recapitulated using a topological expansion in Appendix (C). The proof of the lemma for a simple dissipator is a stepping stone to the first proposition:

Proposition 1** (Simple Dissipator: complex Ginibre ensemble).**

Consider a simple generator (8) with the jump operator $L$ a complex Ginibre random matrix, i.e. $L_{ij}\in\mathbbm{C}$ are i.i.d. complex Gaussian variables with variance $\langle|L_{ij}|^{2}\rangle=v/N$. The dissipative form factor for this generator is

[TABLE]

where $I_{n}(x)$ is the modified Bessel function of the $n^{th}$ kind.

In particular, at long times, the dissipative form factor for the simple dissipator decays as a power-law $F(t)\sim t^{-1}$. Therefore, the spectral gap closes for a simple dissipator with complex Ginibre jump operators. We present numerical evidence that this closing occurs as $\Delta\sim N^{-2}$ in Sec.(2.1), implying a diffusive dynamical exponent for the relaxation time $\tau\sim N^{z}$ with $z=2$ Gier2006 . It is also worth mentioning that $F(\infty)=N^{-2}$, as a consequence of the existence of a unique stationary state $\mathcal{D}_{L}[\rho_{ss}]=0$ Baumgartner2008 .

A close relative of the complex Ginibre ensemble is the ensemble of random normal matrices Zabrodin2004 . Both ensembles have the same eigenvalue joint probability distribution function, but their effect on the phenomenology of the dissipator is quite distinct. While the dissipator with Ginibre jump operators has a unique steady state, there are exactly $N$ stationary states for normal matrices. The dissipative form factor also exhibits different power-law relaxation in the two cases, which we state presently.

Proposition 2** (Simple Dissipator: Random normal matrix).**

Consider a simple generator (8) with the jump operator $L$ a random normal matrix $[L,L^{\dagger}]=0$. The exact dissipative form factor for finite $N$ is

[TABLE]

At $N=\infty$, the second term takes over and $F(t)=(1-e^{-t})^{2}/t^{2}$.

The derivation of (2) appears in Appendix (D), and does not require Lemma (1), using instead exact results for the eigenvalue density correlation functions in the random normal matrix ensemble Zabrodin2004 .

The only difference between the Ginibre and normal matrix ensembles is the eigenvector statistics, which is trivial for the latter. It is natural that non-orthogonality of eigenvectors in the Ginibre case leads to a slower rate relaxation. Furthermore, the relaxation time with normal jump operators scales as $\tau\sim\sqrt{N}$, which means the steady state is achieved much sooner than in the Ginibre case.

For the sake of comparison, we also discuss the simple dissipator with GUE jump operators in Sec. (2.3). The spectral gap here also closes, and the form factor scales as $F(t)\sim t^{-1/2}$ at long times. Nevertheless, the relaxation time $\tau\sim N^{2}$, as in the Ginibre case.

The simple dissipator turns out to be rather exceptional in its algebraic relaxation to equilibrium. As we see next, considering even a single additional jump operator in (8) immediately leads to exponential relaxation, indicative of the opening of a thermodynamic spectral gap.

Proposition 3** (Multiple Jump Operators).**

Let $H=0$ and $m$ be finite such that $m/N\to 0$ as $N\to\infty$. For $a=1,...,m$, let the jump operators $L_{a}$ be independent complex Ginibre matrices with identical variance as in Prop.(1). The dissipative form factor in the long-time limit is given by

[TABLE]

for a constant $C=O(1)$ given in Eq.(90). Here, the asymptotic decay rate $\Delta$ is determined by

[TABLE]

The proof requires Lemma (1) and appears in Sec.(3) (see also Appendix (C)). The exponential relaxation leads to a nonzero average spectral gap, per Eq.(17).

The existence of a unique steady state for the dissipator with multiple jump operators implies again the asymptotic limit $F(\infty)=N^{-2}$ Baumgartner2008a . For finite systems, the exponential relaxation indicates that this expression is valid on relaxation time scales which are roughly $\tau\sim\frac{2}{\Delta}\log N$.

Finally, we consider a simple Lindbladian which consists of a Hamiltonian and a single jump operator.

Proposition 4** (Simple Lindbladian).**

Let the Hamiltonian $H$ be drawn from the Gaussian unitary ensemble, i.e. $H_{ij}\in\mathbbm{C}$ are i.i.d. complex Gaussian variables such that $H_{ij}=\bar{H}_{ji}$, with zero mean and variance $\langle|H_{ij}|^{2}\rangle=1/N$. Let $m=1$ with $L$ a complex Ginibre matrix with variance $\langle|L_{ij}|^{2}\rangle=1/N$. Consider the simple Lindbladian $\mathcal{L}=\mathcal{C}_{H}+\gamma\mathcal{D}_{L}$, where the first and second terms are defined in Eq. (5) and Eq. (8), respectively, with $\gamma\in\mathbbm{R}^{+}$ a positive real-valued constant independent of $N$.

The large $N$ form factor at long times is bounded above by

[TABLE]

with the asymptotic decay rate $\Delta=2\delta$, where $\delta$ is the spectral gap of the corresponding non-Hermitian Hamiltonian $H-i\gamma L^{\dagger}L$.

The inequality for the DFF implies a lower bound on the average spectral gap which is finite in the thermodynamic limit. That the random non-Hermitian Hamiltonian should have a spectral gap is a remarkable fact first discovered in Haake1992 . This proposition shows that the leading order behavior of the form factor is determined by the non-unitary dynamics induced by the effective non-Hermitian Hamiltonian.

Lemma (1) is necessary to show that the asymptotic decay rate is determined by the non-Hermitian Hamiltonian, and the derivation of the inequality for the resulting form factor is given in Sec.(4). The functional form of this decay rate was found in Haake1992 ; Lehmann1995 , and is re-derived in Appendix (F) using straightforward algebraic techniques. In particular, the decay rate has the asymptotic scaling $\Delta\sim 2\gamma$ for small $\gamma$, and $\Delta\sim 2^{1/3}(\gamma)^{-1/3}$ at large dissipative coupling. We conjecture that the spectral gap for the simple Lindbladian will obey the same scaling behavior. The observation of exponential decay at large dissipative coupling is reminiscent of the quantum Zeno effect (see e.g. Mensky2013 ; Popkov2018 ), though the non-analytic scaling of $\Delta(\gamma)$ is perhaps unusual in this context.

The form factor for the simple Lindbladian is in fact a damped oscillating function of time. In the $\gamma=0$ limit, the unitary evolution produces a dynamical scaling $t^{-3}$, well known from the spectral form factor averaged over the GUE (see (159)). In the opposite limit, Prop.(1) shows that the power-law becomes $t^{-1}$.

In Appendix (G) we comment briefly on how the asymptotic decay rate, and thus the mean spectral gap, is expected to change when there are $m$ jump operators and a Hamiltonian. In this case, we find $\Delta\sim 2m\gamma$ for small $\gamma$, and $\Delta\sim 2\gamma(1-\sqrt{m})^{2}$ for large $\gamma$.

Our methods do not grant access to the detailed structure of the eigenvalue distribution of the Lindblad superoperator. The reason is that the moment generating function (DFF) is the Laplace transform of the holomorphic Green’s function, whose domain of validity is outside the support of the eigenvalue density. Nevertheless, the holomorphic Green’s function is sufficient to extract one of the most basic properties of Lindblad evolution, which is its mixing behavior, and can thus identify the existence of a spectral gap.

1.2 Discussion and Background

Our main results support the claim that generic Lindblad evolution relaxes exponentially, and is thus rapidly mixing. This behavior gives strong evidence for the existence of a spectral gap in the spectrum of the Lindblad superoperator. Exponential relaxation of the form factor in the limit of large $N$ implies that the bulk spectrum (sometimes called the essential spectrum) is separated from the imaginary axis by a finite gap. However, it does not exclude the possibility of a subextensive $o(N^{2})$ number of isolated eigenvalues residing within this bulk spectral gap (see Dellnitz2000 for examples of such isolated eigenvalues of the Frobenius-Perron operator in classical dynamics). Such isolated eigenvalues were indeed found to determine the asymptotic decay rate in a companion paper Ref.can , which considers a simple Lindbladian with $H$ and $L$ drawn from the Gaussian orthogonal ensemble (GOE). In particular, we found that while the spectral gap is always finite in the thermodynamic limit, for large values of dissipative coupling the spectral gap is determined by a single isolated eigenvalue which appears to split off from the bulk spectrum. While preliminary investigations do not reveal such an isolated mode for the ensembles considered in this paper, its occurrence is not strictly ruled out.

On the other hand, algebraic relaxation of the dissipative form factor is a clear signature of the spectral gap closing in the thermodynamic limit. We have seen that this occurs for simple dissipators (Props.(1) and (2)), whether the jump operators are non-hermitian, normal, or hermitian, and have checked this for these three cases with full random matrices. The power-law is different for all three, scaling as $t^{-1}$ (Ginibre), $t^{-2}$ (normal) and $t^{-1/2}$ (GUE). These disparate relaxing behaviors lead to relaxation times of the form factor which scale as $\tau\sim N^{z}$, with a diffusive dynamical exponent $z=2$ for the Ginibre and GUE, and $z=1/2$ for the normal jump operators.

For the Hermitian jump operators, the relaxation time is in fact the Heisenberg time, since the mean spacing of eigenvalues of the dissipator is $N^{-2}$. This interpretation is not directly applicable to the dissipator with non-hermitian jump operators, since in this case the eigenvalues are complex valued. However, as we show later, the complex Ginibre relaxation is controlled by its “non-crossing truncation", which consists of Wishart matrices. In this setting, the level spacing once again scales as $N^{-2}$. Such relaxation times are observed in closed quantum systems as well Dymarsky2018a ; Schiulaz2018 .

As for the random normal matrix, the unusual scaling of the relaxation time suggests that the asymptotic relaxation rate is controlled by the boundary of the distribution of eigenvalues, along which the eigenvalue density (which is a proxy for level spacing) is expected to scale like $1/\sqrt{N}$.

Finally, we remark that the time scales we have access to in our asymptotic analysis are likely semi-classical in nature. Here we give some heuristic justification for this statement. In Ref.Casati1997 , it was argued that for open chaotic systems, the survival probability, which is known classically to exhibit exponential decay, experiences significant quantum fluctuations on an intermediate time scale that can be parametrically smaller than the Heisenberg time, but still diverging with $N$. The classical dynamics then holds for times up to this quantum time scale, which is set by the mean spacing between scattering resonances on the complex plane (as opposed to the energy level spacing which sets the Heisenberg time). In our asymptotic analysis, we work with the limit in which the mean spacing of the Lindblad eigenvalues tends to zero. Furthermore, the effective Hamiltonian describing the simple Lindbladian (i.e. $K$ in Eq.(92)) can be interpreted in the context of scattering theory, and describes $N$ bound states coupled via $N$ decay channels to the continuum. An extensive number of decay channels corresponds to a semi-classical limit in the chaotic scattering problem, and generically leads to exponential decay in e.g. bound state survival probability Lewenkopf1991 ; Savin1997 . Finally, our results are qualitatively similar to those appearing in classical chaotic scattering Motter2002 . For these reasons, we speculate that large $N$ random Lindbladians display dynamics which is morally classical, despite not necessarily having a clear classical limit. This conclusion is also supported by the argument in Ref.xu2019extreme that random dissipators (specifically with Hermitian jump operators) have a maximal decoherence rate (at short times), and thus should naturally lead to classical evolution at late times.

Rapid Mixing and Spectral Gap in Open Quantum and Classical Systems

Having summarized and discussed the main results of our paper, we now proceed to attempt giving our work meaning by properly fleshing out the context. There is a deep and extensive literature on related questions in the setting of classical and quantum dynamics that we now review.

In the broadest context, our results concern the property of mixing in dynamical systems, which appears commonly in ergodic theory openbook2006 . The mixing time is the characteristic time scale required to reach the stationary state or equilibrium. For classical and quantum Markov processes, exponential relaxation toward steady state is referred to as rapid mixing Kastoryano2013 ; Lucia2015 .

For a quantum Markov process, the rapid mixing property is controlled by the spectral gap (15) of the Lindblad generator, sometimes referred to as the dissipative gap Znidaric2015 . For finite-dimensional systems, the spectral gap determines the slowest decay rate. In the thermodynamic limit, the spectral gap might close, leading to the possibility of observing algebraic relaxation of certain observables or correlations functions Medvedyeva2014 ; Cai2013 , and also potentially signaling a dissipative quantum phase transition Prosen2008 ; Kessler2012 . The role of the spectral gap in relaxation times of Lindblad dynamics was studied extensively for a particular class of boundary-driven Lindblad evolution in Znidaric2015 .

The gap also has implications for the stability of the resulting steady state to perturbations of the generator of time evolution Cubitt2015 ; Lucia2015 , as well as for the clustering of correlations between local observables Poulin2010 ; Kastoryano2014 . In Ref.Temme2013 , the existence of a spectral gap is proven for a certain class of Davies generators, which are generators of Lindblad form constructed explicitly to relax to thermal equilibrium.

Random quantum channels are also known to exhibit a spectral gap. In this setting, the steady state is an invariant state under the map, and has eigenvalue equal to unity. The rest of the eigenvalues are contained in the unit disk. If the spectral radius (i.e. largest eigenvalue in absolute value) is strictly less than unity, these operators are said to have a spectral gap. Using the Kraus representation, Ref.Bruzda2008 ; Bruzda2010 proposed a random ensemble of quantum channels based on their action on the orthonormal basis of $SU(N)$ generators. The result for a quantum channel with $m$ Kraus operators was a spectrum resembling the real Ginibre matrix ensemble Ginibre1965 , supported on a circle with radius $1/\sqrt{m}<1$.

Random quantum channels induced by Haar-distributed random isometries were introduced in Hayden2008 and studied extensively in the quantum information context Collins2015 . Recently in Ref.Gonzalez-Guillen2018 , such maps were shown to have a spectral gap in the limit of large Hilbert space dimension and for a fixed number of Kraus operators. This follows on previous work Hastings2007 which proves a spectral gap for a quantum channel with a fixed number of Haar distributed unitaries.

A similar approach considers maps describing open quantum systems obtained by projecting an $N$ dimensional random unitary onto a smaller $M$ dimensional subspace, resulting in a “truncated" unitary matrix. The spectrum of truncated random unitaries was studied in Zyczkowski2008 , and also shown to exhibit a spectral gap in the $N\to\infty$ limit when $M/N$ is kept fixed.

The spectral gap appears to be a manifestation of quantum chaos, whereas its closing indicates the emergence of integrability. In Ref.Novaes2013 , the projection technique was used to quantize a classical open chaotic map. They observed a spectral gap in the quantum map, indicating a minimal resonance lifetime. Ref.Nonnenmacher2011 determines a condition for the existence of a spectral gap in quantized open chaotic maps (i.e. quantum maps which have a clear chaotic classical limit) related to a quantity known as the topological pressure, characterizing the dynamics that is recovered in the classical limit Gaspard1989 .

In the classical setting, the dynamical map acting on the probability density function is the Frobenius-Perron (FP) operator, and the eigenvalues are known as the Pollicott-Ruelle (PR) resonances. Refs. Nonnenmacher2003 and Saraceno2004 consider quantum maps with well-defined classical limits, and found that the PR resonances governed the asymptotic decay rates of the quantum system. In Nonnenmacher2003 , it was argued that a spectral gap in the quantum map is symptomatic of ergodicity and mixing of its classical limit, whereas a closing spectral gap appeared to coincide with integrability. Classical chaotic scattering also exhibits a gap in the distribution of resonances, related to the topological pressure, also observed experimentally Barkhofen2013 .

Indeed, the dynamical consequences of proximity to integrability have been studied in the classical setting. Ref.Vivaldi1983 observe a crossover from exponential to power-law decay of auto-correlation functions in the stadium billiard, connected to the presence of arbitrarily long regular motion in the evolution of stochastic (ergodic) orbits.

Classical chaotic scattering is also known to show a transition from exponential to power law decay of survival probability of a particle in a scattering region, depending on the existence of KAM tori in phase space (see e.g. Barra2002 or the monograph Gaspard2005 ). Our results are curiously consonant with the conclusions of Ref. Motter2002 on classical chaotic scattering, which argues that algebraic decay is unstable to the presence of dissipation, and exponential decay takes over for any amount of dissipation. In the context of open quantum systems, we find that Hamiltonian evolution with the addition of dissipation leads to exponential decay of fidelity (13), which becomes the survival probability for pure Hamiltonian evolution. These results suggest that random Lindblad dynamics in the large $N$ limit are somehow semi-classical in nature, and capture qualitative features of fully chaotic classical systems.

There do exist some noteworthy connections between dissipative chaotic systems and random matrix theory. In dissipative quantum systems with a chaotic classical limit, the eigenvalues of the dynamical map were argued to exhibit a universal cubic level repulsion Grobe1989 , reviewed in the book Haake2001 (see also Braun2001 ). A semiclassical expression involving periodic orbits was found for the trace of a discrete time dynamical map for the dissipative kicked rotor (analogous to our dissipative form factor) in Braun1999 , (for a nice review, see the monograph Braun2001 ). The spectral gap in the dissipative kicked rotor model was shown to be determined by the leading Pollicot-Ruelle resonance of the FP operator Braun2001 ; Braun1999b .

Quantum dynamics in the presence of random environments or interactions has also been studied in the past ko1976one ; lebowitz2004random . Notably, in Gorin2004 the fidelity is studied in the linear response regime assuming a random matrix perturbation of the time-reversed Hamiltonian. Their result shows the universal relaxation behavior of quantum chaotic systems. It would be interesting to connect these results to the mixed state fidelity (13) for random Lindblad evolution.

Finally, we mention a conjecture for the Lindbladian spectral gap motivated by the comparison to chaotic scattering theory. It is known that for a small number of decay channels (non-extensive in Hilbert space dimension), the survival probability asymptotically decays algebraically in time Lewenkopf1991 . Translating to the Lindblad case, this corresponds to the setting in which the jump operators are low-rank matrices. At present, we are not able to extend our diagrammatic expansion to low-rank random matrices, but we can make the following conjecture: the spectral gap of the simple Lindbladian will close in the large $N$ limit for jump operators with a fixed rank that does not grow with $N$.

2 Simple Dissipator with a Random Jump Operator

2.1 Complex Ginibre Ensemble

In this section we provide part of the proof of Lemma (1), and as a first application prove Prop. (1). The techniques which are developed in this section will be called upon later to provide the proofs of Props. (3) and (4). Our approach to studying the dissipative form factor is to develop a diagrammatic calculus for calculating moments of the dissipator (for now dropping the subscript appearing in (8) for notational convenience)

[TABLE]

when the jump operator is a complex Ginibre random matrix with i.i.d. complex Gaussian entries with zero mean and variance

[TABLE]

We begin with some descriptive numerical results to help us anticipate the structure of our final result. The eigenvalues are entirely contained in the left-hand complex plane, and form a droplet which has maximal density close to the origin (see Fig.(1)). The hermiticity preserving property of the QDS generator $(\mathcal{L}[\rho])^{\dagger}=\mathcal{L}[\rho^{\dagger}]$ also implies that the eigenvalues are either real or come in complex conjugate pairs. For this reason, the density is symmetric under reflection across the real axis. The distribution of eigenvalues is highly reminiscent of that for real Markov generators Timm2009 ; Bordenave2014 , which are real Ginibre matrices whose columns sum to zero 444This similarity becomes essentially exact for a dissipator with $N^{2}-1$ jump operators, as shown recently in Denisov2018 .

The spectral gap exhibits a clear $N^{-2}$ scaling with Hilbert space dimension in Fig.(2), indicating that in the thermodynamic limit the spectral gap will close. We prove this assertion below.

It is not possible to describe the spectrum of (25) in a simple way in terms of eigenvalues or singular values of $L$, and we have not found an expression for the joint probability distribution function of its eigenvalues. Such shortcomings rule out the possibility of pursuing many classic approaches in RMT such as the method of orthogonal polynomials or the Coulomb plasma mapping. Instead, we appeal to another popular approach which involves developing a systematic large $N$ diagrammatic expansion of the moments.

The method of moments was first used by Wigner to derive the semi-circle distribution of eigenvalues for full random Hermitian matrices Wigner1955 , and later for k-body random matrix ensembles pandey1979binary ; mon1975statistical . The large $N$ diagrammatics for computing moments goes back to Refs. t1974planar ; brezin1978planar , which develop a topological expansion capturing the $N$ scaling of graphs by their Euler characteristic. Similar techniques arise in considering moments of more structured operators, such as k-body random matrices mon1975statistical , and propagators for systems with spatial locality constraints (for some recent examples see chan2018solution ; kukuljan2016corner ). Our ensembles of random Lindbladians with fullly random jump operators and Hamiltonians bear similarities to, but have important differences from, these previously studied models. While the planar approximation for traditional matrix models does not quite seem to hold due to the tensor product structure of Lindbladians, the lack of additional structure such as spatial locality allows for a simple revision of the classic topological expansion.

To begin, we discuss some basic symmetry properties of the moments. Expressed as a fourth-order tensor, the dissipator has the property

[TABLE]

This is a consequence of the dissipator being trace-preserving, and it can be seen directly from the representation

[TABLE]

Furthermore, this property holds for all matrix moments of the dissipator

[TABLE]

where the meaning of $\mathcal{D}^{n}$ is the $n$-fold composition of the map with itself. In the tensor notation, for instance,

[TABLE]

where summation over repeated indices here is implied. As a consequence of this symmetry, the propagator $K_{t}=e^{t\mathcal{D}}$ satisfies

[TABLE]

which is another statement of the conservation of probability, since it implies $\mbox{tr}(\rho(t))=\mbox{tr}(\rho(0))$. This constraint will be preserved upon ensemble-averaging over the jump operators. The most general form for the average matrix moment of the dissipator is

[TABLE]

The structure of the matrix moment will necessarily be a sum of products of Kronecker delta functions, due to Wick contraction of jump operators with the free indices $\{i,j,k,l\}$. The coefficients $A_{n}$, $B_{n}$, are $C_{n}$ are then combinatorial factors which depend non-trivially on $n$ and $N$, and are proportional to $v^{n}$.

First, we will show that the coefficient $C_{n}=0$. The matrix moment is the sum of many terms, all of the form

[TABLE]

where $X$ is a string of $L$ and $L^{\dagger}L$, and $Y$ is a string of $L^{\dagger}$ and $L^{\dagger}L$. Ensemble averaging over the complex Gaussian matrices $L$ requires Wick contraction between pairs $L$ and $L^{\dagger}$. Therefore, in order for a contraction to produce $\delta_{il}$, we need terms of the form $L_{iq}L_{pl}^{\dagger}$, or conversely $L_{qi}L_{lp}^{\dagger}$. However, the structure of the moments implies that these combinations are impossible. From (33), all we can hope to find are combinations $L_{iq}L^{\dagger}_{lp}$ or $L_{iq}^{\dagger}L_{lp}$. Therefore, it is impossible in general to find a contraction producing $\delta_{il}$ (an identical argument applies to $\delta_{kj}$) and we must have $C_{n}=0$. Note that these arguments fail for GOE jump operators, for which $C_{n}$ does not vanish.

Next, the constraint (29) implies that

[TABLE]

which is only possible if $A_{n}+NB_{n}=0$. Using this, we see that the averaged trace moments of the dissipator are given by $A_{n}$,

[TABLE]

The factor $N^{2}-1$ which appears in the denominator is the maximal number of non-vanishing eigenvalues of the dissipator, since there is a single unique zero mode Baumgartner2008 . A simple counting implies that for a variance which scales like $N^{-1}$ as in Eq.(26), the moment $A_{n}$ should be order $O(1)$ with $N^{-1}$ corrections. Thus, we find the general formula for the matrix moment in terms of the trace-moments

[TABLE]

As a consequence, the average propagator has the power series expansion

[TABLE]

where

[TABLE]

is the moment generating function. Accordingly, the disorder-averaged density matrix obeys

[TABLE]

This expression is somewhat misleading, since it seems to imply that the steady state is the uniform state proportional to the identity. However, it is known that the steady state is given by Baumgartner2008

[TABLE]

which is clearly not the uniform state. The discrepancy clearly comes from the order in which we do things. If we take an ensemble average before sending $t\to\infty$, then indeed the evolution tends to take the state to the infinite temperature, maximal entropy state. However, if we quench “disorder" and run $t\to\infty$, the nature of the steady state is far from uniform, and in fact appears to have relatively low entropy (since $M$ tends to have small eigenvalues, the steady-state density matrix will generically have a single large eigenvalue).

With the propagator in hand, we can calculate the dissipative form factor by taking the trace $\sum_{i,j}K_{t,ijij}$ to find

[TABLE]

Developing the generating function as an asymptotic series in $N^{-1}$, $G(t)=G_{0}(t)+N^{-1}G_{1}(t)+...$, we have in the large $N$ limit

[TABLE]

Our primary concern will be to understand the time evolution of $G_{0}(t)$. If this exhibits power-law decay, then we can say confidently that the eigenvalue spectrum of the dissipator is gapless. Having isolated the main object of study, we develop now the graphical calculus which is used to calculate the moments, and thus the generating function, to leading order.

We begin with a simple example - the mean-field dissipator. It can be evaluated exactly to yield

[TABLE]

The mean-field dissipator is the generator of a depolarizing channel Kastoryano2013 ,

[TABLE]

which has a spectral gap, and is thus rapidly mixing. We now show that taking fluctuations (higher moments) correctly into account will collapse the spectral gap. In order to do this, we introduce our diagrammatic thinking. First of all, since we are dealing with fourth-order tensors, we need four external legs. The tensor product structure of the dissipator in (25) allows us to organize these four external legs into two edges. The identity matrix in this notation is $\delta_{ik}\delta_{jl}$ and is represented by two parallel lines connecting $i$ to $k$ and $j$ to $l$. From this starting point, we indicate the insertion of a jump operator by a pair of colored dots, as in Fig. (3a), representing the indices of the matrix $L$. Then we have $L^{\dagger}L\otimes_{t}\mathbbm{1}$ represented in Fig. (3b); $L\otimes_{t}L^{\dagger}$ in Fig. (3c); and $\mathbbm{1}\otimes_{t}L^{\dagger}L$ represented in Fig. (3d).

In taking the expectation value, we must perform a Wick contraction to evaluate the correlators of complex Gaussian fields. This is accomplished by connecting red to black, and filled to open. Figure (4) shows how (44) would arise diagrammatically.

Equipped with these diagrammatics, we can compare the exact second moment of the dissipator with its diagrammatic representation. A direct, if tedious, calculation reveals that the second matrix moment is

[TABLE]

The $O(1)$ second-moment diagrams are shown in Fig.(5) (all second-order diagrams are presented in Fig.(15)). Since these must all be multiplied by the variance squared $(v/N)^{2}$, the infinite $N$ limit will preserve only diagrams of order $N^{2}$. This is only possible if under Wick contraction, the diagram includes two closed loops. It is clear that order $N^{2}$ terms in these second-order diagrams are only possible when the contractions occur on either the top or bottom edge, and not across the edges. Furthermore, such separated diagrams are not possible if there is any insertion of the recycling term $L\otimes L^{\dagger}$, for the simple reason that it only ever introduces an odd number of contractible Gaussian variables on a single edge, requiring a contraction that joins the edges.

Having found the structure of the leading order diagrams, we are left with the challenge of correctly counting them. This is accomplished if we simply remove the recycling term when calculating the moment. Note that for the second moment, we can recover the leading order contribution by looking at the simpler problem

[TABLE]

More generally, these considerations lead us to conclude that the leading order diagrams are completely accounted for by the “non-crossing trunctation" of the dissipator

[TABLE]

This concludes the proof of a special case of Lemma (1) for simple dissipators. We finish the proof in Sec.(3) for multiple jump operators (see also Appendix (C) for a more topologically oriented proof). Note that while our diagrammatic picture is exact, it does not reduce in any simple way to the familiar self-consistent Born approximation, an observation also made in Ref.Oganesyan2003 in attempting such an approach to calculate the self-energy for the Liouville-von Neumann operator. The culprit in their case and ours is clearly the double-legged structure of the “bare propagator", which spoils the traditional topological expansion in terms of planar diagrams t1974planar ; brezin1978planar ; bessis1980quantum . For more on the topological nature of the large $N$ limit see Appendix (C). Now, we proceed to evaluate the leading order form factor.

The non-crossing truncation in (48) significantly simplifies the problem, since now we can use the eigenvalues $u_{i}$ of $L^{\dagger}L$ (i.e. the squared singular values of $L$) to compute

[TABLE]

The distribution of these eigenvalues follows the Marčenko-Pastur law for complex Wishart matrices

[TABLE]

where we have $\int\rho_{w}(x)=1$. We also note here that this distribution follows from the resolvent

[TABLE]

which we will make use of later on. The leading order generating function then follows immediately, and is given by

[TABLE]

Along with (43), this proves our first proposition (1). The asymptotic behavior of the modified Bessel function implies that for long time, the exponential factor will disappear and the dynamics follows a power-law $F(t)\sim t^{-1}$. Since there is a unique steady state, the asymptotic limit of the form factor is $F(\infty)=N^{-2}$. Therefore, we expect the large $N$ behavior to be relevant on time scales $t^{-1}>>N^{-2}$, i.e. $t<<N^{2}$. However, finite size numerics indicates that a second regime takes over at times $t\sim O(N)$ interpolating between large $N$ behavior (52) and the asymptotic limit $N^{-2}$. This is an indication that the sub-leading $O(N^{-1})$ correction to $F(t)$.

2.2 Complex Random Normal Matrices

Here we consider a simple generator with the jump operator taken to be a random normal matrix satisfying $[L,L^{\dagger}]=0$. The eigenvectors of a general normal matrix form an orthonormal basis $\langle\tau_{i}|\tau_{j}\rangle=\delta_{ij}$ where $L|\tau_{i}\rangle=\tau_{i}|\tau_{i}\rangle$. The eigenmodes of the dissipator are then $|\tau_{i}\rangle\langle\tau_{j}|$ for all $i,j=1,...,N$, with eigenvalues

[TABLE]

Since $z_{ii}=0$ identically, there exist exactly $N$ zero eigenvalues. The mean density of eigenvalues is defined by the ensemble average of the density function

[TABLE]

Apart from the $N$ exact zero modes, the eigenvalues lie inside a shifted circle defined parametrically by

[TABLE]

where $(x,y)\in\mathbbm{R}^{2}$ are coordinates on the real plane (see Fig.(7)).

This demonstrates that the spectral gap is zero in the thermodynamic limit, which leads to the algebraic relaxation of Prop.(2). This result follows by direct computation, which is given in full detail in Appendix (D), and accomplished by utilizing the exact result for the $n$-point eigenvalue distribution functions

[TABLE]

where the kernel function is given by

[TABLE]

In particular, the dissipative form factor is

[TABLE]

and the two-point density correlation function is

[TABLE]

In the large $N$ limit, and setting $v=1$, the form factor behaves as

[TABLE]

The asymptotic expansion up to order $O(1)$ is presented in Appendix (D). The power-law relaxation as $t^{-2}$ is consistent with the closing of the spectral gap.

In Fig.(8), we compare the exact form factor with the leading order asymptotic formula, to illustrate the regime in which the leading order behavior is relevant. The form factor appears to saturate around $\tau\sim\sqrt{N}$, before which it is well described by the large $N$ behavior. Since the asymptotic limit is known to be $F(\infty)=1/N$, we expect that the asymptotic value takes over when $t^{-2}<<N^{-1}$, i.e. $t>\sqrt{N}$. Note that, unlike the simple dissipator with Ginibre jump operators, the asymptotic value enters at the same order as the sub-leading correction to $F(t)$. Consequently, the cross-over regime observed in Fig.(6) is not observed in Fig.(8).

2.3 Gaussian Unitary Ensemble

For Hermitian jump operators, the simple dissipator can be written as the square of the Liouville-von Neumann operator

[TABLE]

which acts on the density matrix as a nested commutator $\mathcal{D}\rho=-[L,[L,\rho]]$. We work exclusively with the Gaussian unitary ensemble, in which $L_{ij}\in\mathbbm{C}$ with mean zero and variance $\langle\bar{L}_{ij}L_{kl}\rangle=\frac{v}{N}\delta_{ik}\delta_{jl}$. As in (53), the eigenvalues of the dissipator $z_{ij}$ can be written in terms of the eigenvalues of the jump operators, and are given by

[TABLE]

now being real-valued since $\tau_{i}\in\mathbbm{R}$ are the eigenvalues of an hermitian $L$. As with (53), there exist $N$ exact zero eigenvalues of the dissipator, along with a distribution on the real line which can be extracted from the resolvent

[TABLE]

which in the large $N$ limit evaluates to

[TABLE]

with $N^{-1}$ corrections (see Appendix (E)). Obviously, this limit completely misses the $N$ zero eigenvalues, which will appear as a simple pole at $z=0$ with residue $N^{-1}$. The density of eigenvalues is obtained from the resolvent $\rho(x)=-(1/\pi){\rm Im}\,G(x+i0^{+})$, and plotted in Fig. (9).

These results follow directly from the method of moments Wigner1955 . Writing the resolvent as a Laurent series for large $z$ away from the real axis gives

[TABLE]

where we define the $k$-th trace moment of the GUE dissipator

[TABLE]

In the large $N$ limit, the two-point density correlation function can be approximated by the product of one-point densities. This approximation neglects the eigenvalue correlations. Defining $M_{k}=(-v)^{-k}\lim_{N\to\infty}\mu_{k}$, we have the explicit expression

[TABLE]

where the asymptotic mean density is the well-known semi-circle distribution

[TABLE]

The moments (68) clearly obey the recursion relation

[TABLE]

Inserting this into the Laurent series (66) and taking derivatives in order to massage the sum back to its original form, we arrive at the hypergeometric equation for the rescaled resolvent $F(\zeta)=-vG(1/\zeta)$,

[TABLE]

which is solved by (65).

The method of moments is also the most effective way to compute the dissipative form factor, which in this case can be compactly represented as as a hypergeometric function

[TABLE]

which relaxes algebraically at long times

[TABLE]

The power-law in this case is expected since the spectral gap closes in the thermodynamic limit. However, the rate of relaxation is slower than both the complex Ginibre and random normal dissipators above. The technical reason for this is the $1/\sqrt{x}$ (integrable) divergence of the spectral density close to the origin. To see this, we may use the Pfaff transformation (156) to find the resolvent (65) near the origin

[TABLE]

which closely resembles the resolvent for the Marčenko-Pastur distribution. The $1/\sqrt{x}$ divergence of the density thus becomes manifest, see Fig.(9).

If we repeat the arguments in the previous section to estimate the relaxation time in terms of the Hilbert space dimension, we expect the large $N$ form factor to be valid for times $t^{-1/2}\geq N^{-1}$, which gives $\tau\sim N^{2}$ as approximately the relaxation time after which the dissipative form factor can be expected to saturate, identical to the case with complex Ginibre jump operators. The same algebraic decay was found in Ref. Tameshtit1993 for the decay of purity under Lindblad evolution in which the jump operator and the Hamiltonian are identical and drawn from the Gaussian orthogonal ensemble.

3 Dissipators with Multiple Jump Operators

In this section, we are concerned with the proof of Proposition (3) for the long-time behavior of the propagator for the dissipator with multiple jump operators

[TABLE]

where $L_{a}$ are independent complex Ginibre matrices whose elements have zero mean and variance $\langle|L_{a,ij}|^{2}\rangle=v_{a}/N$. We assume $m$ is finite in the large $N$ limit. In Sec.(2.1), we proved Lemma (1) for a simple dissipator. Here, we extend the argument to multiple jump operators, completing the proof. Let us introduce the Hermitian matrix

[TABLE]

In terms of which the dissipator can be written

[TABLE]

The leading order term in the $k^{th}$ moment of $\mathcal{D}^{k}$ will be order $N^{k+2}$. This obviously requires the diagrams to be planar, since a contraction with intersecting lines will reduce the number of closed loops. However, according to our graphical calculus, this requirement is not enough. A planar diagram which has contractions that connect the two edges must necessarily be subleading. We can see this inductively. For a single crossing, the edges which accomplish the crossing cannot participate in a closed loop, making it impossible to produce $k$ loops. Thus, such diagrams are at most order $N^{k+1}$. If there are two crossing contractions, then a closed loop can form which connects the edges; however, this closed loop then must involve two contractions, making it also order $N^{k+1}$ at most. With three crossing contractions, we have $N^{k}$, since out of the three crossing contractions we can construct at most one closed loop. Proceeding, we see clearly that for $p$ crossing contractions on the $k$-th moment, the diagrams are at most order $N^{k+2-p+\lfloor p/2\rfloor}$, where $\lfloor p/2\rfloor$ is the largest integer $q\in\mathbbm{Z}$ such that $q\leq p/2$ (in other words, if $p$ is even, $q=p/2$, otherwise $q=(p-1)/2$). We conclude that the leading order diagrams are both planar and non-crossing. The insertion of any number of recycling terms will require a crossing contraction, and therefore be subleading to $N^{k+2}$. These arguments concern the topology of the diagrams which are possible for Ginibre jump operators, and are insensitive to the number of jump operators, and therefore apply to (75) for all $m$ (see Appendix (C) for an elaboration of this topological argument). The diagrammatics thus tells us that we will successfully capture the leading order trace moment by neglecting the recycling term and taking the $N\to\infty$ limit of the non-crossing truncation of the dissipator

[TABLE]

This completes the proof of Lemma (1) for the case of multiple jump operators. The limit which appears on the right-hand side is necessary, since the moments of the dissipator are not strictly equivalent to the moments of the non-crossing truncation of the dissipator at subleading order. As we emphasized before, this is essentially a trick which allows us to count the planar diagrams for which the two edges are disconnected.

As the sum of $m$ independent random Wishart matrices, we may find the eigenvalue density of $M$ using the methods of free probability theory. This requires first computing the self-energy $\Sigma$ (known as the R-transform in free probability theory, see Theorem 18 of Ref Mingo2017 ) of the single Wishart matrix $L^{\dagger}L$ (which we do not normalize by $N$). By inverting the resolvent Eq.(51), we find for the single Wishart matrix

[TABLE]

from which we find the self-energy by definition

[TABLE]

The self-energy for a sum of $m$ Wishart matrices with different $v_{a}$ is given by the sum of self-energies (which is the famous result that the $R$-transform is additive for free random variables)

[TABLE]

The resolvent for $M$ then follows implicitly by the equation

[TABLE]

We consider two scenarios for simplicity below.

Two Distinct Jump operators

Here we consider the case with two jump operators with variances $v_{1}=1$ and $v_{2}=\gamma$. The self-energy is

[TABLE]

The resulting resolvent has a spectral gap which scales linearly in $\gamma$ for small $\gamma$ and tends to a constant as $\gamma\to\infty$. This follows from an obvious duality in the eigenvalue density for $M=L_{1}^{\dagger}L_{1}+L_{2}^{\dagger}L_{2}$

[TABLE]

The small $\gamma$ gap follows easily from first order perturbation theory, and scales linearly in $\gamma$, which implies that $\rho_{\gamma}(x)=0$ for $x<\delta\sim\gamma$. Therefore, at large $\gamma$, the dual density $\rho_{1/\gamma}(x/\gamma)=0$ for $x/\gamma<\delta\sim 1/\gamma$. Therefore, $\rho_{\gamma}(x)$ vanishes for $x<\Delta\sim 1$ for large $\gamma$. A plot of the leading order eigenvalue density for various $\gamma$ is presented in Fig.(10).

The existence of a spectral gap in $M$ implies exponential decay of the dissipative form factor. We explore this in more detail in the next section.

Identical Jump Operators

Now we consider $m$ jump operators with the same variance $v_{a}=v$. The self-energy is

[TABLE]

which allows us to find a simple expression for the holomorphic resolvent

[TABLE]

The density of eigenvalues in terms of the dimensionless variable $\xi=x/v$ is

[TABLE]

The leading order contribution to the moment generating function for the dissipator is

[TABLE]

At large times, the integral is dominated by values of $ut\sim O(1)$. Defining the new variable $u^{\prime}=ut$, we get the asymptotic expression for large $t$

[TABLE]

which completes the proof of Prop.(3), and gives the explicit formula for the constant coefficient $C$.

4 Simple Lindbladian

In this section, we provide the proof of Proposition (4) and elaborate on the structure of the spectral gap. We consider here the simple Lindblad superoperator

[TABLE]

where $H$ is a complex hermitian random matrix drawn from the Gaussian unitary ensemble with zero mean and variance $\langle\bar{H}_{ij}H_{kl}\rangle=N^{-1}\delta_{ik}\delta_{jl}$, the jump operator $L$ is a complex Ginibre random matrix with variance $\langle\bar{L}_{ij}L_{kl}\rangle=N^{-1}\delta_{ik}\delta_{jl}$, and the dissipative coupling $\gamma>0$. In terms of the non-Hermitian Hamiltonian $K=H-i\gamma L^{\dagger}L$, the Lindblad superoperator takes the form

[TABLE]

We will refer to $\tilde{\mathcal{L}}$ as the non-unitary LvN generator (which is the non-crossing truncation per Def.(3)), while the second term is the recycling term. Lemma (1) implies that the moments of the Lindbladian are to leading order given by the diagrams with disconnected legs, which follows in a way identical to (78)

[TABLE]

One can construct a biorthogonal basis for $K$ using right-eigenvectors satisfying $K|R_{i}\rangle=z_{i}|R_{i}\rangle$, and left eigenvectors $\langle L_{i}|K=z_{i}\langle L_{i}|$, such that $\langle L_{i}|R_{j}\rangle=\delta_{ij}$. The conjugate transpose of the right eigenvectors satisfy $(K|R_{i}\rangle)^{\dagger}=\langle R_{i}|K^{\dagger}=\bar{z}_{i}\langle R_{i}|$. Using these facts, it is easy to see that the eigenmodes of the non-unitary LvN superoperator are given by $|R_{i}\rangle\langle R_{j}|$ with the spectrum

[TABLE]

From this, we find that the moment generating function for the simple Lindbladian is given to leading order by

[TABLE]

where the mean density for the non-hermitian Hamiltonian is defined as usual by

[TABLE]

The form factor (96) was also considered in Ref. Gudowska-nowak1998 , where it was interpreted as a survival probability for non-unitary LvN evolution. In Alhassid1998 , it was calculated for a related effective Hamiltonian with a lower-rank non-hermitian deformation in the weakly open limit, and related to the spectral autocorrelation function describing the induced photodissociation of molecules. Using the results in Appendix (F), the moment generating function for the density can be written

[TABLE]

where the boundary of the support of the density is described by the curve

[TABLE]

and the density contained in this curve is independent of $y$ and given by

[TABLE]

Note that we have switched to the convention used in Appendix (F), where we computed the density for the “rotated" matrix $iK$. Since the spectrum of $iK$ is symmetric about the $x$ axis, the spectral gap $\delta$ is given by the smallest real-valued $x_{min}$ which solves (99) when $y=0$. The resulting cubic equation for $x$ has explicit solutions, though their complicated form is not illuminating. In Fig.(12), $x_{min}(\gamma)$ is plotted revealing a non-monotonic function of the dissipative coupling. The spectral gap is known to close in both the $\gamma=0$ and $\gamma\to\infty$ limits, and obtains a maximal value when $\gamma=O(1)$. The asymptotic behavior of the gap can be easily determined in the small $\gamma$ limit from the solution to

[TABLE]

whereas for large $\gamma$, we can insert the scaling ansatz $x=\gamma^{-\alpha}$, and find

[TABLE]

To obtain the large $t$ asymptotics, we first perform the $y$ integration for the moment generating function

[TABLE]

where the lower limit is set by the spectral gap, while the upper limit is the largest real-valued solution to (99) with $y=0$. This has the obvious upper bound

[TABLE]

which leads directly to an upper bound for the form factor

[TABLE]

The integral can be evaluated directly using the exact form of the resonance width distribution (100)

[TABLE]

where the indefinite integral is given by

[TABLE]

At large time, we may use the definition of the exponential integral and its asymptotic behavior

[TABLE]

to find the leading behavior of $I_{b}(t)$ at long times

[TABLE]

This shows that $I_{x_{min}}(t)$, which has a slower decay rate, will be dominant at late times. Then plugging this into the expression (105) leads to the result stated in Prop.(4),

[TABLE]

where $\delta=x_{min}$ is the spectral gap of $K$. This result is strictly valid only for $\gamma\neq\{0,\infty\}$, since the limits $t\to\infty$ and $\gamma\to 0\,(\infty)$ do not commute in the expression for $I_{b}(t)$.

We compare our asymptotics and upper bound with the exact numerical form factor in Fig.(13). We observe that the exact form factor never seems to dip below its asymptotic value at any intermediate times, even though the large $N$ asymptotic formula (103) does hit zero. Increasing $N$ in Fig.(14) shows a clear trend of the first minimum in the form factor getting deeper with increasing $N$. This gives support to the conclusion that in the asymptotic limit, the form factor does vanish at finite time of order $\Delta^{-1}$ and determined solely by the dissipative coupling. The figure also indicates the regimes where the large $N$ asymptotics fails in finite-sized systems. This is approximately when the first subleading $O(N^{-1})$ correction takes over, and occurs before the leading asymptotic formula dips fully below $F(\infty)$. The asymptotic formula is an excellent fit to the exact numerical $F(t)$ for $F(t)>1/N$, and compares favorably for intermediate times when $N^{-1}>F(t)>N^{-2}=F(\infty)$.

5 Conclusion

In this paper we have introduced the dissipative form factor and determined its asymptotic long-time behavior for various ensembles of random Lindblad generators. The picture that emerges is that exponential relaxation to steady state is generic and to be expected when the dynamics is governed by multiple non-commuting operators. We see this in particular for a dissipator with multiple jump operators and a simple Lindbladian with a Hamiltonian and a single jump operator. We also find that random simple dissipators are gapless, and lead to algebraic relaxation of observables and correlation functions in the thermodynamic limit.

Our work is in the spirit of Wigner’s original proposal for a statistical description of complex atomic nuclei. Our setting extends this work to the time domain and asks how complex open systems equilibrate. An interesting follow-up to our work would be to characterize the structure of the steady state to which the random Lindblad equation flows, and how this depends on the spectral gap. A thermodynamic closing of the spectral gap has been implicated in second-order dissipative phase transitions Minganti2018 . In this context, it would be interesting to study the simple Ginibre dissipator which has a non-trivial steady state and is also gapless in the thermodynamic limit.

The closing of the spectral gap for the simple dissipator requires deeper understanding. Even though the evolution is dissipative for such a generator, the closing of the spectral gap indicates the existence, in the thermodynamic limit, of very long lived modes. This suggests that some sort of integrability emerges. The relaxation times involved in such systems (except for the random normal jump operators) implicate simple diffusion, though the lack of spatial structure perhaps produces more questions than answers along this line of reasoning.

As it stands, we believe the random Lindbladian can serve a benchmark for studies of integrability and thermalization in open quantum systems. In a closed chaotic quantum system, the spectral correlations are reproduced by random matrix theory, indicating universality in quantum chaotic systems. To understand how these features might come to bear on questions of universality of Lindblad evolution, we should understand how to interpret the jump operators. One possibility which we find promising is to invoke the eigenstate thermalization hypothesis (ETH) for thermalizing systems. In particular, it was originally argued that observables in closed thermalizing systems can be treated as random matrices Srednicki1999 . Since the jump operators can be constructed from the observables on the Hilbert space, the ETH would posit the jump operators as being described by random matrices in the energy eigenstate basis. This motivates our speculation that random Lindblad equations provide a universal dynamical description of open chaotic quantum systems. Presumably the signatures of universality will appear in the subleading contributions to the dissipative form factor, as they do for the spectral form factor for closed systems. This presents a promising direction for future work.

*Note Added - * While this paper was in preparation, the preprint Denisov2018 appeared which considered the spectrum of pure dissipators with the maximal $N^{2}-1$ number of jump operators. They similarly found a spectral gap, as well as an explicit expression for the limiting large $N$ distribution of eigenvalues.

Acknowledgements I thank my collaborators on can , out of which this work grew: S. Gopalakrishnan, V. Oganesyan and D. Orgad, with whom I have shared many fruitful discussions. The main results here and in Ref.can were announced on October 4, 2018 in the workshop “Random Matrices, Integrability, and Complex Systems" in Yad Hashmona, Israel. I have also benefited from discussions with A. Abanov, V. Albert, J. Feinberg, T. Seligman, W. Tarnowski, and V. K. Varma.

Appendix A Asymptotic Decay Rate and Spectral Gap

In this appendix, we provide an argument to support the connection between our results on the asymptotic decay rate and the existence of a spectral gap. For quenched “disorder", the spectral gap can be obtained formally from the propagator $K_{t}$ from the procedure

[TABLE]

This definition is sufficient for our purposes since the steady state is unique and purely imaginary eigenvalues occur with vanishing probability. Note that the form factor coincides with the ensemble average of $\hat{F}(t)$ in the large $N$ limit. The extreme value statistic would then require an ensemble averaging of this quantity

[TABLE]

We can only proceed by postulating that the infinite time limit commutes with ensemble averaging. This leaves us to evaluate $\langle\log\hat{F}\rangle$, which can be done formally by a replica-like limit

[TABLE]

The expectation value will now involve moments of the trace of the propagator. To understand the structure of this, let us consider $\langle\hat{F}^{2}\rangle$. Expanding the exponential, we see that this object will involve moments of the form

[TABLE]

Since a product of traces is involves, the basic starting point for constructing the diagrams for this will be four disconnected loops. In other words, replacing both instances of $\mathcal{L}$ with the identity produces four disconnected loops, each giving a factor of $N$. Insertions of the Hamiltonian or jump operators and subsequent Wick contraction will cut these loops and either join two loops, reducing the total order of the diagram, or split one loop into two, increasing the order. Following the arguments given in Appendix C, this implies that the leading order diagrams will be given by neglecting Wick contractions which bridge the two moments of $\mathcal{L}$.

We can put this in slightly different language, more familiar in the statistical mechanics setting. Let $\langle\rho(z_{1})\rho(z_{2})\rangle$ be the two-point eigenvalue density correlation function. Then

[TABLE]

In the large $N$ limit, our diagrammatic argument essentially boils down to the statement that the connected two-point function is subleading in $N$, i.e.

[TABLE]

Therefore, returning to the object of interest, we have that to leading order

[TABLE]

This suggests that the trace of the propagator is self-averaging. Note that now we have only to take the logarithm of the dissipative form factor. Therefore, we may extract the average spectral gap via

[TABLE]

Though not a rigorous proof, this argument supports the identification made in the paper between the asymptotic decay rate and the spectral gap.

Appendix B Complete Diagrams for Second Moment of the Simple Dissipator

The diagrams in (15) can be used to compute the second moment (46) exactly. These results can also be used to compute the third moment

[TABLE]

The leading order coefficient multiplying $(\delta_{ik}\delta_{jl}-N^{-1}\delta_{ij}\delta_{kl})$ can be expressed using moments of the Wishart matrix $L^{\dagger}L$ via

[TABLE]

which can be computed using the Marčenko-Pastur distribution (50)

[TABLE]

where $B(x,y)$ is the Euler beta function.

Appendix C Topological Expansion and Proof of Lemma (1)

The moments of the simple dissipator enjoy a topological expansion

[TABLE]

The factor $n_{b_{0}}$ is a combinatorial factor which counts the number of closed loops $b_{0}$, equal to the zeroth Betti number of the resulting graph. This coefficient must also account for the factors of $2$ or $-1$ coming from the different terms in the dissipator. Let us consider integer partitions of $n=p+q+r$. Then we may write this as

[TABLE]

which has the following meaning: a given partition $(p,q,r)$ accounts for the diagrams with $p$ insertions of $L\otimes_{t}L^{\dagger}$, $q$ insertions of $L^{\dagger}L\otimes_{t}\mathbbm{1}$, and $n-p-q$ insertions of $\mathbbm{1}\otimes_{t}L^{\dagger}L$. A particular partition will have a total number of diagrams given by

[TABLE]

However, since the insertions are not commutative, some of these can be distinct. This means there can be graphs with different topological index for the same partition. This is illustrated on a simple example in Fig.(17).

At a given order $n$, the number of connected components respects the bound $1\leq b_{0}\leq n+2$. The lower bound is assured since upon Wick contracting and taking the trace, one is guaranteed to have at least one connected component. The upper bound can be understood by the following argument: the identity “propagator" consists of no insertions, which means that upon taking the trace we have two closed loops. An insertion on one of the loops corresponds, topologically, to cutting the loop. Thus at order $n$, with $n$ insertions bilinear in $L$ and $L^{\dagger}$, before Wick contraction there will be a total of $2n$ connected components. The rules for Wick contraction can be visualized by inserting a short tube into each gap, and requiring tubes to meet. This is just a pictorial way to say that the indices of one matrix must be contracted with the indices of only one other matrix; hence the tubes must connect to tubes.

To understand the maximal possible $b_{0}$ for a given graph, we can proceed inductively and completely generally. The identity graph is two disconnected loops, and thus $b_{0}=2$. A single insertion will produce two cuts, which will either be on the same loop or on different loops. Thus, after Wick contraction, we may have $b_{0}=1$ or $b_{0}=3$. A second insertion will now produce two more cuts, which may occur on the same loop or on different loops. Repeating the argument for the single insertion, we end up with the possibilities $b_{0}\in\{2,4\}$. We may formally prove this by induction, but it is also clear by inspection that we can expect

[TABLE]

A new insertion can either join two loops, or split a single loop in two pieces. This recursive structure implies that to reach the maximal $b_{0}$ requires that one never join two loops. And since any number of insertions of the recycling term will result in a diagram that joins two loops, they do not participate in the leading order $N^{n+2}$ diagrams. Thus, we arrive at a proof of Lemma (1) which is completely agnostic about the number of jump operators in the problem.

The topological expansion we have outlined in this appendix seems distinct from the classical topological expansion proposed for large $N$ gauge theories t1974planar and matrix models brezin1978planar . We speculate that the difference lies in the fact that our random fields are fourth-order tensors (as opposed to matrices), due to the tensor product structure of the Lindbladian. Since the proof of Lemma (1) required nothing more sophisticated than Betti numbers, we did not seek a deeper connection to geometry by embedding our graphs in Riemann surfaces. However, this goal seems entirely attainable and could very well provide additional insights into the problem, possibly making the hermitian dissipator amenable to graphical enumeration techniques.

We should also comment here that another natural setting in which the tensor product structure will arise is when there exist (spatial) locality constraints on an operator (see e.g. chan2018solution ; kukuljan2016corner ). In this case, the spatial structure of the Hilbert space enforces the tensor product structure of operators (e.g. the Hamiltonian) acting in this space. It is possible that in this setting as well, one might benefit from the topological expansion discussed in this appendix.

Appendix D Random Normal Matrices: Correlation Functions and Exact Form Factor

The joint probability distribution function of eigenvalues for a random normal matrix with Gaussian distributed entries follows the same law as the complex Ginibre ensemble, and is consequently of determinantal form. It is given by Zabrodin2004

[TABLE]

where $Z_{N}$ is a normalization constant. All correlation functions follow from the kernel

[TABLE]

For instance, the density is given by

[TABLE]

This can be expressed more economically in terms of the incomplete Gamma function

[TABLE]

which in the $N\to\infty$ limit describes a droplet with constant density inside the disk of radius $R=\sqrt{v}$, and zero density outside. The two-point function is

[TABLE]

where

[TABLE]

The simple dissipator with random normal jump operators has eigenvalues

[TABLE]

In terms of the vectors on the 2D plane ${\bf r}_{1}$ and ${\bf r}_{2}$, this can be written

[TABLE]

A simple geometric argument using the fact that $|r_{i}|\leq\sqrt{v}$ implies that the support of the eigenvalues of the dissipator is described by a shifted circle law (see Fig.(7) )

[TABLE]

which clearly shows that the spectral gap vanishes.

The dissipative form factor is

[TABLE]

The exact result for finite $N$ (setting $v=1$) is

[TABLE]

For $t<<N$, we may develop this expression as an asymptotic series for large $N$,

[TABLE]

Curiously, this indicates that for times much smaller than $t\sim N$, there exists quadratic growth at order $O(1)$, reminiscent of the sub-leading growth of the spectral form factor in closed chaotic systems Berry1985 , and similarly due to eigenvalue correlations. The evidence for the closing of the spectral gap appears early on, however, at order $N^{2}$, which displays long-time power-law decay.

We split now the evaluation of the second integral in (139) into two parts using (133). The first is

[TABLE]

To get to the second expression in the first line, we use the fact that the density depends only on the absolute value of its argument. Upon expanding $\exp(2z_{1}\bar{z}_{2}t)$ in a power series, and performing the angular integration, only the first term in the expansion survives. The following steps follow by direct calculation using the explicit form of the density (130).

The second part of the integral in (139) is

[TABLE]

This is also evaluated by using a power series expansion of $\exp(2\tau x\bar{y})$ and noting that the angular integration will kill any terms with $p>q$, hence the restricted summation in the second line. The third line shows the terms in the integrand that do not vanish upon taking the angular integrals. The rest are straightforward calculations.

Combining (142) and (143) leads to the exact form factor (140).

Appendix E Random Liouville-von Neumann Generator

For easy reference, we discuss here the spectral properties of the Liouville-von Neumann operator when the Hamiltonian is a random matrix. These results are used in Sec.(2.3) to discuss the simple dissipator with GUE jump operators. The primary objects of interest are the resolvent and spectral form factor. We do not dwell on the details and refer to the literature for proofs and derivations (e.g. Mehta ; forrester2010log ).

Let $S=H\otimes_{t}\mathbbm{1}-\mathbbm{1}\otimes_{t}H$, with $H$ a random matrix drawn from the Gaussian unitary ensemble. Studying the statistical properties of this operator gives us access to both the Liouville-von Neumann operator $\mathcal{C}_{H}=-iS$ (5), and the simple dissipator with a hermitian jump operator $\mathcal{D}_{H}=-S^{2}$ (62). We fix the variance of the matrix elements to be independent of the symmetry class of the ensemble,

[TABLE]

The eigenvalue density function is

[TABLE]

where $\lambda_{i}$ are the eigenvalues of $H$. The expectation value of the density in the large $N$ limit is a smooth function given by the semi-circle law

[TABLE]

The variance $v$ determines the bandwidth to be $4\sqrt{v}$. The $n^{th}$ trace moment of the Liouville-von Neumann operator is defined as the expectation value of $\mbox{tr}\,S^{n}$. Since all odd powers will vanish due to the Gaussian statistics, we only have the even moments

[TABLE]

This can be conveniently expressed in terms of the two-point eigenvalue density correlation function

[TABLE]

The trace moments enjoy a large $N$ expansion whose leading order $O(1)$ term is given by using the disconnected two-point correlation function in (148). The asymptotic value of the moment $M_{2k}=\lim_{N\to\infty}\mu_{2k}$ is then

[TABLE]

We can confirm this by computing the first three exact moments at finite $N$

[TABLE]

Note the common factor of $(N^{2}-1)/N^{2}$ in these first few moments. This is due to the general structure of the ensemble averaged matrix moments discussed in Sec. (2.1).

The resolvent for the Liouville-von Neumann operator can be expressed as an asymptotic Laurent series for large $z$ with the coefficients given by the moments above

[TABLE]

whereas the resolvent for the hermitian simple dissipator is

[TABLE]

Substituting the exact large $N$ expression for the moments, we find that the resolvent can be written in terms of the Gauss hypergeometric function as

[TABLE]

with the dissipator resolvent quoted in Eq.(65). We present a proof of this in Sec.(2.3). The large $N$ limit here has the unfortunate feature that it completely misses the $N$ exact zero eigenvalues of both $\mathcal{C}_{H}$ and $\mathcal{D}_{H}$, since these would appear as a subleading correction $\frac{1}{Nz}$.

The Gauss hypergeometric function has singularities when its argument takes the values [math], $1$ and $\infty$. Since the power series is convergent only when the argument $16v/z^{2}$ is outside the unit circle, analytical continuation inside the unit disk is required. To this end, we can employ the Pfaff transformation of the hypergeometric functions Yoshida1997 ; Koepf1998

[TABLE]

which we use to evaluate the resolvent when $z$ is close to zero. To this end, let $a=-1/2$, $c=2$, $b=3/2$, and $x=16v/z^{2}$. Then we may use the Pfaff transformation to get that as $x\to\infty$,

[TABLE]

which implies for the resolvent (155)

[TABLE]

From this we see that the eigenvalue density has finite support in the domain $[-4\sqrt{v},4\sqrt{v}]$. The total length of the branch cut gives the bandwidth of the eigenvalue distribution, which in this case is $8\sqrt{v}$. As expected, this is twice the bandwidth of the spectrum of the Hamiltonian.

The leading order contribution to the spectral form factor of the Liouville-von Neumann operator can also be obtained by the method of moments, although a direct integration using the semi-circle law is much simpler. Comparing these, we find two representations of the spectral form factor

[TABLE]

The first expression, in terms of the Bessel function of the first kind, is very well known, while perhaps the latter, in terms of the generalized hypergeometric function, is less known. Using well-known asymptotics for the Bessel function one can verify the $t^{-3}$ power-law decay for the leading order contribution to the spectral form factor at long-times.

Appendix F Random Non-Hermitian Hamiltonian: Resolvent and Eigenvalue Density

Here we use the Hermitization trick Feinberg1997 ; Janik1997PRE (nicely reviewed in Feinberg2006 ; Janik1999a ) combined with a linearization procedure Belinschi2018 to obtain the resolvent for the non-Hermitian hamiltonian. Our results here reproduce those of Haake1992 ; Lehmann1995 , which were later significantly elaborated upon in Fyodorov1997a .

Let $H$ be a GUE matrix with variance $\langle|H_{ij}|^{2}\rangle=N^{-1}$, and $L$ be a complex Ginibre matrix with variance $\langle|L_{ij}|^{2}\rangle=\gamma/N$, where $N$ is the dimension of the matrices. Let $z$ be a complex coordinate (the argument of our Green’s function), and $\eta$ an auxiliary complex variable introduced to analytically continue the Green’s function “off" the complex plane Feinberg1997 . Consider the Green’s function obtained by inverting the following $4N\times 4N$ matrix

[TABLE]

Let $A=z+iH-L^{\dagger}L$, and define

[TABLE]

Then the elements of the full Green’s function here becomes

[TABLE]

We are ultimately interested in obtaining an expression for $G_{21}$, which in the $\eta\to 0$ limit becomes the resolvent $G_{21}=A^{-1}$ of the non-Hermitian Hamiltonian.

By symmetry arguments, we can see that the off-diagonal blocks of ${\bf G}$ will be zero after ensemble averaging. Since the matrix is linear in Gaussian fields, the self-consistent Born approximation, which consists of keeping only the second cumulant in the self-energy, must be exact in the large $N$ limit. In order to do this, we write

[TABLE]

where

[TABLE]

The self-energy functional is then

[TABLE]

This will involve expectation values of the following form

[TABLE]

Defining

[TABLE]

we wind up with the self-energy

[TABLE]

Note that each block is actually a matrix proportional to the $N\times N$ identity matrix; for ease of notation, we write only the proportionality constants, which are functions of $(z,\bar{z})$.

By inspection, we see that upon disorder averaging, we must have

[TABLE]

which simplifies the self-energy considerably

[TABLE]

When solving the Dyson equation for the Green’s function with this self-energy functional in the self-consistent Born approximation $\langle{\bf G}\rangle^{-1}={\bf G}_{0}^{-1}-\Sigma[\langle{\bf G}\rangle]$, we can take $\eta=0$ and invert the following matrix

[TABLE]

This implies

[TABLE]

which gives us

[TABLE]

Similarly, for the upper block we get the equations

[TABLE]

To find the solution inside the support of the eigenvalue density, we first consider the solution with $g_{11}\neq 0$. At this point, we have a straightforward, if extremely tedious, algebraic problem on our hands. Its solution can be found by hand or using symbolic manipulation software (e.g. Mathematica), resulting in the sought after trace of the resolvent

[TABLE]

from which follows the mean density

[TABLE]

which is independent of the horizontal direction $y$. The resolvent for the singular values $g_{11}$ can be written

[TABLE]

The boundary of the eigenvalue support is described parametrically by the condition $g_{11}=0$, giving

[TABLE]

or

[TABLE]

Appendix G Non-Hermitian Hamiltonian with Multiple Jump Operators

Here we provide an argument for the generalization of the spectral support curve (218) for a non-Hermitian Hamiltonian with multiple jump operators

[TABLE]

where $L_{a}$ are independent random complex Ginibre matrices whose elements have variance $1/N$. By following the steps taken in the replica approach used in Haake1992 , it is clear that the curve describing the support of the eigenvalues is simply modified to read

[TABLE]

In Haake1992 , $mN$ is the number of scattering channels, and $m<1$ is implicitly assumed in their work. However, this assumption does not appear necessary, and it is simple to see that it should generalize to our case in which $m\in\mathbbm{N}$ is an integer greater than one.

At small $\gamma$, we can drop the second term in brackets in (220), and solve

[TABLE]

The unique solution is $x_{min}=m\gamma$. However, the notable difference between $m=1$ is seen at large $\gamma$. This is perhaps not surprising given the fact that in the $\gamma\to\infty$ limit, the spectrum remains gapped, according to (87). To proceed, we can utilize the scaling ansatz again $x_{min}\sim x_{0}\gamma^{\alpha}$, and assuming $\alpha>0$, we are left to solve

[TABLE]

Self-consistency requires $\alpha+1={\rm min}(2,2\alpha)$, which is only compatible with $\alpha=1$. Now plugging this back in we can solve for the prefactor

[TABLE]

which has the solutions $x_{0}=(1\pm\sqrt{m})^{2}$. The correct solution comes by requiring $x_{0}$ vanish at $m=1$. Summarizing, we have for the asympototic behavior of the spectral gap of (219)

[TABLE]

Using that the spectral gap of the Lindbladian is $\langle\Delta\rangle=2x_{min}$, we recover the results found recently in sa2019spectral . Note that the spectral gap for $m>1$ appears to coincide for large $\gamma$ with the pure dissipator.