Spectral and Steady-State Properties of Random Liouvillians

TL;DR

This paper analyzes the spectral and steady-state behaviors of random Liouvillian operators in open quantum systems, revealing how dissipation strength influences spectral gaps, distributions, and scaling properties across different regimes.

Contribution

It introduces a comprehensive analysis of spectral and steady-state properties of stochastic Liouvillians, identifying regimes and scaling laws based on dissipation strength and system size.

Findings

Spectral gap increases with system size for multiple dissipation channels

Steady-state spectral distribution transitions from Poissonian to random matrix statistics with increasing dissipation

Different regimes exhibit distinct scaling exponents for spectral features

Abstract

We study generic open quantum systems with Markovian dissipation, focusing on a class of stochastic Liouvillian operators of Lindblad form with independent random dissipation channels (jump operators) and a random Hamiltonian. We establish that the global spectral features, the spectral gap, and the steady-state properties follow three different regimes as a function of the dissipation strength, whose boundaries depend on the particular quantity. Within each regime, we determine the scaling exponents with the dissipation strength and system size. We find that, for two or more dissipation channels, the spectral gap increases with the system size. The spectral distribution of the steady state is Poissonian at low dissipation strength and conforms to that of a random matrix once the dissipation is sufficiently strong. Our results can help to understand the long-time dynamics and steady-state properties of generic dissipative systems.

Tables

Table 1. Table 1: Scaling exponents for r = 1 𝑟 1 r=1 and r > 1 𝑟 1 r>1 , as defined in the main text. ν 𝜈 \nu exponents give the finite-size scaling of the spectral and steady-state quantities within one of the regimes P, C, D and their primed counterparts; λ 𝜆 \lambda exponents give the extra g eff subscript 𝑔 eff g_{\mathrm{eff}} -scaling besides g eff 2 superscript subscript 𝑔 eff 2 g_{\mathrm{eff}}^{2} in regime C and its primed counterparts; κ 𝜅 \kappa exponents give the finite-size scaling of the boundaries of the regimes.

$r$	${\nu }_P$	${\nu }_C$	${\nu }_D$	${\lambda }_X$		${\kappa }_X^{>}$	${\nu }_{P^{\prime }}$	${\nu }_{C^{\prime }}$	${\nu }_{D^{\prime }}$	${\lambda }_{\Delta }$		${\kappa }_{\Delta }^{>}$
$=1$	$0$	$1/2$	$1/2$	$2$	$-1/4$	$0$	$1/2$	$1/2$	$-3/2$	$-8/3$	0	$3/4$
$>1$	$0$	$1/2$	$1/2$	$2$	$-1/4$	$0$	$1/2$	$1/2$	$1/2$	—	$0$	$0$

Equations

$$\mathcal{L}={\mathcal{L}}_H+\sum _{\ell }{\mathcal{D}}_{W_{\ell }},$$

$${\mathcal{D}}_{W_{\ell }}\left(\rho \right)=W_{\ell }\rho W_{\ell }^{†}-\frac{1}{2}{W}_{\ell }^{†}{W}_{\ell }\rho -\frac{1}{2}\rho W_{\ell }^{†}{W}_{\ell }$$

$$\mathcal{L}\left(\rho \right)=-i\left[H,\rho \right]+g^2\sum _{j,k=1}^{N^2-1}{d}_{jk}\left\{G_j\rho G_k^{†}-\frac{1}{2}\left[\rho G_k^{†}{G}_j+G_k^{†}{G}_j\rho \right]\right\},$$

$$P_N\left(H\right)\propto \exp \left\{-\frac{1}{2}\mathrm{Tr}\left(H^2\right)\right\},$$

$$P_{\left(N^2-1,r\right)}\left(w\right)\propto \exp \left\{-\frac{1}{2}\mathrm{Tr}\left(w^{†}w\right)\right\}.$$

$$g_{\mathrm{eff}}={\left(2r\beta N\right)}^{1/4}g.$$

$$X\propto {\left[(\beta N{)}^{1/4}{g}_{\mathrm{eff}}\right]}^2{f}_X^{>}(g_{\mathrm{eff}}),$$

$$X\simeq g_{\mathrm{eff}}^2{f}_X^{<}\left[(\beta N{)}^{1/4}{g}_{\mathrm{eff}}\right],$$

$$\underset{g_{\mathrm{eff}}\to 0}{\lim }(\beta N{)}^{1/2}{f}_X^{>}(g_{\mathrm{eff}})\simeq \underset{(\beta N{)}^{1/4}{g}_{\mathrm{eff}}\to \infty }{\lim }{f}_X^{<}\left[(\beta N{)}^{1/4}{g}_{\mathrm{eff}}\right].$$

$${\varrho }_{\mathrm{I}}\left({\Lambda }_I\right)\simeq \int \mathrm{d}{E}_1\mathrm{d}{E}_2\delta \left({\Lambda }_I-E_1+E_2\right){\varrho }_{\mathrm{W}}\left(E_1\right){\varrho }_{\mathrm{W}}\left(E_2\right),$$

$${\sigma }_{{\rho }_0}^2\propto N^{-3}{f}_{{\rho }_0}^{<}\left(N^{1/2}{g}_{\mathrm{eff}}\right),$$

$${\sigma }_{{\rho }_0}^2\propto N^{-2}{f}_{{\rho }_0}^{>}(g_{\mathrm{eff}}),$$

$$Q\propto g_{\mathrm{eff}}^2(\beta N{)}^{{\nu }_Q^P}{f}_Q^{<}[(\beta N{)}^{-{\kappa }_Q^{<}}{g}_{\mathrm{eff}}],$$

$$Q\propto g_{\mathrm{eff}}^2(\beta N{)}^{{\nu }_Q^D}{f}_Q^{>}[(\beta N{)}^{-{\kappa }_Q^{>}}{g}_{\mathrm{eff}}].$$

$$f_Q^{<}(x)\propto x^0\mathrm{if}x\to 0,$$

$$f_Q^{<}(x)\propto x^{{\lambda }_Q^{<}}\mathrm{if}x\to \infty ,$$

$$f_Q^{>}(x)\propto x^0\mathrm{if}x\to \infty ,$$

$$f_Q^{>}(x)\propto x^{{\lambda }_Q^{>}}\mathrm{if}x\to 0.$$

$$\underset{g_{\mathrm{eff}}(\beta N{)}^{-{\kappa }_Q^{<}}\to 0}{\lim }(\beta N{)}^{{\nu }_Q^P}{f}_Q^{<}[(\beta N{)}^{-{\kappa }_Q^{<}}{g}_{\mathrm{eff}}]\simeq \underset{g_{\mathrm{eff}}(\beta N{)}^{-{\kappa }_Q^{>}}\to \infty }{\lim }(\beta N{)}^{{\nu }_Q^D}{f}_Q^{>}[(\beta N{)}^{-{\kappa }_Q^{>}}{g}_{\mathrm{eff}}],$$

$$g_{\mathrm{eff}}^{2+{\lambda }_Q^{<}}{\left(\beta N\right)}^{{\nu }_Q^P-{\kappa }_Q^{<}{\lambda }_Q^{<}}=g_{\mathrm{eff}}^{2+{\lambda }_Q^{>}}{\left(\beta N\right)}^{{\nu }_Q^D-{\kappa }_Q^{>}{\lambda }_Q^{>}}$$

$${\lambda }_Q\left({\kappa }_Q^{>}-{\kappa }_Q^{<}\right)={\nu }_Q^D-{\nu }_Q^P.$$

$$G_{\mathcal{L}}(z)=\sum _{n=0}^{\infty }\frac{1}{z^{n+1}}{\left(-\frac{1}{2}\right)}^2\sum _{k=0}^n\left(\genfrac{}{}{0ex}{}{n}{k}\right)C_k{C}_{n-k},$$

$${\varrho }_{\mathcal{L}}(x)=2{\int }_{{\nu }_m}^{{\nu }_M}\mathrm{d}\nu {\varrho }_{\mathrm{MP}}(\nu ){\varrho }_{\mathrm{MP}}(-\nu -2x),$$

$${\varrho }_{\Gamma }(x)=\frac{1}{2\pi x}\sqrt{({\xi }_{+}-x)(x-{\xi }_{-})},{\xi }_{-}<x<{\xi }_{+},$$

$$⟨\Delta ⟩=\beta N{g}^2{\xi }_{-}=\beta Nr{g}^2=\sqrt{\beta N}{g}_{\mathrm{eff}}^2\frac{(1-\sqrt{r}{)}^2}{\sqrt{2r}},$$

$$X\propto \beta N{g}^2({\xi }_{+}-{\xi }_{-})=4\beta N\sqrt{r}{g}^2\propto \sqrt{\beta N}{g}_{\mathrm{eff}}^2,$$

$$A_{nm}\equiv\,\left\langle\left\langle nn\right|\right|\sum_{\ell}\mathcal{D}_{W_{\ell}}\left|\left|mm\right\rangle\right\rangle=\left\{\begin{array}[]{r@{\quad}cr}\sum_{\ell}W^{(\ell)}_{nm}{W^{(\ell)}_{nm}}^{*},&\mathrm{if}&m\neq n\,,\\ -\sum_{\ell}\sum_{k\neq m}W^{(\ell)}_{km}{W^{(\ell)}_{km}}^{*},&\mathrm{if}&m=n\,,\end{array}\right.$$

$$A_{nm}\in \mathbb{R}(\mathrm{for}\mathrm{all}n),$$

$$A_{nm}>0(n\ne m),$$

$$\sum _{n=1}^N{A}_{nm}=0(\mathrm{for}\mathrm{all}n),$$

Full text

1 Introduction

Although open quantum systems are strongly influenced by their environments, a description of all environment degrees of freedom is often impossible and unnecessary. Master equations describing the evolution of the system’s reduced density matrix, $\partial_{t}\rho=\mathcal{L\left(\rho\right)}$, where the Liouvillian acquires the Lindblad form [1]

[TABLE]

provide a simple approach to model open quantum dynamics. Here, $\mathcal{L}_{H}\left(\rho\right)=-i\left[H,\rho\right]$ represents the unitary evolution under the Hamiltonian $H$, and

[TABLE]

the contribution of each dissipation channel by the action of the jump operator $W_{\ell}$.

This approach assumes either a weakly coupled or strongly coupled environment [1], with memory times much shorter than all other characteristic energy scales. Due to this so-called Markovian assumption, Lindblad dynamics fails to capture certain processes, such as the ones responsible for coherent low-temperature transport [2]. Nevertheless, the method is widely used in areas ranging from thermodynamics of quantum engines [3] to the description of quark-gluon plasma [4]. Perhaps, its most important application is to model quantum optic setups, which typically fulfill the required conditions as driving fields shift the relevant energies of the environment to a region where the spectral density is large [5].

Although the Lindblad form significantly simplifies the problem, obtaining dynamic and steady-state properties given a Lindblad operator remains a major theoretical challenge. In one dimension, efficient numerical approaches have been developed based on matrix product operator ideas [6, 7]. An alternative strategy, that has recently been extremely successful, consists of studying exactly solvable (or integrable) models [8, 9, 10, 11, 12, 13]. Albeit enlightening, integrable Lindbladians, as their Hamiltonian counterparts, are expected to have very peculiar properties and remain a set of measure zero among all possible Lindbladian dynamics. For generic cases, exact diagonalisation of relatively small systems remains the only option.

Tools to study generic Hamiltonians are available for non-integrable closed quantum systems. They rely on the widely supported conjecture [14] that universal features of spectral and eigenstate properties of quantum systems with a well-defined chaotic classical limit follow those of random matrix theory (RMT) [15, 16, 17, 18, 19, 20]. This is to be contrasted with integrable systems where level spacings typically follow Poisson statistics [21]. RMT has been reported to also describe complex many-body systems without well defined classical correspondents [22, 23, 24]. These predictions are particularly appealing as they are insensitive to the microscopic details of particular models, relying solely on symmetry properties of the Hamiltonian, i.e. to which: Gaussian unitary (GUE), Gaussian orthogonal (GOE), or Gaussian symplectic (GSE) ensemble, it belongs to. In view of the tremendous success of RMT, it is natural to ask if a similar approach can be followed in the case of Lindbladian dynamics. While open chaotic scattering [25, 26, 27] and quantum dissipation and decoherence [28, 29, 30, 31, 32, 33] have been studied using RMT in the past, this route has remained essentially unexplored up to very recently [34, 35, 36].

Following these general ideas,  [34] studies an ensemble of random Lindblad operators with a maximal number of independent decay channels and finds that, in this case, the spectrum acquires a universal lemon-shaped form. In a similar spirit, [35] considers an ensemble of random Lindbladians consisting of a Hamiltonian part and a finite number of Hermitian jump operators and finds a sharp spectral transition as a function of the dissipation strength. Finally, [36] studies in detail the spectral gap for several different Liouvillians. Yet, there is a number of open questions related to the nature of the spectrum, specifically when the dissipative and the Hamiltonian components are comparable. Additionally, when the jump operators are Hermitian, the steady state is the infinite-temperature thermal state (i.e. proportional to the identity operator). The properties of the non-trivial steady state, ensuing in the presence of non-Hermitian jump operators, are completely unexplored.

The aim of this paper is to study the universal properties of the spectrum and the steady state of Liouvillian operators consisting of a Hamiltonian component and a set of dissipation channels, taken from appropriate ensembles of random matrices. We consider a finite number of non-Hermitian jump operators allowing for a non-trivial steady state. As a function of the effective dissipation strength, $g_{\mathrm{eff}}$, we find a rich set of regimes regarding global spectral features such as the spreading of the decay rates $X$, the spectral gap $\Delta$, and the spectral properties of the steady-state density matrix $\rho_{0}$. Each regime is characterised by a scaling of the corresponding quantity with the system size $N$ and with $g_{\mathrm{eff}}$. The finite-size scaling of the boundaries between different regimes is also determined.

We expect these results to model a very broad class of nonregular Markovian open quantum systems. Indeed, in view of the quantum chaos conjecture, systems effectively modeled by random Lindbladians should be the rule rather than the exception. The predictive power of the present analysis relies on the fact that, once a regime is determined, the system’s universal properties can be obtained solely by symmetry arguments.

The paper is organised as follows. First, we describe an unbiased construction of a random Liouvillian in Section 2. Second, we numerically study its spectral (Section 3) and steady-state (Section 4) properties. Finally, we end with a short summary of our findings and their implications for determining the properties of generic Markovian dissipative systems in Section 5. Several appendices contain details of analytical derivations and additional numerical results.

2 Random Liouvillian Ensembles

To obtain a suitable set of random Liouvillians of the Lindblad form, we define a complete orthogonal basis, $\left\{G_{i}\right\}$ with $i=0,\dots,N^{2}-1$, for the space of operators acting on an Hilbert space of dimension $N$111 We consider only systems with finite-dimensional local Hilbert spaces. Below, we are interested in taking the thermodynamic limit, which is to be understood as $N\to\infty$., respecting $\mathrm{Tr}\left[G_{i}^{\dagger}G_{j}\right]=\delta_{ij}$, with $G_{0}=\mathbbm{1}/\sqrt{N}$ proportional to the identity. Each jump operator can be decomposed as $W_{\ell}=g\sum_{j=1}^{N^{2}-1}G_{j}w_{j\ell}$. Note that $W_{\ell}$ is taken to be traceless, i.e. orthogonal to $G_{0}$, to ensure that the dissipative term in (1) does not contribute to the Hamiltonian dynamics. In the $\left\{G_{i}\right\}$ basis, the Liouvillian is completely determined by matrices $H$ and $w$,

[TABLE]

where $d_{jk}=\sum_{\ell=1}^{r}w_{j\ell}w_{k\ell}^{*}=(ww^{\dagger})_{jk}$ is an $(N^{2}-1)\times(N^{2}-1)$ positive-definite matrix. We denote by $r$ the number of jump operators in (1) (i.e. $\ell=1,\dots,r$) which counts the number of independent system operators coupled to independent environment degrees of freedom. To obtain a random Liouvillian, we draw $H$ from a Gaussian ensemble with unit variance [37, 18], i.e.

[TABLE]

and $w$ from a Ginibre ensemble [38, 18]222For each realisation of the system, our operators are time-independent. The statistical approach then corresponds to ensemble-averaging over random matrices. This is to be contrasted with the different approach of time-averaging perturbations that evolve randomly in time, but have a fixed direction in matrix-space., i.e.

[TABLE]

The coupling constant $g>0$ parameterizes the dissipation strength (when the typical scale of the frequencies of $H$ is of order one, then the typical scale of the decay rate of $\mathcal{D}_{W_{\ell}}$ is of order $g^{2}$). We consider two cases: real matrices, $H_{ij}=H_{ji},w_{ij}\in\mathbb{R}$ (labeled by $\beta=1$ in the following), and complex matrices, $H_{ij}=H_{ji}^{*},w_{ij}\in\mathbb{C}$ ($\beta=2$).

We study the statistical properties of the spectrum of $\mathcal{L}$ drawn from an ensemble of random Lindblad operators parameterised by $N$, $r$, $\beta$ and $g$. The right eigenvectors of $\mathcal{L}$, respecting $\mathcal{L}\left(\rho_{\alpha}\right)=\Lambda_{\alpha}\rho_{\alpha}$, with $\alpha=0,\dots,N^{2}-1$, are denoted by $\rho_{\alpha}$, with $\Lambda_{\alpha}$ the respective eigenvalue. By construction, $\mathrm{Re}\left(\Lambda_{\alpha}\right)\leq\Lambda_{0}=0$, and $\rho_{0}$, if unique333The steady state is unique in the absence of any additional symmetries. Since we are considering the less structured Liouvillian possible the steady states we find are unique by construction., is the asymptotic steady state, left invariant by the evolution. If an experiment probes a dissipative system for long-enough timescales, it is this state it studies.

In Figure 1 $(b)$, we represent several relevant energy- (or inverse time-) scales. The spectral gap, $\Delta=\min_{\alpha>0}\mathrm{Re}\left(-\Lambda_{\alpha}\right)$ describes the typical time it takes the system to reach the steady state, i.e. the duration of transient effects. The variance along the real axis, $X^{2}=\sum_{\alpha}\left[\mathrm{Re}\left(\Lambda_{\alpha}-R\right)\right]^{2}/N^{2}$, sets the spread of decay rates and also the typical minimum time to observe the onset of dissipation and decoherence. The variance along the imaginary axis, $Y^{2}=\sum_{\alpha}\left(\mathrm{Im}\Lambda_{\alpha}\right)^{2}/N^{2}$, gives the timescale for the oscillations of the states’ phases. We also depict the center of mass of the spectrum, $R=\sum_{\alpha}\Lambda_{\alpha}/N^{2}$, which can usually be trivially shifted away.

The spectrum and eigenvectors of $\mathcal{L}$ are obtained by exact diagonalisation.

3 Spectral properties

Figures 1 $(a)$–$(c)$ show the spectrum of a random Lindblad operator in the complex plane computed for different values of $g$. The boundaries of the spectrum evolve from an ellipse, for small $g$, to a lemon-like shape at large $g$. In Figure 1 $(c)$ we plot (solid line) the spectral boundary, explicitly computed in [34] for the case $r=N^{2}-1$444 The spectral boundary from [34] is centered at the origin and has to be displaced by a shift $-(1+r)/(2\sqrt{r})$ and then rescaled by $\sqrt{2}$ to give the line of Figure 1 $(c)$.. These results, obtained here for $r=2$, indicate that the lemon-shaped spectral boundary is ubiquitous in the strong dissipation regime.

Figures 1 $(d)$ and $(e)$ show the scaling of $X$ and $Y$ with $N$ as a function of

[TABLE]

The three representative points (P, C, D) of Figure 1 $(d)$ correspond to regimes for which $X$ depicts a qualitatively different behaviour. The perturbative regime P corresponds to the weak-coupling limit of the Lindblad equation, while the dissipative regime D relates to the singular-coupling limit. The observed scaling collapse shows that, for large $g_{\mathrm{eff}}$, the standard deviation along the real axis behaves as

[TABLE]

with $f_{X}^{>}$ an unknown scaling function satisfying $f_{X}^{>}(x)\propto x^{0}$ for $x\!\to\!\infty$ and $f_{X}^{>}(x)\!\propto\!x^{2}$ for $x\to 0$. A similar scaling collapse can be obtained for small $g_{\mathrm{eff}}$ (not shown, see B) yielding

[TABLE]

for which the scaling function satisfies $f_{X}^{<}(x)\propto x^{0}$ when $x\!\to\!0$ and $f_{X}^{<}(x)\propto x^{2}$ when $x\!\to\!\infty$. The asymptotic power law behaviour of the scaling functions is hard to determine for the available values of $N\leq 120$, however, it is fully compatible with the large-$N$ extrapolation both in the large- and small-$g_{\mathrm{eff}}$ regimes (see black crosses in Figure 1 $(d)$ and also B). This result is further corroborated by the compatibility conditions (for further details see A)

[TABLE]

We can thus identify the three regimes: P, for $g_{\mathrm{eff}}\lesssim\left(\beta N\right){}^{-1/4}$; C, for $(\beta N)^{-1/4}\lesssim g_{\mathrm{eff}}\lesssim(\beta N)^{0}$; and D, for $(\beta N)^{0}\lesssim g_{\mathrm{eff}}$, corresponding to each representative point. Note that the rescaled quantities plotted in Figure 1 $(d)$ and $(e)$ do not seem to depend on the index $\beta$. There is a small $r$-dependence, which converges rapidly for increasing $r$ (see insets). For the standard deviation along the imaginary axes we find, $Y\simeq\sqrt{\beta N}$ for regimes P and C, and $Y\propto\sqrt{\beta N}g_{\mathrm{eff}}^{2}$ within regime D. Therefore, the variance in regime P can be explained by a perturbative treatment of the dissipative term—the value of $Y$ corresponds to the unit variance of the random Hamiltonian, and $X\propto g_{\mathrm{eff}}^{2}$ is expected from a degenerate perturbation theory.

Figures 2 $(a)$–$(f)$ and $(g)$–$(l)$ show the spectral density, $\varrho\left(\Lambda\right)=\sum_{\alpha}\delta^{2}\left(\Lambda-\Lambda_{\alpha}\right)$, along the cuts depicted in Figure 1 $(c)$ for regimes P and D, respectively, for various numbers of jump operators. Results for $\beta=1$ show the same limiting behaviour in the large-$N$ limit, after proper rescaling. For P, the density along $\mathrm{Re}\left(\Lambda\right)=R$ and $R+X$, does not depend on the number of jump operators, $r$. The integrated distribution of the imaginary parts, $\varrho_{\mathrm{I}}\left(\Lambda_{\mathrm{I}}\right)=\int\mathop{}\!\mathrm{d}^{2}\Lambda\,\delta\!\left(\Lambda_{\mathrm{I}}-\mathrm{Im}\Lambda\right)\varrho\left(\Lambda\right)$, is well described by a convolution of Wigner’s semicircle laws, $\varrho_{\mathrm{W}}(E)\simeq\sqrt{E_{*}^{2}-E^{2}}$, i.e.

[TABLE]

except at $\Lambda_{\mathrm{I}}=0$ where there is an increase of spectral weight, which is depleted from the immediate vicinity of the real axis. This last feature holds for general hermiticity-preserving operators, which can be brought to a real representation by a trivial similarity transformation. They are therefore related to the real Ginibre ensemble, where it can be explicitly shown [39]. The spectral weights along the cuts $\mathrm{Im}\left(\Lambda\right)=0$ and $Y$ scale with $\sqrt{r}$ and the integrated distribution of the real parts, $\varrho_{\mathrm{R}}\left(\Lambda_{\mathrm{R}}\right)=\int\mathop{}\!\mathrm{d}^{2}\Lambda\,\delta\!\left(\Lambda_{\mathrm{R}}-\mathrm{Re}\Lambda\right)\varrho\left(\Lambda\right)$, is well approximated by a Gaussian. As for the variance, in the P regime, these results can be derived from a perturbative small-$g$ expansion. In the strongly dissipative regime, D, the spectral density along $\mathrm{Re}\left(\Lambda\right)=R$ and $R+X$ depends on $r$ but converges rapidly to the $r\to\infty$ limit.

The case $r=1$ is qualitatively different from $r>1$. This can also be observed in the $\mathrm{Im}\left(\Lambda\right)=0$ and $Y$ cuts, and in $\varrho_{\mathrm{R}}$, where for $r=1$ the spectral weight is finite for $\Lambda_{\mathrm{R}}\to 0^{-}$, in contrast with the $r>1$ results that develop a spectral gap. We also found different scaling properties for $r=1$ compared to $r>1$, see below. Presently, the physical reason behind this discrepancy between $r=1$ and $r>1$ is unknown.

We now turn to the study of $\Delta$, which is a particularly important spectral feature since it determines the long-time relaxation asymptotics. Figure 3 shows the average spectral gap, $\langle\Delta\rangle$, for $\beta=2$ and $r=2$ as a function of $g_{\mathrm{eff}}$ for different values of $N$. In C, we checked that in the thermodynamic limit the distribution of the gap becomes sharply peaked around its mean; hence, the latter accurately describes the long-time dynamics.

Here, we also find three qualitatively different regimes (P*′, C′, D′) whose boundaries do not coincide with those in Figure 1. Here the boundaries $g_{\mathrm{eff}}\simeq g_{\mathrm{P^{\prime}C^{\prime}}}$ and $g_{\mathrm{eff}}\simeq g_{\mathrm{C^{\prime}D^{\prime}}}$, separating the P′* and C*′, and the C′* and D*′* regimes, respectively, are independent of $N$ for $N\to\infty$. For regime P*′, the average gap behaves as $\langle\Delta\rangle\propto(\beta N)^{1/2}g_{\mathrm{eff}}^{2}$; for C′, the gap varies as $\langle\Delta\rangle\propto(\beta N)^{1/2}g_{\mathrm{eff}}$; for D′*, we observe again $\langle\Delta\rangle\propto(\beta N)^{1/2}g_{\mathrm{eff}}^{2}$. As before, these results are only possible to establish by performing a large-$N$ extrapolation of the available data due to the presence of large finite-size corrections (see C).

The limiting $N\to\infty$ values of the spectral gap for small and large $g_{\mathrm{eff}}$ can be determined either by perturbative arguments or by computing the holomorphic Green’s functions, respectively (see C for details; see also [40, 36] for a computation for arbitrary $g_{\mathrm{eff}}$). One finds that $\langle\Delta\rangle=\beta Ng^{2}(1-\sqrt{r})^{2}$ for large $g_{\mathrm{eff}}$, which is compatible with [36]. For small $g_{\mathrm{eff}}$, employing degenerate perturbation theory to generalize the results of [41] (see C for details), we find that for large $r$ the average gap is given by $\langle\Delta\rangle=\beta Ng^{2}r$. The predictions describe the gap increasingly well for growing $r$; for $r=2$, although not exact, they give a good estimate, see inset in Figure 3. Notwithstanding that the above two results are derived in the limits $g_{\mathrm{eff}}\to 0$ and $g_{\mathrm{eff}}\to\infty$, they provide a remarkable description for the whole P*′* and D*′* regimes, respectively. For the special case $r=1$, although three regimes are also present (see E), the scaling of the gap with $N$ changes in the strongly dissipative regime. Finally, none of the regimes above show the mid-gap states reported in [35], where only Hermitian jump operators were considered.

4 Steady-state properties

We next characterize the steady state $\rho_{0}$. First, we consider the variance, $\sigma_{\rho_{0}}^{2}$, of the eigenvalues of $\rho_{0}$. This quantity is related to the difference between the purity of the steady state, $\mathcal{P}_{0}=\mathrm{Tr}\left(\rho_{0}^{2}\right)$, which quantifies the degree of mixing of $\rho_{0}$, and that of a fully-mixed state $\mathcal{P}_{\mathrm{FM}}=1/N$, $\mathcal{P}_{0}-\mathcal{P}_{\mathrm{FM}}=N\sigma_{\rho_{0}}^{2}$. Figure 4 $(a)$ shows the variance $\sigma_{\rho_{0}}^{2}$ as a function of $g_{\mathrm{eff}}$ for $\beta=2$ and $r=2$. Here again, three different regimes can be observed whose boundaries do not coincide with those given for previous quantities. In regime P*′′*, $g_{\mathrm{eff}}\lesssim(\beta N)^{-\frac{1}{2}}$, we observe

[TABLE]

with $f_{\rho_{0}}^{<}(x)\propto x^{0}$, for $x\!\to\!0$, and $f_{\rho_{0}}^{<}(x)\!\propto\!x^{2}$ for $x\!\to\!\infty$. In regime D*′′*, for $N^{0}\lesssim g_{\mathrm{eff}}$, we have

[TABLE]

with $f_{\rho_{0}}^{>}(x)\!\propto\!x^{0}$ for $x\!\to\!\infty$ and $f_{\rho_{0}}^{<}(x)\!\propto\!x^{2}$ for $x\!\to\!0$. The crossover regime C*′′* can be accessed by both asymptotic expansions and corresponds to $\sigma_{\rho_{0}}^{2}\propto N^{-2}g_{\mathrm{eff}}^{2}$. The asymptotic matching again agrees with data extrapolation (see D). These scalings imply that, up to subleading $1/N$ corrections, the steady state is fully mixed in regime P*′′, $\langle\mathcal{P}_{0}\rangle=\mathcal{P}_{\mathrm{FM}}+\mathcal{O}\left(1/N^{2}\right)$, while in regimes C′′* and D*′′*, $\mathcal{P}_{0}$ is only proportional to $\mathcal{P}_{\mathrm{FM}}$. At large dissipation, the steady state can be very well described by a random Wishart matrix, in agreement with general results of the entanglement spectrum of random bipartite systems [42, 43]. Purity can then be computed straightforwardly. This will be discussed for a more tractable model in a future publication [44]. The case with $r=1$ is again qualitatively different from $r>1$ at strong dissipation, the purity $\langle\mathcal{P}_{0}\rangle$ being of order $N^{0}$ (see E), signaling a steady state closer to a pure state than to a fully-mixed one.

Next, we investigate the steady state’s spectrum. For that, it is useful to introduce the effective Hamiltonian $\mathcal{H}=-\log\rho_{0}$. We denote the eigenvalues of $\mathcal{H}$ by $\varepsilon_{i}$ and present their spectral density $\varrho(\varepsilon)$ in Figures 4 $(b)$ and $(c)$, corresponding to a point in regime P*′′* and D*′′, respectively. In the weak coupling regime, P′′, $\varrho(\varepsilon)$ is well described by a Gaussian, while at strong coupling, D′′, it acquires a non-Gaussian shape. For large $r$, $\varrho(\varepsilon)$ is increasingly well described by a Wigner semicircle distribution. A remarkable agreement can already be seen for $r=10$ in the example of Figure 4 $(c)$. For small $r$ (see $r=2$ in Figure 4 $(c)$) there is a systematic skewing of the spectrum to the right. Figure 4 $(d)$ presents the probability distribution, $P(r)$, of adjacent spacing ratios, $r_{i}=s_{i}/s_{i-1}$, with $s_{i}=\varepsilon_{i+1}-\varepsilon_{i}$, which automatically unfolds the spectrum of $\mathcal{H}$ [45]. The analytic predictions for the GUE and for the Poisson distribution [45] (full lines) are given for comparison. The agreement of the numerical data of the points in the P′′* and D*′′* regimes, respectively, with the Poisson and GUE predictions is remarkable. Within regime C*′′*, we observe a crossover between these two regimes.

To illustrate the crossover in the spectral properties of $\mathcal{H}$ with $g_{\mathrm{eff}}$, we provide in Figure 4 $(e)$ the ratio of the first two moments of the distribution of $r$, $(\sigma_{r}/\langle r\rangle)^{-1}$ which can distinguish between Poisson and GUE statistics. Since in the Poissonian case $P(r)=1/(1+r)^{2}$, the $n$-th moment of the distribution diverges faster than the $(n-1)$-th and thus $\sigma_{r}/\langle r\rangle\to\infty$. On the other hand, for the GUE this ratio is given by a finite number of order unity, $\sigma_{r}^{2}/\langle r\rangle^{2}=256\pi^{2}/(27\sqrt{3}-4\pi)^{2}-1\simeq 1.160$. Figure 4 $(e)$ shows that the Poissonian statistics are only attained in the dissipationless limit $N^{1/2}g_{\mathrm{eff}}\to 0$. On the other hand, the GUE values are attained for $g_{\mathrm{eff}}N^{1/2}\simeq 1$. Thus, in the thermodynamic limit, the effective Hamiltonian $\mathcal{H}$ is quantum chaotic for all finite values of $g_{\mathrm{eff}}$.

5 Conclusions

In summary, by analysing an ensemble of stochastic Liouvillians of Lindblad form, where unitary dynamics coexists with $r$ independent dissipation channels, we find that the dispersion of the decay rates, the spectral gap, and the steady-state properties are divided into weak-dissipation, crossover and strong-dissipation regimes, whose boundaries do not necessarily coincide for the different observables. For a given observable, each regime and its boundaries are characterised by a set of scaling exponents ruling the dependence on the effective dissipation strength, $g_{\mathrm{eff}}$, and system size $N$. We determine these exponents and find them to be independent of the universality index $\beta$ of random matrices. It is important to note that the different scaling regimes refer to the value of the effective coupling constant $g_{\mathrm{eff}}$. Since $g_{\mathrm{eff}}=(2r\beta N)^{1/4}g$, our large-$g_{\mathrm{eff}}$ regime is, in fact, what should be typically observed in the thermodynamic limit, for any constant finite $g$ [46].

As dissipation increases, the support of the spectrum passes from an ellipse to the lemon-like shape reported in [34] for the case of $r=N^{2}-1$. For fixed $g_{\mathrm{eff}}$ and $r>1$, the spectral gap increases with $N$ and its distribution becomes peaked in the $N\to\infty$ limit. The case of a single dissipation channel, $r=1$, is qualitatively different. For $r=1$, $\Delta$ vanishes with increasing $N$ for infinite dissipation. Finally, with increasing system size and $r>1$, the steady-state purity approaches that of the maximally mixed state, $\mathcal{P}_{\mathrm{FM}}=1/N$, in the week dissipative regime, while for strong dissipation it attains a value larger than, though proportional to, $\mathcal{P}_{\mathrm{FM}}$. Interestingly, the steady-state spectral statistics exhibit a crossover from Poissonian to GUE as a function of $g_{\mathrm{eff}}N^{1/2}$.

The generic features identified can be contrasted to those of integrable (regular) systems. The lemon-like shape at large dissipation is generic ([34] discusses the independence from specific sampling schemes), with a universal crossover to an ellipse-shaped spectral support. For integrable systems, the spectral density is highly model-dependent; furthermore, there are, in general, lines of eigenvalues of high degeneracy emanating from the bulk spectrum. The correlations inside the spectrum are also markedly different, according to a generalisation of the quantum chaos conjecture [46]. The steady-state of integrable systems is, furthermore, expected to be Poisson-like for an extensive range of parameters, while in random states it is highly suppressed in the thermodynamic limit. We leave a more detailed comparison for future work.

A natural question our results raise is—which regime characterises a particular physical system? As in the case of chaotic Hamiltonian dynamics, this has to be determined on a case-by-case basis. One must also check whether any non-Markovian effects arise. We delegate these questions to future works.

Finally, an important direction of further research is the extension of the results on random Liouvillians to other symmetry classes. We note that the Altland-Zirnbauer tenfold classification of (closed) Hamiltonian dynamics [47] has to be extended for non-Hermitian Hamiltonians, leading to the 38-fold classification of Bernard and LeClaire [48, 49]. The symmetries of Liouvillian dynamics restrict the allowed symmetric classes for quadratic Liouvillians back to ten, agreeing with the Altland-Zirnbauer classes in the dissipationless limit [50]. A Liouvillian in the new non-standard symmetry class AI† was recently reported in [51]. This opens the door for novel effects in open quantum systems to arise in these symmetry classes, with potentially high impact in condensed matter and optical setups.

Appendices

In the Appendices we provide a summary of the various scaling exponents (for scaling with $N$ and $g_{\mathrm{eff}}$) for $r>1$ (given in the main text for $r=2$) and for $r=1$ (A). We also give further details for the evolution of global spectral properties (B), the spectral gap (C) and steady-state properties (D) for $r>1$ and compare with the results for $r=1$ (E).

Appendix A Summary of exponents

The compatibility-of-scaling-function arguments given in the main text for $X$ and $\sigma_{\rho_{0}}^{2}$ can be given in general, for any quantity which has the same qualitative behaviour as the ones described in this work. This procedure also allows us to systematize the exponents found in the main text, by defining three types of exponents: $\nu$, $\kappa$, and $\lambda$, related to the finite-size scaling of the spectral and steady-state quantities (say, $X$), the finite-size scaling of the boundaries of the multiple regimes (say, P, C, D) and the $g_{\mathrm{eff}}$-scaling of the quantities in the crossover regime (say, of $X$ in C), respectively.

We define three regimes with $N$-dependent boundaries in which some quantity $Q$ has qualitatively distinct behaviours: PQ, for $g_{\mathrm{eff}}\lesssim(\beta N)^{\kappa_{Q}^{<}}$; CQ for $(\beta N)^{\kappa_{Q}^{<}}\lesssim g_{\mathrm{eff}}\lesssim(\beta N)^{\kappa_{Q}^{>}}$; and DQ for $(\beta N)^{\kappa_{Q}^{>}}\lesssim g_{\mathrm{eff}}$. In PQ, $Q$ behaves as

[TABLE]

while in DQ,

[TABLE]

Note that some extra factors of $g_{\mathrm{eff}}$ may exist (such as the extra $g_{\mathrm{eff}}^{2}$ for $X$) but they do not modify the argument, as long as they are the same in regimes PQ and DQ. If they are different, a straightforward modification of (21) below is required, but this issue did not arise for the quantities studied in this work. We further assume that the scaling functions $f_{Q}$ have asymptotic power-law behaviours, that is,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The two limiting behaviours of $Q$ should match in the intermediate regime CQ, that is,

[TABLE]

whence the equality

[TABLE]

follows. This equality implies that $\lambda_{Q}^{<}=\lambda_{Q}^{>}\equiv\lambda_{Q}$ and establishes a relation between $\nu$, $\kappa$ and $\lambda$ exponents, which are thus not all independent but constrained by

[TABLE]

Table 1 shows the values of the various exponents for $r=1$ and $r>1$, as defined above. It can be checked that all satisfy (21). Note that the exponent $\nu_{Q}^{C}$ (giving the collapse of curves of different $\beta N$ in regime $C_{Q}$) is not defined in the above argument and thus does not enter the constraint.

Appendix B Global Spectrum

The scaling of spectral quantities can be motivated as follows. The spectral density of $H$ is given by the Wigner semicircle law, with zero mean and standard deviation $\sqrt{\beta N/2}$. The dissipation matrix $d$ is a sum of $r$ Wishart matrices, hence its spectral density follows a Marchenko-Pastur law, with mean $\beta Nrg^{2}$ and standard deviation $\sqrt{r}\beta Ng^{2}$. Since the eigenvalues of $H$ and $d$ are real, when considered separately, i.e. at the non-dissipative ($g=0$) or at the fully dissipative ($g=\infty$) limits, respectively, they only contribute to the imaginary or to the real parts of $\Lambda$. Although at finite $g$ these considerations are no longer exact, we expect them to yield the leading scaling behaviour. Thus, only the dissipative term contributes to the mean of $\Lambda$ and we have $R\propto\beta Nrg^{2}$, see Figure 5. At large $g$, the main contribution to $X$ comes from the dissipative term and $X\propto\sqrt{r}\beta Ng^{2}$ (see Figure 1 $(d)$ in the main text). For weak dissipation the $Y$ scaling is dominated by the Hamiltonian term, hence $Y\propto\sqrt{\beta N}$ (see Figure 1 $(e)$ in the main text). The passage from the weak dissipation scalings to the strong dissipation scalings should occur when the two terms in the Liouvillian are of the same order, $\sqrt{\beta N/2}\simeq\beta N\sqrt{r}g^{2}$, i.e. $g_{\mathrm{eff}}\simeq 1$.

By fitting the global spectral data to power-laws, we find very good agreement with the exponents in Table 1. $\nu_{C}=\nu_{D}$ is completely independent of $r$. $\nu_{P}$ seems to decrease slightly with $r$, although the dependence is weak, and is most likely due to the reduced sizes possible to attain numerically. Furthermore, the exponents $\lambda_{X}=2$ and $\kappa_{X}^{<}=-1/4$ are compatible with data extrapolation for $N\to\infty$, see Figure 6. For both large $g_{\mathrm{eff}}$ (Figure 6 $(a)$) and small $g_{\mathrm{eff}}$ (Figure 6 $(b)$), the extrapolated data to $N\to\infty$ (black crosses), taken as the $y$-intercept of a linear fit in $(\beta N)^{-1}$, agrees very well with the power-law $g_{\mathrm{eff}}^{\lambda_{X}}$ (gray line) in regime C. Upon entering regimes P and D, the extrapolation points naturally start deviating from the power-law curve.

Appendix C Spectral Gap

We now consider the spectral gap for $r>1$. For all regimes, the variance of the gap decreases with $N$ and the value of $\Delta$ becomes sharply defined around its mean, see the gap distribution functions depicted in Figure 7. Hence, the mean gap accurately describes the long-time dynamics.

The mean of the different distributions in Figure 7 is not constant due to the finite-size effects, which we address next. Fitting directly the numerical average spectral gap for several $N$ and multiple $r>1$ seems to indicate that there is a dependence of the scaling exponents $\nu_{P^{\prime}}$, $\nu_{C^{\prime}}$, $\nu_{D^{\prime}}$ with $r$. In particular, $\nu_{P^{\prime}}=\nu_{C^{\prime}}=\nu_{D^{\prime}}$ only for large $r$ but $\nu_{P^{\prime}}>\nu_{C^{\prime}}>\nu_{D^{\prime}}$ for small $r>1$. However, this is incompatible with both the picture drawn in A and with the finiteness of the gap in the thermodynamic limit for $r>1$, discussed below. Since we can only access relatively small $N$, this apparent $r$-dependence is likely a finite-size effect suppressed by $r$. Indeed, values of the exponents obtained by data extrapolation to large $N$ are compatible with $\nu_{P^{\prime}}=\nu_{C^{\prime}}=\nu_{D^{\prime}}=1/2$, independently of $r$.

We next provide the analysis yielding the results presented in the main text on the limiting behaviour of the spectral gap at small and large $g_{\mathrm{eff}}$. For strong dissipation, we have $\mathcal{L}\simeq\sum_{\ell}\mathcal{D}_{W_{\ell}}$. The endpoints of the spectrum of $\mathcal{L}$ along the real direction in the complex plane (that is, the spectral gap and left-most point of the lemon-shaped support, which is the gap plus some order-unity multiple of $X$) can be accessed directly by the use of the holomorphic Green’s function of $\mathcal{L}$ [52, 53], $G=\mathrm{Tr}\langle(z-\mathcal{L})^{-1}\rangle$. This was exploited in [36] to compute the spectral form factor. Unfortunately, the information along the imaginary axis (for instance $Y$) and the spectral density inside the support cannot be accessed by this method.

When computing the holomorphic Green’s function, it can be shown [36, 40] that the terms $W_{\ell}\otimes W^{*}_{\ell}$ do not contribute in the large-$N$ limit. For a single jump operator, the superoperator under consideration is then of the form $\mathcal{L}\simeq-\frac{1}{2}(\Gamma\otimes\mathbbm{1}+\mathbbm{1}\otimes\Gamma$), where $\Gamma=W^{\dagger}W$ is a square Wishart matrix. Since the two sectors of the Liouvillian tensor product representation do not see each other, we can do independent expansions for each Wishart matrix. Then, $n^{\mathrm{th}}$-order terms of the expansion of the Liouvillian Green’s function are given by all possible combinations of $n_{1}^{\mathrm{th}}$-order terms of one of the Wishart matrices with all $n_{2}^{\mathrm{th}}$-order terms of the other Wishart matrix, restricted by $n_{1}+n_{2}=n$ and there are ${n}\choose{n_{1}}$ such combinations for each fixed $n$. By summing all terms, the holomorphic Green’s function of $\mathcal{L}$ is then

[TABLE]

where $C_{k}$ counts the number of Wishart planar diagrams at $k^{\mathrm{th}}$ order and is given by the Catalan numbers, whose real-integral representation is $C_{n}=(1/2\pi)\int_{0}^{4}\mathop{}\!\mathrm{d}x\,x^{n-1}\sqrt{4x-x^{2}}$, i.e. the moments of the Marchenko-Pastur distribution, $\varrho_{\mathrm{MP}}(x)=\sqrt{x(4-x)}/(2\pi x)$, $0<x<4$. Inserting the integral representation of the Catalan numbers into (22), performing the sum and comparing with the relation between the holomorphic Green’s function and the spectral density $\varrho_{\mathcal{L}}$, $G_{\mathcal{L}}(z)=\int\mathop{}\!\mathrm{d}\nu\,\varrho_{\mathcal{L}}(\nu)/(z-\nu)$, we arrive at

[TABLE]

with $\nu_{m}=\max\{0,-2x-4\}$ and $\nu_{M}=\min\{4,-2x\}$, whence it follows that $-4<x<0$. Hence, the computation of the holomorphic Green’s function has lead us to a convolution of two Marchenko-Pastur laws, which is the spectral density of half the sum of two uncorrelated Wishart matrices. Although the spectral density $\varrho_{\mathcal{L}}$ does not correctly describe the spectrum of the Liouvillian, it does predict its endpoints and the interval $[-4,0]$ in which $\mathrm{Re}(\Lambda)$ is supported. Note, in particular, that the spectrum is gapless for a single decay channel and infinite dissipation.

For more than one jump operator, $\Gamma$ is a sum of $r$ identical (i.e. independent but sampled from the same distribution) Wishart matrices. In [36] it was shown, using free-addition of random matrices, that the spectral density of $\Gamma$ is again a Marchenko-Pastur distribution, but with endpoints dependent on $r$,

[TABLE]

where $\xi_{\pm}=(1\pm\sqrt{r})^{2}$. But we have just seen that the computation of the spectral density from the holomorphic Green’s function is equivalent to computing it for the half-sum of two uncorrelated matrices $\Gamma$, whence it immediately follows that $\mathcal{L}$ is supported in $[-\xi_{+},-\xi_{-}]$. We thus see, in agreement with [36], that

[TABLE]

while

[TABLE]

with the proportionality constant of order unity, and where we have reintroduced the variances $\beta Ng^{2}$. The latter result was already verified in the main text and motivated in an alternative way in the previous section.

The prediction for the gap at large $g_{\mathrm{eff}}$ for different $r$ is depicted in Figure 8, together with extrapolated data to $N\to\infty$ for several values of $r$, and in Figure 3 (main text) for $r$ up to $60$. We see that even at small $r$, the theoretical prediction correctly describes the spectral gap for a wide range of $g_{\mathrm{eff}}$ (that is, in the whole D*′* regime). Note also that $\langle\Delta\rangle\to 0$ at large $g_{\mathrm{eff}}$ for $r=1$ in the thermodynamic limit, consistent with our results in E.

For small $g_{\mathrm{eff}}$, we can proceed perturbatively. If $H$ is diagonal in the basis $\left|k\right\rangle$, $H=\sum_{k}\varepsilon_{k}\left|k\right\rangle\left\langle k\right|$, then the eigenstates of $\mathcal{L}$ at $g=0$ are $\rho_{kk^{\prime}}=\left|k\right\rangle\left\langle k^{\prime}\right|$ with eigenvalue $\Lambda_{kk^{\prime}}=-i(\varepsilon_{k}-\varepsilon_{k^{\prime}})$ ($k,k^{\prime}=1,\dots,N$). We thus have a $N$-fold degeneracy of stationary eigenstates, while the other eigenvalues form $N(N-1)/2$ complex conjugate pairs. Any amount of dissipation lifts the degeneracy of the zero-eigenvalue subspace but for small enough dissipation the two sectors will not mix and we expect that the long-time dynamics (the gap) and the (now unique) steady state are determined from the states in the (previously-) degenerate subspace. We lift the degeneracy by diagonalising the perturbation $\sum_{\ell}\mathcal{D}_{W_{\ell}}$ in that subspace. We obtain the eigenvalues of $\mathcal{L}$ to first order,

[TABLE]

with the superoperator notation $\left|\left|mn\right\rangle\right\rangle=\left|m\right\rangle\otimes\left\langle n\right|^{\sf T}$.

From (27) it follows that

[TABLE]

[TABLE]

[TABLE]

which are the conditions for $A$ to be the generator of a classical stochastic equation [41, 34], $\partial_{t}\mathbf{P}=A\mathbf{P}$, where $\mathbf{P}$ is a probability vector. The diagonal entries $A_{mm}=-\sum_{k\neq m}A_{km}$ are determined by the off-diagonal entries, hence we must only determine the distribution of the latter. For $n\neq m$, the distribution of the entry $A_{nm}$ depends on the nature of the jump operator as follows. $W^{(\ell)}_{nm}$ is a Gaussian random variable with variance $g^{2}$ with $\beta$ real degrees of freedom. Then $A_{nm}=\sum_{\ell}{W^{(\ell)}_{nm}}^{*}W^{(\ell)}_{nm}=\sum_{\ell}\left|W^{(\ell)}_{nm}\right|^{2}$ is (in terms of real degrees of freedom) a sum of $r\beta$ squared Gaussian iid variables and hence follows a $\chi^{2}$-distribution with $k=r\beta$ degrees of freedom. The distribution function of the entries $A_{nm}$ is

[TABLE]

Ensembles of real exponentially-distributed matrices (corresponding to $k=2$) with the diagonal constraint (30) were studied in [41]. In particular, it was shown that, in the large-$N$ limit the spectral density depends only on the second moment of the distribution, $\tau^{2}_{k}\equiv\langle A_{nm}^{2}\rangle-\langle A_{nm}\rangle^{2}=2kg^{4}$, and can therefore be calculated from any distribution of entries having this variance, particularly for a Gaussian distribution, which was done in [54] for symmetric $A$. Furthermore, [41] numerically found that the average smallest nonzero eigenvalue of $A$ (the symmetric value of the spectral gap here) is $N\langle R\rangle$ for large $N$, where $\langle R\rangle$ is the mean of the exponential ($k=2$) distribution of $A_{nm}$. Generalising to our $\chi^{2}$-distribution, with mean $\mu_{k}=kg^{2}$, we find that $\langle\Delta\rangle=\beta Nrg^{2}=\sqrt{\beta Nr/2}g_{\mathrm{eff}}^{2}$.

The comparison of this result with numerical data for the average spectral gap for various $r$ is given in Figure 8 as a function of $g_{\mathrm{eff}}$ and for extrapolated $N\to\infty$ and in the inset of Figure 3 (main text) for $r$ up to $r=60$ for fixed $N=60$ and two fixed $g$ (one for strong dissipation, one for weak). We see that although there is some deviation from the theoretical curves at small $r$, they describe the numerical data increasingly well with growing $r$.

Finally, we note that the exponent $\lambda_{\Delta}$ is not defined. The arguments of A do not provide the value of $\lambda_{\Delta}$ since $\nu_{P}^{\prime}=\nu_{D}^{\prime}$ and $\kappa_{\Delta}^{<}=\kappa_{\Delta}^{>}$. This stems from the fact that, for the gap, different scaling functions for small and large dissipation, which asymptotically match in an intermediate regime of size growing with $N$, do not exist. Instead, the gap is described by a single scaling function for all $g$.

Appendix D Steady state

Next, we analyse the steady-state properties for $r>1$. The finite-size scaling exponents $\nu_{P^{\prime\prime}}=-3$ and $\nu_{C^{\prime\prime}}=\nu_{D^{\prime\prime}}=-2$ are in very good agreement with fits of the data to the respective power-laws, for various $r>1$ ($r=2$, $3$, $5$, $10$) and show no dependence on $r$ as expected. The exponents $\lambda_{\rho_{0}}$ and $\kappa_{\rho_{0}}^{<}$ are also compatible with data extrapolation, see Figure 9, by proceeding in the same way as for $X$.

Appendix E Spectral and steady-state properties for $r=1$

For weak dissipation, the unitary contribution to the Liouvillian dominates and, therefore, the spectral and steady-state properties do not differ qualitatively between $r=1$ and $r>1$. In fact, in regime P*′′*, the exponents are always the same in both cases, see Table 1. However, there are important differences at large $g_{\mathrm{eff}}$, where dissipation dominates, for both the spectral gap and the steady state.

Contrary to $r>1$, for $r=1$ the scaling of the gap at large $g_{\mathrm{eff}}$ and at small $g_{\mathrm{eff}}$ is not the same. In particular, since $\nu_{D^{\prime}}$, $\lambda_{\Delta}<0$, in the thermodynamic limit the gap starts to close as dissipation increases (strictly closing only for $g_{\mathrm{eff}}=+\infty$). This can be seen in Figure 10, where we plot the evolution of $\langle\Delta\rangle$ as a function of $g_{\mathrm{eff}}$: choosing the small-$g_{\mathrm{eff}}$ scaling (which coincides with the scaling for all $g_{\mathrm{eff}}$ for $r>1$), at large $g_{\mathrm{eff}}$ the gap follows the dashed line (power-law $g_{\mathrm{eff}}^{2-\lambda_{\Delta}}$) and hence closes at $g_{\mathrm{eff}}=\infty$. Once again, the exponents listed in Table 1 are compatible with data extrapolation, see the black crosses in Figure 10. The procedure for data extrapolation is the same as before.

Furthermore, for the $r=1$ case, contrary to $r>1$, the values of $g_{\mathrm{eff}}$ that define the boundaries between different regimes are $N$-dependent. Hence, $\lambda_{\Delta}$ can be defined for $r=1$. The gap has a nonzero exponent $\kappa_{\Delta}^{>}$ and the matching procedure of A is applicable. In particular, with $\kappa_{\Delta}^{>}=3/4$ (see Figure 10 $(b)$), we can use (21) to determine $\lambda_{\Delta}=-8/3$. Hence, we find, at large dissipation, $\langle\Delta\rangle\propto g_{\mathrm{eff}}^{-2/3}$ in agreement with [36].

Regarding the steady state at large dissipation and $r=1$, in Figure 11 we show the evolution of the variance of steady-state eigenvalues as a function of $g_{\mathrm{eff}}$. For $r=1$ the condition $\nu_{C^{\prime\prime}}=\nu_{D^{\prime\prime}}$ is no longer verified, instead $\nu_{C^{\prime\prime}}=-2$. Using the values of $\nu_{P^{\prime\prime}}$ and $\kappa_{\rho_{0}}^{<}$ obtained from Figure 11 $(a)$, and of $\nu_{D^{\prime\prime}}$ and $\kappa_{\rho_{0}}^{>}$ obtained from Figure 11 $(c)$, (21) gives the value $\lambda_{\rho_{0}}=8/5$. Note that in this case, both boundaries of regime C*′′* scale with $N$ and the scaling function does not have a single power-law behaviour throughout the intermediate regime. However, the power-law $g_{\mathrm{eff}}^{\lambda_{\rho_{0}}}$ describes accurately the behaviour of $\sigma_{\rho_{0}}^{2}$ both when entering regime C*′′* from regime P*′′* and when entering regime C*′′* from regime D*′′, see the dashed lines in Figure 11 $(b)$. The two power-laws are shifted with respect to each other and are linked by a crossover regime in the center of regime C′′*.

Finally, we have $\sigma^{2}_{\rho_{0}}\propto 1/N$ for strong dissipation, which gives a purity $\mathcal{P}_{0}\propto N^{0}$ in the large-$N$ limit. Thus, contrarily to the $r>1$ case, the steady-state purity is not proportional to that of the maximally mixed state but remains finite for large $N$.

We gratefully acknowledge valuable discussions with Tankut Can. LS acknowledges support by FCT through PhD Scholarship SFRH/BD/147477/2019. PR acknowledges support by FCT through the Investigador FCT contract IF/00347/2014 and Grant No. UID/CTM/04540/2019. TP acknowledges ERC Advanced Grant 694544-OMNES and ARRS research program P1-0402.