João Costa

Pedro Ribeiro

Andrea De Luca

Tomaž Prosen

Lucas Sá

Spectral and steady-state properties of fermionic random quadratic Liouvillians

come

Introduction

The vast majority of systems in nature have their own time evolution deeply influenced by the interaction with their environment. Under the assumption of either a very weakly or very strongly coupled environment, with memory times much shorter than all other characteristic time scales (Markovian approximation), the time evolution equation for the system's reduced density matrix ρ, assumes the Lindblad form [START_REF] Belavin | Relaxation of quantum systems with equidistant spectra[END_REF][START_REF] Gorini | Completely positive dynamical semigroups of N-level systems[END_REF][START_REF] Lindblad | On the generators of quantum dynamical semigroups[END_REF]:

dρ d t = L [ρ] = -i Ĥ, ρ + g M µ=1 2 Lµ ρ L † µ -{ L † µ Lµ , ρ} , (1)
where the superoperator L is known as the Liouvillian and g is a parameter that quantifies the dissipation strength. The M independent jump operators, Lµ , represent channels of interaction with the environment that act, e.g., as sources of dephasing and dissipation. In the absence of these operators, the evolution is that of a closed system, and the Liouvillian becomes just the von Neumann generator, -i Ĥ, ρ . Note that, although for strictly zero dissipation the ensuing unitary time evolution is completely determined by the Hamiltonian, Ĥ, for any finite dissipation strength there is a (generically unique) steady-state the system relaxes to at large times. Although the Markovian approximation leads to a considerable simplification to the time-evolution equation, it leaves out the possibility of studying, for instance, some quantum transport setups,for which different approaches need to be taken [START_REF] Ribeiro | Non-Markovian effects in electronic and spin transport[END_REF][START_REF] Ribeiro | Steady-state properties of a nonequilibrium Fermi gas[END_REF]. Nevertheless, it still finds many valuable applications in various subjects, namely in quantum optics and quantum computation [START_REF] Gardiner | Quantum noise: a handbook of Markovian and non-Markovian quantum stochastic methods with applications to quantum optics[END_REF].

Despite the clear simplification the Markovian approximation brings to the problem, the determination of the spectral and steady-state properties of the Liouvillian for generic Hamiltonian and jump operators remains a formidable task and is the object of intense ongoing research. A further simplification can be achieved if we restrict our analysis to quadratic systems, which are characterized by a quadratic Hamiltonian and linear jump operators in bosonic or fermionic creation and annihilation operators [START_REF] Prosen | Third quantization: a general method to solve master equations for quadratic open Fermi systems[END_REF][START_REF] Prosen | Quantum Phase Transition in a Far-from-Equilibrium Steady State of an X Y Spin Chain[END_REF][START_REF] Prosen | Spectral theorem for the Lindblad equation for quadratic open fermionic systems[END_REF][START_REF] Prosen | Quantization over boson operator spaces[END_REF][START_REF] Dzhioev | Super-fermion representation of quantum kinetic equations for the electron transport problem[END_REF][START_REF] Guo | Solutions for bosonic and fermionic dissipative quadratic open systems[END_REF][START_REF] Lieu | Tenfold way for quadratic Lindbladians[END_REF][START_REF] Barthel | Solving quasi-free and quadratic lindblad master equations for open fermionic and bosonic systems[END_REF][START_REF] Kawasaki | Topological phases protected by shifted sublattice symmetry in dissipative quantum systems[END_REF]. More precisely, and focusing on a system with N F complex fermions satisfying {c i , c † j } = δ i j , the Liouvillian of Eq. ( 1) is said to be quadratic if

Ĥ = 1 2 2N F i, j=1 C † i H i j C j and Lµ = 2N F j=1 l µ j C j , (2)
with C = {c 1 , . . . , c N F , c † 1 , . . . , c † N F } T a vector of fermionic creation and annihilation operators that satisfies {C i , C † j } = δ i j . The quadratic Lindblad operator obtained from this construction ensures that the dynamics preserve the Gaussian form of an initial density matrix. Thus, the time evolution of the 2 N F × 2 N F density matrix can be encoded by its second moments matrix-the correlation matrix-of size 2N F × 2N F . Analogously to quadratic Hamiltonian systems, it is possible to construct a single-particle basis whose dimension scales linearly with the number of fermionic modes, N F . Many-body observables, such as the Liouvillian's many-body spectrum and steadystate correlators, can be straightforwardly computed from single-particle quantities. Moreover, the single-body spectrum can be identified with that of a non-Hermitian Hamiltonian [START_REF] Prosen | Third quantization: a general method to solve master equations for quadratic open Fermi systems[END_REF], leaving the determination of the Liouvillian's spectral properties only dependent on the specification of the single-particle Hamiltonian H and jump operators l.

However, most systems of interest are extremely complex, with many degrees of freedom and exhibiting very complicated dynamics, rendering the task of determining these operators impossible in practice. We thus resort to Jayne's principle of maximal entropy [START_REF] Jaynes | Information Theory and Statistical Mechanics[END_REF][START_REF] Jaynes | Information Theory and Statistical Mechanics[END_REF], constraining these operators to a manifold compatible with their symmetries and randomizing all other degrees of freedom. This principle has proven immensely successful for complex closed quantum systems, as pioneered by Wigner, who proposed to approximate the Hamiltonian of a compound nucleus by a random matrix [START_REF] Wigner | On the statistical distribution of the widths and spacings of nuclear resonance levels[END_REF]. The success of the approach led to the formulation of the celebrated Bohigas-Giannoni-Schmit conjecture [START_REF] Bohigas | Characterization of Chaotic Quantum Spectra and Universality of Level Fluctuation Laws[END_REF] that states that the spectral correlations of quantum chaotic Hamiltonians coincide with those of a random matrix of the appropriate symmetry class. Besides level correlations, also level densities are captured by many-body random matrix models with few-body interactions (quartic in creation and annihilation operators), the random embedded ensembles [START_REF] French | Validity of random matrix theories for many-particle systems[END_REF][START_REF] French | Some random-matrix level and spacing distributions for fixedparticle-rank interactions[END_REF][START_REF] Bohigas | Two-body random Hamiltonian and level density[END_REF][START_REF] Bohigas | Spacing and individual eigenvalue distributions of two-body random Hamiltonians[END_REF][START_REF] Mon | Statistical properties of many-particle spectra[END_REF][START_REF] Brody | Randommatrix physics: spectrum and strength fluctuations[END_REF] and the Sachdev-Ye-Kitaev (SYK) model [START_REF] Sachdev | Gapless spin-fluid ground state in a random quantum Heisenberg magnet[END_REF][START_REF] Kitaev | Hidden correlations in the Hawking radiation and thermal noise[END_REF][START_REF] Kitaev | KITP Program: Entanglement in Strongly-Correlated Quantum Matter[END_REF][START_REF] Kitaev | KITP Program: Entanglement in Strongly-Correlated Quantum Matter[END_REF][START_REF] Sachdev | Bekenstein-Hawking entropy and strange metals[END_REF]. Finally, the interplay of single-body chaos and many-body integrability in random quadratic fermionic Hamiltonians has also been studied [START_REF] Cotler | Black holes and random matrices[END_REF][START_REF] Magán | Random Free Fermions: An Analytical Example of Eigenstate Thermalization[END_REF][START_REF] Łydżba | Eigenstate Entanglement Entropy in Random Quadratic Hamiltonians[END_REF][START_REF] Liao | Many-Body Level Statistics of Single-Particle Quantum Chaos[END_REF][START_REF] Winer | Exponential Ramp in the Quadratic Sachdev-Ye-Kitaev Model[END_REF].

More recently, the random matrix theory (RMT) approach has been extended to generic open quantum systems with Markovian dissipation [START_REF] Denisov | Universal Spectra of Random Lindblad Operators[END_REF][START_REF] Can | Random Lindblad dynamics[END_REF][START_REF] Can | Spectral Gaps and Midgap States in Random Quantum Master Equations[END_REF][START_REF] Sá | Spectral and steady-state properties of random Liouvillians[END_REF][START_REF] Sá | Spectral transitions and universal steady states in random Kraus maps and circuits[END_REF][START_REF] Lange | Random-matrix theory for the Lindblad master equation[END_REF][START_REF] Tarnowski | Random generators of Markovian evolution: A quantum-classical transition by superdecoherence[END_REF]. A random Liouvillian was shown to have a lemon-shaped spectral support [START_REF] Denisov | Universal Spectra of Random Lindblad Operators[END_REF][START_REF] Tarnowski | Random generators of Markovian evolution: A quantum-classical transition by superdecoherence[END_REF] and its spectral gap was extensively studied [START_REF] Can | Random Lindblad dynamics[END_REF]. The dependence of spectral and steady-state properties of the random Liouvillian with system size, dissipation strength, and number of jump operators was addressed in Ref. [START_REF] Sá | Spectral and steady-state properties of random Liouvillians[END_REF]. Considering jump operators with few-body interactions has clarified the role of locality in the separation of dissipative timescales [START_REF] Wang | Hierarchy of Relaxation Timescales in Local Random Liouvillians[END_REF][START_REF] Sommer | Many-body hierarchy of dissipative timescales in a quantum computer[END_REF] and metastability [START_REF] Li | Random matrix theory for quantum and classical metastability in local Liouvillians[END_REF] and allowed for the analytic computation of the spectral gap of the strongly-coupled SYK Liouvillian [START_REF] Sá | Lindbladian dissipation of strongly-correlated quantum matter[END_REF][START_REF] Kulkarni | Lindbladian dynamics of the Sachdev-Ye-Kitaev model[END_REF][START_REF] García-García | Keldysh Wormholes and Anomalous Relaxation in the Dissipative Sachdev-Ye-Kitaev Model[END_REF][START_REF] Kawabata | Dynamical quantum phase transitions in SYK Lindbladians[END_REF]. An important open question is to establish the universality of these results. Encouraging first steps showed that, besides local level statistics [START_REF] Sá | Complex Spacing Ratios: A Signature of Dissipative Quantum Chaos[END_REF], the steady-state of fully random Kraus maps and Liouvillians coincide [START_REF] Sá | Spectral transitions and universal steady states in random Kraus maps and circuits[END_REF] and that global spectral features of quartic Liouvillians are qualitatively similar to the fully random case [START_REF] Sá | Lindbladian dissipation of strongly-correlated quantum matter[END_REF][START_REF] Kulkarni | Lindbladian dynamics of the Sachdev-Ye-Kitaev model[END_REF].

In this paper, we extend this ongoing effort and study the single-body spectral and steadystate properties of fermionic random quadratic Liouvillians. The rest of the paper is organized as follows. In Sec. 1.1 we summarize the main results of this work, which are then worked out in detail in the following sections. In Sec. 2 we review the formalism of quadratic Liouvillians and explain our random sampling. Spectral properties (spectral boundary, phase transition, and spectral gap) are discussed in Sec. 3 and steady-state properties (distribution, purity, and statistics) in Sec. 4. In Sec. 5, we present concluding remarks and possible further directions. Technical calculations and proofs are presented in a series of six appendices.

Main results

In the thermodynamic limit N F → ∞, the single-particle properties of random quadratic Liouvillian, specified by Eqs. ( 1) and [START_REF] Gorini | Completely positive dynamical semigroups of N-level systems[END_REF], are determined by two parameters: the dissipation strength g and the ratio of the number of jump operators to the number of fermions m = M /(2N F ). In Fig. 1(a), we plot the phase diagram of this system in the 1/g versus m plane. For large enough m and small enough g, the system is in phase I, characterized by a single-body spectrum supported on a simply-connected region of the complex plane, see Fig. 1(b). When g is increased or m decreased across a critical value, a phase transition occurs and, in phase II, the single-body spectrum splits into two disconnected components, see Fig. 1(c). The existence of these two regions signals the existence of an intermediate period of metastability, during which an extensive number of modes coexist without (considerable) decay. The critical line separating the two phases [dashed line in Fig. 1(a)] and the boundaries of the single-body spectral support can be computed analytically [START_REF] Haake | Statistics of complex levels of random matrices for decaying systems[END_REF][START_REF] Lehmann | Chaotic scattering: the supersymmetry method for large number of channels[END_REF][START_REF] Janik | Non-hermitian random matrix models[END_REF][START_REF] Janik | Non-Hermitian random matrix models: Free random variable approach[END_REF].

For very strong dissipation g → ∞, the transition occurs at m = 1/2 and corresponds to the decoupling of some fermionic degrees of freedom from the dynamics. Indeed, the Hamiltonian contribution vanishes when g → ∞ and there is an insufficient number of jump operators (M < 2N F ) to couple all 2N F fermionic creation and annihilation operators to the environment. Below the transition [m < 1/2, red line in Fig. 1(a)] the decoupled fermions exhibit nonergodic features, discussed in detail below. As g is lowered to a finite value, the Hamiltonian starts to couple the dissipatively-decoupled fermions and the critical value of m decreases. At weak dissipation g < 1/ 2, the Hamiltonian contribution is strong enough to couple all fermions and there is no transition. The spectral gap, which sets the (inverse) timescale of relaxation to the steady-state, coincides for the single-and many-body spectra. It can also be obtained analytically and assumes a very simple form in the limits g → 0 and g → ∞. For weak dissipation, the spectral gap closes linearly with g and m (for all m), as expected from perturbation theory. On the other hand, for large dissipation, the spectral gap acts as an order parameter of the transition at m = 1/2. For m < 1/2, in the limit g → ∞ the gap closes like 1/g leading to a gapless Liouvillian, signalling a slow approach to the steady-state. For m > 1/2, the gap has the linear scaling in g, typical of a dissipation-driven relaxation, growing with m as ( 2m-1) 2 , as dictated by the Marchenko-Pastur law. At the critical point m = 1/2, the gap closes as (m/g) 1/3 . For large but finite g, the gap is nonzero for any value of m, but still exhibits qualitatively different behaviours above and below the transition.

The steady-state, to which the system relaxes in the long-time limit, is Gaussian for quadratic Liouvillians and is thus fully characterized by its single-particle properties. In the limit g → 0 (where there is no decoupling transition), the steady-state single-body spectrum is Gaussian, fullymixed, and displays Poisson statistics irrespectively of the value of m. As g increases, there is a perturbative crossover well-below the threshold g = 1/ 2 for the appearance of the decoupling transition and the steady state becomes distributed according to RMT [see Fig. 1(d)], is mixed but not fully mixed, and exhibits GUE statistics [see Fig. 1(e)]. The purity interpolates monotonously between the two limits g → 0 and g → ∞.

When g is further increased, the steady state is also influenced by the decoupling transition. In the limit g → ∞, the results can be obtained analytically through perturbation theory. Above the transition, m > 1/2, the properties attained after the perturbative crossover do not change. However, for m < 1/2, the spectrum of the steady state also splits into two independent sectors, see Fig. 1(f). The sector of the fermions coupled to the environment retains the properties of m > 1/2, except for the spectral statistics which, remarkably, crossover from GUE to GOE as m → 0, see Fig. 1(g). The spectrum of the decoupled sector, on the other hand, is composed of uncorrelated Gaussian random variables displaying Poisson statistics (nonergodic behavior), see Fig. 1(g). These two sectors are well-separated for small-enough m, see Fig. 1(f), but overlap for larger m.

We emphasize that the nonergodic features of the steady state were obtained in the limit g → ∞. However, they leave strong imprints in the dynamics at large but finite g and N F . Indeed, the interplay of the two sectors leads to a decrease of the purity, with a nonmonotonic behavior as a function of g, and to an interpolation between RMT and Poisson statistics.

Finally, we note that the properties of random quadratic Liouvillians above the transition are quantitatively similar to those identified in previous studies of fully random Liouvillians with unconstrained interactions [START_REF] Sá | Spectral and steady-state properties of random Liouvillians[END_REF][START_REF] Sá | Spectral transitions and universal steady states in random Kraus maps and circuits[END_REF] 1 . The nonergodic behaviour below the transition is, however, not accessible to these fully-random Liouvillians.

Random quadratic Liouvillians

Liouvillian dynamics

The dynamics of the system can be entirely specified by looking at the Liouvillian eigenvalue problem. In fact, after determining a complete set of Liouvillian eigenmodes ρ j with eigenvalue Λ j [L(ρ j ) = Λ j ρ j ], we can write the evolution of any initial state ρ(0) = j c j ρ j as2 ρ(t) = j c j e Λ j t ρ j .

(

) 3
All Λ j have a non-positive real part as ensured by the complete positivity of the Lindblad equation. Moreover, the Hermiticity preservation property guarantees that eigenvalues are real or come in complex-conjugated pairs. Generically, because of trace preservation, there is a single eigenstate, the steady state ρ NESS , with corresponding zero eigenvalue. ρ NESS is invariant under time evolution, as

dρ NESS d t = L(ρ NESS ) = 0. ( 4
)
If there are no other eigenvalues with Re Λ j = 0, the system will relax to the steady state as t → ∞, since all the other eigenmodes of the Liouvillian decay to zero. The rate at which the system relaxes to the steady state is dictated by the Liouvillian spectral gap, defined as Gap = -max

Λ j =0
Re(Λ j ).

(5)

Vectorization and adjoint fermions

To study the eigenvalue problem of the quadratic Liouvillian superoperator, it is convenient to recast the it as a matrix acting in an enlarged Hilbert space, a procedure known as vectorization.

In the Fock space of N F fermions, pure states are represented by 2

|α〉 〈β| → |α〉 ⊗ 〈β| T , (6)
Ô1 |α〉 〈β| Ô † 2 → Ô1 ⊗ Ô † T 2 |α〉 ⊗ 〈β| T . ( 7
)
Following this procedure, the Liouvillian [Eq. [START_REF] Belavin | Relaxation of quantum systems with equidistant spectra[END_REF]] is mapped to:

L = -i Ĥ ⊗ -⊗ Ĥ T + g M µ=1 2 Lµ ⊗ L * µ -L † µ Lµ ⊗ -⊗ ( L † µ Lµ ) T . ( 8
)
It is possible to generalize the notion of creation and annihilation operators to this space while keeping the Liouvillian quadratic, resembling the Hamiltonian of a free theory [START_REF] Prosen | Third quantization: a general method to solve master equations for quadratic open Fermi systems[END_REF]. To enforce the canonical anticommutation relations in the vectorized representation also, we define the vector

ã = c 1 ⊗ e iπN T , . . . , c N F ⊗ e iπN T , c † 1 ⊗ e iπN T , . . . , c † N F ⊗ e iπN T , ⊗ c † T 1 e iπN T , . . . , ⊗ c † T N F e iπN T , ⊗ e iπN T c T 1 , . . . , ⊗ e iπN T c T N F T , (9)
where N = i c † i c i is the number operator, which satisfies {ã i , ã † j } = δ i j as required. ãi are known as adjoint fermions. An important feature of this vector is that it is particle-hole symmetric, i.e.,

ã † T = S ã, S = S 0 0 S , S = 0 0 . ( 10
)
S and S implement the particle-hole in the original and vectorized spaces, respectively. This vectorization scheme is equivalent to third quantization [START_REF] Prosen | Third quantization: a general method to solve master equations for quadratic open Fermi systems[END_REF]. For an explicit demonstration, we refer the reader to Appendix A.

Single-particle spectrum and diagonalization

Recall that, as mentioned in Sec. 1 [see Eq. ( 2)], the Liouvillian in Eq. (1) for a system of N F fermions is said to be quadratic if

Ĥ = 1 2 2N F i, j=1 C † i H i j C j and Lµ = 2N F j=1 l µ j C j .
The 2N F × 2N F matrix H, the so-called single-particle Hamiltonian, is Hermitian and can always be chosen to satisfy particle-hole symmetry, SH T S = -H. In turn, the dissipative contribution to the Liouvilian is determined by the M × 2N F non-Hermitian (and in general complex) matrix {l µ j } j=1,...,N F µ=1,...,M (recall that M is the number of independent decay channels).

Next, we define the matrices

N jk = µ l µ j l * µk , N S = SN T S, (11)
the particle-hole-symmetric and antisymmetric combinations

Γ = N + N S 2 , Γ B = N -N S 2 , (12)
and the non-Hermitian single-particle effective Hamiltonian

K = H -i gΓ . (13)
In Appendix B.1 we show that

L = - i 2 ã † L ã - i 2 Tr (K) , (14)
where the single-particle Liouvillian is given by

L = H -i gΓ B -i g N S SJ -i g JSN -J (H + i gΓ B ) T J , (15)
with

J = 0 0 - . ( 16
)
Note that the matrix L also satisfies particle-hole symmetry, S L T S = -L.

As shown in Appendix B.2, we can (almost) always find a change of basis that renders the Liouvillian diagonal,

L = - i 2 2N F i β i b i b i , (17)
where b i and b i (note, however, b i = b † i ) satisfy the canonical anticommutation relations

{b i , b j } = 0, {b i , b j } = 0, and {b i , b j } = δ i j , (18)
and {β j } j∈{1,...,2N F } constitutes the single-body spectrum of the Liouvillian. Moreover, the β i coincide with the eigenvalues of the non-Hermitian Hamiltonian K. It is clear that due to the form of Eq. ( 17), the many-body spectrum of the Liouvillian is completely determined by the single-body spectrum (take all possible sums of subsets of {-iβ j /2} j ), which implies that all its properties are encapsulated in the latter. We will thus focus on studying the features of the single-body spectrum in Sec. 3.1. Since the Liouvillian is quadratic and the single-body spectrum is entirely contained in the left-half plane, the Liouvillian spectral gap corresponds to the gap of the single-body spectrum. Thus, the following definition holds:

Gap = - 1 2 max Imβ i . ( 19
)

Steady state

From the diagonal form [START_REF] Jaynes | Information Theory and Statistical Mechanics[END_REF], the steady state ρ NESS is found to be the state annihilated by all b i . For a quadratic Liouvillian, any initial Gaussian state will remain Gaussian under time evolution, which implies that the steady state must be Gaussian as well. We can thus describe the steady state entirely by its correlation matrix 3 ,

χ i j = Tr ρ NESS C i C † j , (20)
which satisfies the particle-hole symmetry

Sχ T S = -χ. ( 21
)
Alternatively, we can define an effective thermal-like Hamiltonian Ω as

ρ NESS = 1 Z e Ω, Z = Tr e Ω. ( 22
)
This parametrization is convenient as it automatically takes care of the normalization and positivedefiniteness of the steady-state density matrix. Moreover, because the steady state is Gaussian, Ω can be fully characterized by the single-particle matrix Ω,

ρ NESS = 1 Z e 1 2 i j C † i Ω i j C j , (23)
with normalization Z = Tr e

1 2 i j C † i Ω i j C j = det (1 + e -Ω ), (24)
which is related to the correlation matrix through the matrix relation:

χ = 1 + e Ω -1 . ( 25
)
Remarkably, the steady-state correlation matrix can also be entirely determined by the singleparticle non-Hermitian Hamiltonian. Indeed, as shown in Appendix B.3, one can solve the steadystate Lyapunov equation for the correlation matrix:

χ K † -Kχ = i g N . ( 26
)
Alternatively, χ can be constructed more efficiently using the right and left eigenvectors of the single-particle matrix L, see Eq. (99) below and Ref. [START_REF] Prosen | Third quantization: a general method to solve master equations for quadratic open Fermi systems[END_REF]. The equivalence of the two methods is established in Appendix B.3.

Majorana basis

At this point, we change to the Majorana basis, which is more convenient for our purposes. It is implemented by the unitary transformation

U = 1 2 i -i . ( 27
)
From it, we define the Majorana operators

f = { f 1 , . . . , f 2N F } T = 2U C, ( 28
)
where C is the Nambu vector defined after Eq. (2), satisfying f † i = f i and the anticommutation relation

{ f i , f j } = δ i j .
In this basis, K is transformed into a more suitable matrix, K (U) ,

K (U) = H (U) -i gΓ (U) = U KU † ( 29
)
where H (U) is an anti-symmetric matrix,

H (U) T = U * H T U T = USH T SU † = -H (U) , (30)
and Γ (U) is a symmetric matrix,

Γ (U) = U l T l * + Sl † lS 2 U † = 1 2 l (U) T l (U) * + l (U) † l (U) = 1 2 N (U) + N (U) T , ( 31
)
with N (U) = l (U) T l (U) * and l (U) = l U T . (32)

Random sampling

The characterization of the Liouvillian's spectrum and steady state is now completely determined by the specification of matrices H (U) and l (U) , which can only be obtained from the knowledge of the system's Hamiltonian and interactions with the environment. However, assuming that the dynamics is generic, one can argue based on Jayne's principle of maximal entropy that they are well described by a random matrix of a symmetry class consistent with the symmetries of the system. In this case, we restrict H (U) to the set of Hermitian matrices satisfying particle-hole symmetry H (U) T = -H (U) and randomize all other degrees of freedom, corresponding to 2N F ×2N F Gaussian random matrices of class D in the Altland-Zirnbauer classification [START_REF] Altland | Nonstandard symmetry classes in mesoscopic normal-superconducting hybrid structures[END_REF]. The simplest way to achieve this is to draw a matrix H int from the Ginibre orthogonal ensemble (GinOE) [START_REF] Ginibre | Statistical Ensembles of Complex, Quaternion, and Real Matrices[END_REF], i.e., sample a real matrix from the probability distribution ∼ exp{-N F Tr H † int H int } and then set

H (U) = i 2 H int -H T int . ( 33
)
On the other hand, since we do not impose any restriction on the jump operators, we sample l (U) from the Ginibre unitary ensemble (GinUE) [START_REF] Ginibre | Statistical Ensembles of Complex, Quaternion, and Real Matrices[END_REF] of rectangular M × 2N F matrices, i.e., sample a complex matrix from the probability distribution ∼ exp -N F Tr l (U) † l (U) . This is equivalent to saying that the 2N F × 2N F matrix N (U) is drawn from the complex Wishart ensemble [also known as the Laguerre unitary ensemble (LUE)]. Now that we have established a procedure to determine matrices H (U) and l (U) , we can turn to the study of the spectral and steady-state properties of random Liouvillians, which are entirely determined by single-body ones, allowing us to focus just on the properties of the latter. We start with the spectral properties.

Spectral properties

We are interested in the support of the single-body spectrum, as it contains the information on the relevant timescales of the problem. The rightmost boundary point gives the spectral gap (recall that the single-and many-particle gaps coincide). The width of the spectrum along the imaginary axis is related to the timescale for the oscillations of the states' phases. Finally, if the spectral support splits into several components there is a hierarchy of decay times [START_REF] Wang | Hierarchy of Relaxation Timescales in Local Random Liouvillians[END_REF][START_REF] Sommer | Many-body hierarchy of dissipative timescales in a quantum computer[END_REF], with separate sets of modes decaying at different rates, interspersed by periods of metastability [START_REF] Li | Random matrix theory for quantum and classical metastability in local Liouvillians[END_REF].

We first show that our random model can be mapped exactly, in the limit N F → ∞, to a slightly different non-Hermitian Hamiltonian whose boundary has been computed using free probability and use it to identify a phase transition in the single-body spectrum, see Sec. 3.1. Then, in Sec. 3.2, we focus on the spectral gap, studying it it detail, both numerically and analytically, as a function of g and m.

Single-body spectrum

Spectral boundary

As mentioned in Sec. 2 and proven in Appendix B.2, the spectrum of the single-body Liouvillian matrix L coincides with that of K = Hi gΓ . In Refs. [START_REF] Haake | Statistics of complex levels of random matrices for decaying systems[END_REF][START_REF] Lehmann | Chaotic scattering: the supersymmetry method for large number of channels[END_REF][START_REF] Janik | Non-hermitian random matrix models[END_REF][START_REF] Janik | Non-Hermitian random matrix models: Free random variable approach[END_REF], the authors studied the spectrum of a related random matrix H eff = H r -i g r Γ r , with H r and Γ r = AA T being 2N F × 2N F Hermitian matrices drawn from the Gaussian Orthogonal Ensemble (GOE) and the real Wishart ensemble [also known as the Laguerre orthogonal ensemble (LOE)], respectively. (A is a real 2N F × M matrix.) Using replicas [START_REF] Haake | Statistics of complex levels of random matrices for decaying systems[END_REF], supersymmetry [START_REF] Lehmann | Chaotic scattering: the supersymmetry method for large number of channels[END_REF], diagrammatics [START_REF] Janik | Non-hermitian random matrix models[END_REF], or free probability [START_REF] Janik | Non-Hermitian random matrix models: Free random variable approach[END_REF], they established that, in the limit N F → ∞, the spectrum of H eff is supported on a bounded region in the complex plane delimited by a boundary that satisfies the equation

x 2 = - m g r y - g r 1 -2g r y + m 2 y - 1 2g r 2 , ( 34
)
where m = M /(2N F ) and x + i y represents a point in the complex plane 4 .

In our case, the problem is slightly different as, to obtain the single-body spectrum, we need to calculate the eigenvalues of K (U) = H (U) -i gΓ (U) rather than H eff = H r -i g r Γ r . Quite remarkably, apart from numerical prefactors and subleading 1/N F corrections, the eigenvalue distribution of K (U) and H eff coincide, as we argue below.

The difference between H (U) and H r stems from the fact that the former is Hermitian and anti-symmetric whereas the latter is Hermitian and symmetric. Although they belong to distinct symmetry classes, their resolvents and, hence, their eigenvalue distributions match to leading order in 1/N F , as we review in Appendix C. Because the resolvent is the only property of H (U) that enters the determination of the eigenvalue distribution of K (U) , we can interchange it with H r .

In addition, Γ (U) is drawn from a symmetrized complex Wishart ensemble, in contrast to Γ r , which is drawn from the real Wishart ensemble. However, we can also simply draw Γ (U) from the real Wishart ensemble, provided we double the number of jump operators (m

= M /2N F → 2M /2N F = 2m): Γ (U) j,k = 1 2 M µ=1 l (U) µ, j l (U) * µ,k + l (U) µ,k l (U) * µ, j = M µ=1 l (U) R µ, j l (U) R µ,k + l (U) I µ, j l (U) I µ,k = 2M µ=1 l (U) r µ, j l (U) r µ,k , (35)
where

l (U) µ, j = l (U) R µ, j + i l (U) I µ, j and l (U)
r is a 2M × 2N F real matrix built from concatenating l (U) R (taken as the first M rows of l (U) r ) and l (U) I (the last M rows). With the equivalence of the spectra of K (U) and H eff established, to obtain the boundary of the single-body spectrum of the Liouvillian in the limit N F → ∞, we must simply replace m by 2m in Eq. ( 34), together with x by -2 y and y by 2x (since the single-body spectrum is obtained from the spectrum of K by multiplication by -i/2):

y 2 = - m 4g x - g/2 1 -4g x + m 4x - 1 4g 2 . ( 36
)
In Fig. 2, we plot the curve parametrized by Eq. ( 36) in the complex plane and compare it with the single-body spectrum of -iK/2 obtained numerically by exact diagonalization (ED), for three different points in the parameter space (m, g). We observe that it adjusts perfectly to the boundary of the spectrum in all cases. The very small number of outliers can be attributed to finite-size effects, since the boundary becomes sharp in the limit N F → ∞. In particular, one can show [START_REF] Haake | Statistics of complex levels of random matrices for decaying systems[END_REF][START_REF] Lehmann | Chaotic scattering: the supersymmetry method for large number of channels[END_REF][START_REF] Janik | Non-hermitian random matrix models[END_REF][START_REF] Janik | Non-Hermitian random matrix models: Free random variable approach[END_REF] that for N F → ∞, no states lie outside the boundary with probability going to 1. All the spectral information, including the phase diagram in the 1/g versus m plane and the spectral gap, can be extracted from Eq. ( 36), as we discuss in the remainder of this section.

Phase diagram

From Fig. 2, we can observe that two different behaviours of the single-body spectrum emerge for different values of m, for a given g > 1/ 2. For large enough m and small enough g, the singlebody spectrum is supported on a simply-connected region of the complex plane, see Fig. 2(a). We call this region of (m, g) space phase I. When g is increased or m decreased across some critical value g c (m) or m c (g), a phase transition occurs [START_REF] Haake | Statistics of complex levels of random matrices for decaying systems[END_REF][START_REF] Lehmann | Chaotic scattering: the supersymmetry method for large number of channels[END_REF][START_REF] Janik | Non-hermitian random matrix models[END_REF][START_REF] Janik | Non-Hermitian random matrix models: Free random variable approach[END_REF], see Fig. 2(b), and the singlebody spectrum splits into two disconnected components. In phase II, for large enough m and small enough g, the two components of the single-body spectrum are well separated, see Fig. 2(c).

In the 1/g versus m plane, phases I and II are separated by a critical line that can be obtained analytically from Eq. ( 36), see Fig. 1(a). Indeed, assuming the spectrum to be the union of convex sets in the complex plane (an assumption verified in all numerical simulations) and since the boundary described by Eq. ( 36) is clearly symmetric under reflections across the real axis, the number of components of the spectrum is determined by the number of real roots of the equation y 2 = 0, which can be rewritten as:

x 4 + (1 + 2m)g - 1 2g x 3 + 1 16g 2 + g 2 (1 -2m) 2 4 -m - 1 4 x 2 + mg 4 1 2g 2 + 1 -2m x + m 2 16 = 0. ( 37
)
Since Eq. ( 37) is quartic in x, it can have at most four real roots. In that case, the spectrum of the Liouvillian is formed by two disconnected components, each bounded between two real roots of Eq. ( 37), i.e., the system is in phase II. Alternatively, there could be only two real roots (and a pair of complex-conjugated roots that do not contribute to the boundary of the spectrum), in which case the system would be in phase I. The phase transition between these two situations occurs when two real roots coalesce into a double root. The number of real roots of Eq. ( 37) is controlled by its discriminant

∆ 4 = - m 3 128 1 -2(1 -2m)g 2 3 + 216mg 4 . ( 38
)
As discussed, the phase transition corresponds to the merger of two distinct real roots into a double root, which occurs if and only if the discriminant vanishes. Setting ∆ 2 4 = 0 thus gives the critical line in the 1/g vs m plane [START_REF] Haake | Statistics of complex levels of random matrices for decaying systems[END_REF][START_REF] Lehmann | Chaotic scattering: the supersymmetry method for large number of channels[END_REF][START_REF] Janik | Non-Hermitian random matrix models: Free random variable approach[END_REF]:

1 g c = 2 1 -(2m) 1/3 3 . ( 39
)
If ∆ 4 > 0, i.e., g > g c , then Eq. ( 37) has four real roots and we are in phase II. In the converse case ∆ 4 < 0, i.e., g < g c , there are only two real roots and the spectrum belongs to phase I. From Fig. 1(a) and Eq. 39), it is clear that a critical point exists at (m, g) = (1/2, ∞). Below m = 1/2, there is always a finite value of g for which the spectrum splits into two distinct regions. However, the closer m gets to 1/2, the larger are the values of g required for the phase transition to occur, tending to infinity as m → 1/2. Above m = 1/2, the spectrum is connected for all values of g and just stretches indefinitely as we increase g. On the other hand, for g < 1/ 2, the system also belongs to phase I for all values of m. As discussed in Sec. 1.1, below g = 1/ 2 the Hamiltonian contribution to the Lindbladian is strong enough to couple all degrees of freedom, preventing the system to split into two decoupled components.

In phase II, when we increase g, the two regions drift further and further apart from each other, both becoming very thin stripes in the limit g → ∞, one of them aligned with the imaginary axis 36). Right: Gap as function of m in limit g → ∞. The red dots correspond to the numerical results for g = e 10 and black curve to the exact result of Eq. ( 40). In both panels, the numerical results were obtained by exact diagonalization, for N F = 500 and 400 realizations, and match perfectly the analytical predictions.

with increasingly small absolute real part and the other aligned with the real axis with increasingly large absolute real part. Physically, the existence of these two distinct regions of eigenvalues means that the group of modes associated to the region with larger absolute real part will decay much faster than the others. In an intermediate time window, the system will evolve to a metastable state in which only the modes that belong to the region with smaller absolute real part are populated. Eventually, those also fade away and the system reaches the the steady state.

Spectral gap

We now turn to the computation of the spectral gap, obtained from the single-body spectrum through Eq. ( 19). In the limit N F → ∞, the boundary of the single-body spectrum is determined by Eq. ( 36) and, thus, the gap is just the largest real root of Eq. [START_REF] Can | Random Lindblad dynamics[END_REF]. In Fig. 3(a), we plot the gap as a function of g for three different values of m obtained both numerically from ED and exactly from Eq. [START_REF] Can | Random Lindblad dynamics[END_REF], showing perfect agreement between the two. While the expression for the gap can be obtained analytically for any g, its precise functional form is rather complicated and not particularly enlightening. In what follows we will see, however, that simple scaling expressions can be obtained in the limits of weak (g → 0) and strong (g → ∞). The latter case is particularly interesting due to the influence of the decoupling transition.

Weak dissipation

Regardless of the value of m, the gap goes to zero in the limit g -→ 0. This behaviour is expected as, for g = 0, there is just unitary evolution (the Liouvillian becomes simply the von Neumann generator) and so all the eigenvalues lie on the imaginary axis, which means that, by definition, the gap vanishes. In fact, since in Fig. 3(a) the slope of all curves approaches 1 for very small g, we see that, in this limit, Gap ∝ g, a result that can be understood perturbatively. Since a Taylor expansion of the gap around g = 0 yields Gap = ∞ n=1 c n g n , the mentioned scaling behaviour holds unless c 1 = 0. Now, H (U) is a random Hermitian and anti-symmetric matrix and, in the space of all such matrices, the set of degenerate matrices has measure zero. Thus, we can safely apply non-degenerate first-order perturbation theory and conclude that c 1 = min {v † i Γ (U) v i } i , where {v i } i is the set of normalized eigenvectors of H (U) . Γ (U) is a positive semi-definite matrix and thus c 1 = 0 only if there is a v i that belongs to the nullspace of Γ (U) . Clearly, Eq. ( 35) implies that, for m < 1/2, Γ (U) has a nullspace of dimension 2N F -2M = 2N F (1 -2m). However, except for the trivial case m = 0, this is always a set of measure zero and, thus, the probability that one of the v i belongs to the nullspace of Γ (U) is 0. For m > 1/2, the nullspace of Γ is empty. As a consequence, c 1 = 0 and Gap ∝ g, when g → 0, for all m.

Strong dissipation

On the other hand, the limit g -→ ∞ has a nontrivial dependence on m as depicted in Fig. 3(a). In fact, for m < 1/2, the spectral gap closes in this limit (Gap ∝ g -1 ), whereas for m > 1/2 it grows linearly with g → ∞, which suggests a phase transition in the gap at m = 1/2, where an intermediate behaviour is observed, Gap ∝ g -1/3 . The quantity lim g→∞ (Gap/g) can be used as an order parameter for the phase transition as it vanishes for m ≤ 1/2 and acquires a nonzero finite value for m > 1/2. In Fig. 3(b), we plot it as a function of m and compare it with the exact result

lim g→∞ Gap g = 1 2 max 2m -1 2 , 0 (40)
that follows from the Marchenko-Pastur law [START_REF] Marčenko | Distribution of eigenvalues for some sets of random matrices[END_REF]. Indeed, we have lim g→∞ iK (U) /g = Γ (U) , where Γ (U) is a 2N F × 2N F matrix drawn from a real Wishart ensemble with rank 2M , and hence the gap coincides the the hard edge of the Marchenko-Pastur distribution.

The critical behaviour can be traced to the existence or not of zero eigenvalues in the spectrum of Γ (U) and understood perturbatively. Since K (U) = g H (U) /g -iΓ (U) , we can expand the eigenvalues of K (U) /g in powers of 1/g and write the gap as Gap = g ∞ n=0 d n g -n . To zeroth order, the spectrum of K (U) /g is the spectrum of -iΓ (U) . If m > 1/2, Γ (U) is a positive-definite matrix and thus d 0 = min(Im({γ i } i )) > 0 , where {γ i } i is the spectrum of Γ (U) . This justifies the linear growth of the spectral gap with g for m > 1/2. However, for m < 1/2, some of the γ i vanish and, consequently, d 0 = 0. We must look at the next term in the expansion and since, to do that, we need to calculate first order corrections to the zero eigenvalues of Γ (U) , we must resort to degenerate perturbation theory. The corrections are thus given by the eigenvalues of (1/g)w † i H (U) w j , where {w j } j is an orthonormal basis of the nullspace of Γ (U) . Since, however, H (U) is Hermitian, all the corrections to these eigenvalues are real and thus do not affect the gap, leading to d 1 = 0. Only at second order in 1/g do non-zero corrections to the gap arise, which means that Gap = g ∞ n=2 d n g -n . In the large g limit, Gap ∝ g -1 as confirmed by Fig. 3(a). Since m = 1/2 is the critical point, the above expansion in powers of 1/g does not hold. We need, therefore, to resort to Eq. [START_REF] Denisov | Universal Spectra of Random Lindblad Operators[END_REF] to determine the gap's scaling behaviour with g. In Appendix D we perform an asymptotic analysis of the solutions of Eq. ( 36) for large g and show that, at m = 1/2, Gap ∝ g -1/3 . The same procedure can be employed as an alternative to the perturbation theory above, in order to obtain the scaling behaviour of the gap for m > 1/2 and m < 1/2 (and the corresponding prefactors), which is also done in Appendix D.

Comparison with non-quadratic models

We conclude this section by comparing the results obtained here with the gap of a fully-random (non-quadratic) Liouvillian [START_REF] Can | Random Lindblad dynamics[END_REF][START_REF] Can | Spectral Gaps and Midgap States in Random Quantum Master Equations[END_REF][START_REF] Sá | Spectral and steady-state properties of random Liouvillians[END_REF]. Because the single-body gap coincides with the many-body gap the results of the two models can be directly compared. In the weak dissipation regime, the same linear growth with the dissipation parameter is found in both cases and it is the expected perturbative result. The strong dissipation regime is, as we have seen, richer. In the fully-random case, the role of the parameter m is played not by the ratio of dissipation channels to the number of degrees of freedom (2N F here), but by the number of channels itself m → M /2. The regime m < 1/2 is, therefore, inaccessible since M is a positive integer. For fully-random Liouvillians with more than one decay channel (corresponding to m > 1/2) we also observed a linear-in-g growth of the gap, which again is the expected perturbative result for the dissipation-dominated dynamics. For the case of a single decay channel (corresponding to m = 1/2 here) shows the same closing of the gap with ∼ g -1/3 (once one accounts for different normalization conventions). This closing was interpreted as a Zeno-like phase, but it was not understood why it is not observed for more than one decay channel (corresponding to m > 1/2). We can now understand the difference between m > 1/2 and m < 1/2 as a transition in which some degrees of freedom become decoupled from the environment and, thus, protected from dissipation, with m = 1/2 corresponding to the critical transition point. Our findings thus shine a new light to the special role played by fully-random Liouvillians with a single jump operator (corresponding here to m = 1/2). Remarkably, the scaling with m also coincides in both cases (see Appendix D for a computation of the prefactors, which can be compared with Ref. [START_REF] Sá | Spectral and steady-state properties of random Liouvillians[END_REF]) and we conclude that the spectral gap coincides in quadratic and fully random Liouvillians in the mutually accessible regimes (m ≥ 1/2), pointing towards a high degree of universality in dissipative quantum chaos. On the other hand, realistic models with constrained interactions, such as quadratic Liouvillians, contain an additional regime (m < 1/2) of suppressed dissipation.

Steady-state properties

We now turn to the characterization of the steady state, to which the system relaxes in the longtime limit. Because of its Gaussian nature (see Sec. 2), we will base our discussion at the level of the single-particle correlation matrix χ or the thermal-like Hamiltonian Ω, which are related through Eq. ( 25) (proven in Appendix B.3). We will start by studying the spectrum of these singlebody operators, which describe the occupation probabilities of different single-particle states in the long-time limit, see Sec. 4.1. We will focus on the limits of very weak (g → 0) and very strong dissipation (g → ∞), which can be studied perturbatively. To infer the behavior of the steady state as a function of g, we will consider, in Sec. 4.2, the first nontrivial moment of the steady-state distribution-the purity-which captures its degree of mixing. Finally, in Sec. 4.3 we probe the ergodicity of the steady state by analyzing the single-particle level statistics of Ω. As was the case for the single-body spectrum and spectral gap, the results are qualitatively different for m > 1/2 and m < 1/2.

Spectral distribution

The correlation matrix χ is the solution of Eq. ( 26), reproduced here for convenience:

χ K † -Kχ = i g N .
Because of the particle-hole symmetry of χ, Eq. ( 21), its eigenvalues come in pairs {λ i , 1 -λ i }, i ∈ {1, ..., N F }. Consequently, the spectrum of Ω (denoted as {ω i } i ) is formed by the pairs {ω i , -ω i }. The statistical behaviour of λ i and ω i depends on the value of the parameters g and m. Equation (26) can be solved perturbatively for g → 0 and g → ∞, allowing us to analytically study the spectral distribution in these two limiting cases. The details of the perturbative expansion are given in Appendix E, while here we state the final results and work out the consequences for the single-particle effective Hamiltonian Ω.

Weak dissipation

In the weak dissipation limit (g → 0), the eigenvalues of χ to first order in g are determined by Eq. (112) of Appendix E, which we rewrite in the Majorana basis by performing the unitary transformation U, Eq. ( 27):

λ i = v † i N (U) v i v † i N (U) + N (U) T v i , (41)
with v i an eigenvector of H (U) . Since H (U) is Hermitian and anti-symmetric, it can always be diagonalized by a unitary matrix of the form OU, for some orthogonal matrix O. Therefore, in the eigenbasis of H (U) , Eq. ( 41) reads

λ i = N (H) i,i N (H) + (U T U) † N (H) T U T U i,i = N (H) i,i N (H) + SN (H) S T i,i , (42)
where U T U = S, N (H) = l (H) T l (H) * and l (H) = l (U) OU * . From Eq. ( 42), we conclude that the spectrum of χ is composed of N F independent pairs {λ i , 1 -λ i }, with

λ i = M µ=1 l (H) µ,i 2 M µ l (H) µ,i 2 + l (H) µ,i+N F 2 =    1 + M µ l (H) µ,i+N F 2 M µ l (H) µ,i 2    -1 . ( 43
)
The eigenvalues of Ω, in turn, are:

ω i = log    M µ l (H) µ,i+N F 2 M µ l (H) µ,i 2    . ( 44
) Each l (H) µ,2i 2 = Re l (H) µ,i 2 + Im l (H) µ,i
2 is a sum of two random variables with average σ 2 = 1 and variance 2σ 4 = 2 (recall that Re l (H) µ,i and Im l (H) µ,i are sampled from a normal distribution with zero mean and standard deviation σ = 1). We can now resort to the central limit theorem to argue that, for sufficiently large M ,

ω i ≈ log 1 + 1 M Z i+N F 1 + 1 M Z i ≈ 2 M Z i , (45)
where Z and Z are, respectively, a set of 2N F and N F random variables following a normal distribution with unit variance. We conclude that the steady-state single-body spectrum in the weak . The histograms were obtained by exact diagonalization for 100 disorder realizations, with the correlation matrix obtained by the method described in Appendix E. The blue line represents a normal distribution with variance 2/M , as predicted by Eq. ( 51).

dissipation limit is composed of a set of uncorrelated Gaussian random variables with variance 2/M . In Fig. 4, we plot this prediction against the numerical results obtained by exact diagonalization and find perfect agreement. It is clear from Eq. ( 45) that lim M →∞ χ = /2, which means that, in the limits N , M → ∞ with m fixed and g → 0, the steady state of a random quadratic Liouvillian is the fully-mixed state. We will elaborate on this in Sec. 4.2 below. Moreover, we expect single-body Poisson spectral statistics because of the uncorrelated nature of different ω i , a prediction confirmed in Sec. 4.3.

Strong Dissipation

In the strong dissipation limit (g → ∞), m plays a significant role in the statistical behaviour of the spectrum of Ω, with two qualitatively different regimes for m > 1/2 and m < 1/2, which can be traced back to the decoupling transition at m = 1/2.

If m > 1/2, χ is determined, in the eigenbasis of Γ (U) , by Eq. (115) of Appendix E:

χ i, j = w † i N (U) w j γ i + γ j , ( 46
)
where {γ i } i and {w i } i are the set of eigenvalues and orthonormal eigenvectors of Γ (U) , respectively. Since Γ (U) is positive semidefinite, γ i ≥ 0. Note that, for m > 1/2, there are no zero eigenvalues of Γ (U) , and hence γ i +γ j = 0 always holds. From Eq. ( 46) it follows that the eigenvalues of χ are, in general, correlated. Since they are completely determined by the eigenvectors and eigenvalues of a Wishart matrix, we conjecture that they can be related to the Marchenko-Pastur law, although we were not able to do so explicitly. When m < 1/2, χ cannot be obtained from Eq. (115) since γ i + γ j = 0 for some pairs (i, j). More precisely, χ becomes a matrix that acts separately in two subspaces (see Appendix E): subspace N spanned by the eigenvectors of Γ (U) with a corresponding non-zero eigenvalue, w N i , and its complement, subspace N , which is also the nullspace of Γ (U) , spanned by the eigenvectors (w N ) i . Correspondingly, the spectrum splits into two independent components. Antecipating the results found below, we call N the RMT or ergodic sector and N the Poisson or nonergodic sector.

The component of χ that acts in N (i.e., in the sector with no vanishing γ i ), χ N N , is directly obtained from Eq. ( 115) by replacing w i with w N i . The only difference to the case m > 1/2 is thus a reduction of dimensionality. On the other hand, in Appendix E, we show that the eigenvalues of χ N N (i.e., in the sector with zero eigenvalues of Γ (U) ) are given by Eq. ( 125), which can be rewritten as

λ i = (w N ) † i H (U) N N Γ (U) N N -1 N (U) N N Γ (U) N N -1 H (U) N N (w N ) i (w N ) † i H (U) N N Γ (U) N N -1 N (U) N N Γ (U) N N -1 H (U) N N (w N ) i + (w N ) † i H (U) N N Γ (U) N N -1 N (U) N N Γ (U) N N -1 H (U) N N T (w N ) i , ( 47
)
where we used that Γ

(U) N N T = Γ (U) N N and H (U) N N T = -H (U)
N N . Given the similarity of Eqs. ( 41) and ( 47), it is possible to replicate the argument we used for the weak-dissipation case to find the eigenvalues λ i and ω i in this sector. To do so, we change coordinates in both subspaces: in N to the eigenbasis of H (U) N N and in N to the eigenbasis of

Γ (U) N N -1 N (U) N N Γ (U) N N -1
. We denote H (U) N

N in this basis by H(U)

N N and the eigenvalues of

Γ (U) N N -1 N (U) N N Γ (U) N N -1
as α µ . Note that both transformations are unitary in their respective subspaces and independent of H (U) N N , which implies that the entries of H(U) N N remain independent and normally distributed after the transformation. As before, it then follows that the eigenvalues of χ and Ω come in pairs {λ i , 1 -λ i } and {ω i , -ω i }, respectively, with

λ i =    1 + 2M µ α µ H(U) N N i+N F ,µ 2 2M µ α µ H(U) N N i,µ 2    -1
and

ω i = log    2M µ α µ H(U) N N i+N F ,µ 2 2M µ α µ H(U) N N i,µ 2    . ( 48
)
Now, assuming that H(U)

N N i,µ 2 and H(U) N N i+N F ,µ 2
are independent random variables with average σ = 1 and variance 2σ 2 = 2 and applying the central limit theorem, we conclude that

2M µ α µ H(U) N N i,µ 2 ∼ 2M µ α µ   1 + Normal   0, 2 µ α 2 µ µ α µ 2     . ( 49
)
In the limit of very large 2M , this relation becomes

2M µ α µ H(U) N N i,µ 2 ∼ 2M 〈α〉 1 + Normal 0, 1 M 1 + Var (α) 〈α〉 2 . ( 50
)
When inserted into Eq. ( 48), it yields

ω i ≈ log 1 + 1 M Y i+N F 1 + 1 M Y i ≈ 2 M Y i , (51)
where Y and Y are, respectively, a set of 2N F and N F random variables following a normal distribution with variance 1 + Var (α) /〈α〉 2 . In summary, for m > 1/2, the eigenvalues of Ω are . The histograms were obtained by exact diagonalization for 100 disorder realizations, with the correlation matrix obtained by the method described in Appendix E. In (a), the spectrum splits into two independent sectors. The central region corresponds to the nonergodic sector (subspace N ), and the distribution is Gaussian, as shown in the left inset, where we compare the eigenvalues ω centered at their mean and normalized by their standard deviation [ ω = (ω -〈ω〉)/σ ω ] against a normal distribution with unit variance (red line). The remaining eigenvalues belong to the ergodic sector and are perfectly described by Eq. ( 46), see right inset (red line). In (b), there is no nonergodic sector, and all eigenvalues are distributed according to Eq. ( 46). correlated according to RMT, as dictated by Eq. [START_REF] Sá | Lindbladian dissipation of strongly-correlated quantum matter[END_REF]. In contrast, for m < 1/2, the steady-state spectrum splits into two sectors: in the first, the eigenvalues of Ω are still distributed according to Eq. ( 46), but in a space of smaller dimension; the second sector is formed by a set of uncorrelated Gaussian random variables (Poisson level statistics are checked in Sec. 4.3). As in the weak dissipation case, it becomes clear from Eq. ( 51) that lim M →∞ χ N N = /2, signalling that the Poisson sector is fully mixed. In Fig. 5, we plot the spectrum of the single-particle effective Hamiltonian obtained numerically form ED. In Fig. 5(a) (phase II, m < 1/2), the spectrum splits into two well-separated sectors. The central eigenvalues are normally distributed, see left inset, while the larger eigenvalues (in absolute value) follow Eq. [START_REF] Sá | Lindbladian dissipation of strongly-correlated quantum matter[END_REF], see right inset. The parameters 〈α〉 and Var(α) appearing in the variance of the nonergodic eigenvalues can be written in terms of the matrices Γ (U) and Γ

(U) B as 〈α〉 = Tr Γ (U) N N -1 and Var(α) = Tr Γ (U) N N -1 2 + Tr Γ (U) N N -1 2 (Γ (U) B ) N N 2 , ( 52
)
where averaging over the appropriate random ensemble is understood. While 〈α〉 and the first term in Var(α) can be re-expressed in terms of the Marchenko-Pastur distribution, we were unable to evaluate the second term in Var(α). While this prevents a parameter-free comparison with the numerical results, we still confirmed perfect Gaussianity of the nonergodic sector of the steady state in the inset of Fig. 5 53). These results were obtained by exact diagonalization, with the correlation matrix obtained by the method described in Appendix E, for N F = 60 and by averaging over 400 different realizations.

Purity

To study the steady-state spectral distribution away from the limits g → 0, ∞, we look at its moments as a function of g. The first moment is identically one, because of the normalization of probability. The purity, P = Tr ρ 2 NESS , is the lowest nontrivial moment and quantifies the degree of mixing of the steady state.

Since it is possible to express the steady-state density matrix as a function of its correlation matrix χ, the same same applies to the purity:

P = Tr ρ 2 NESS = Tr e C † ΩC Tr e 1 2 C † ΩC 2 = det 1 + e -2Ω 1 2 det (1 + e -Ω ) = det (1 -χ) 2 + χ 2 , ( 53
)
where we applied Eqs. ( 23)-( 25) to obtain the third and fourth equalities.

In the following, instead of the purity itself, we consider the quantity

1 N F log P P FM = log 2 + 1 N F log P, (54)
which we dub reduced purity, and where P FM = 2 -N F is the purity of the fully-mixed state. The reduced purity, coincides with the shifted and rescaled second Rényi entropy, is finite in the limit N F → ∞ and its lower bound (corresponding to the fully-mixed state) is zero. In Fig. 6, we plot the reduced purity as function of g for different values of m and N F = 60. In the limit g → 0, it tends to a fixed value close to zero, for all m. In fact, it is expected to converge to zero in the limit N F → ∞ regardless of m, since, as we proved in the last subsection, lim M →∞ χ = /2 (note that this proof also justifies the slower convergence for smaller values of m depicted in Fig. 6).

As we increase g, two different behaviours of the purity emerge independently of N F in the large-N F limit. If m ≥ 1/2, it increases monotonically with g, stabilizing to a constant in the limit of very large g. The value of this plateau decreases with m. On the other hand, for m < 1/2, the purity initially increases, attaining a maximum at a finite value of g, and then it starts decreasing, converging in the limit g → ∞ to a value smaller than for m > 1/2, and which increases with m. The nonmonotonic behavior of the purity is a consequence of the splitting of the steady-state spectrum into two independent sectors at m = 1/2. In the RMT sector, the steady state has a finite reduced purity, which follows the same functional form as for m > 1/2. On the other hand, in the Poisson sector, the steady-state is fully mixed, as shown in the previous section, and hence has vanishing reduced purity. The competition of the two sectors determines the total purity of the steady-state. Since the nullspace of Γ has dimension 2N F (1 -2m), the Poisson sector becomes increasingly dominant as m decreases. In the limit m → 0, almost all eigenvalues of χ are obtained from the component χ N N and we conclude that the reduced purity converges to zero.

Spectral statistics

We conclude our study of the steady state by studying the single-particle spectral statistics of Ω, which characterize the ergodicity, or lack thereof, of the steady state. We will focus on the distribution of consecutive spacing ratios, rn [START_REF] Oganesyan | Localization of interacting fermions at high temperature[END_REF][START_REF] Atas | Distribution of the ratio of consecutive level spacings in random matrix ensembles[END_REF]. Let s n be the sequence of differences between consecutive eigenvalues of Ω. Then, rn is defined as [START_REF] Oganesyan | Localization of interacting fermions at high temperature[END_REF]

rn = min s n+1 s n , s n s n+1 . ( 55
)
This quantity has the advantage of being independent of the local level density and bounded (0 < rn < 1). If the steady state is ergodic, the spectral correlations of Ω coincide with those of a random matrix of the appropriate symmetry class. On the other hand, if it is nonergodic, it will display Poisson statistics characteristic of uncorrelated random variables. The spacing ratio distributions for all three classes of level repulsion-the Gaussian Unitary Ensemble (GUE), Gaussian Orthogonal Ensemble (GOE), and Gaussian Symplectic Ensemble (GSE)-are well known 41), for 100 different realizations. The colored lines correspond to the analytical results of Eqs. ( 56) and ( 57). We find perfect agreement with the Poisson prediction, confirming nonergodic behavior. and approximated by [59]

P β (r) = 2 Z β (r + r2 ) β (1 + r + r2 ) 3β/2+1 , (56)
with Z 1 = 8/27 (β = 1, GOE), Z 2 = 81 3/4π (β = 2, GUE), and Z 4 = 729 3/4π (β = 4, GSE), while for Poisson statistics it is

P Poi (r) = 2 (1 + r) 2 . ( 57
)
The average spacing ratio 〈r〉 = 1 0 dr P(r) has become a popular measure of ergodicity and of the presence of time reversal in the ergodic phase. Its value for the GOE, GUE, and Poisson statistics is given by, respectively [START_REF] Atas | Distribution of the ratio of consecutive level spacings in random matrix ensembles[END_REF]:

〈r〉 1 = 4 -2 3 ≈ 0.536, 〈r〉 2 = 2 3 π - 1 2 ≈ 0.603, and 〈r〉 Poi = 2 log 2 -1 ≈ 0.386. ( 58
)
(In the following, the GSE will play no role and will not be referred to further.) As a first indication of the influence of m and g on the spectral statistics of Ω, we plot 〈r〉 as a function of g for different values of m in Fig. 7. We observe that, in the limit g → 0, 〈r〉 converges to 〈r〉 Poi (the Poisson value) for all m. Figure 8 corroborates this result by clearly showing that the full distribution of r in the limit g → 0 perfectly matches P Poi (r). In fact, this result follows from Eq. ( 45), where we showed that the spectrum of Ω is Gaussian and composed of N F uncorrelated pairs of eigenvalues.

On the other hand, in the limit g → ∞, two distinct behaviors are observed in Fig. 7: if m > 1/2, r converges to 〈r〉 2 (the GUE value), whereas, if m < 1/2, it converges to different values depending on m. The GUE ratio statistics for m > 1/2 are confirmed in Fig. 9(d). The results for m < 1/2, however, show that there is a crossover from GUE statistics at m = 1/2 to Poisson statistics as m → 0. Again, this behaviour can be understood from our results in Sec. 4.1. There, we showed that, in the Poisson sector (N ), the spectrum of Ω N N is composed of uncorrelated pairs of eigenvalues [Eq. [START_REF] Haake | Statistics of complex levels of random matrices for decaying systems[END_REF]]. This prediction is confirmed in Fig. 9(c). In contrast, in the ergodic sector ( N ), the eigenvalues of Ω N N are in general correlated. In our generic setting, χ (and hence Ω) has no symmetries besides particle-hole symmetry and, consequently, we expect GUE statistics in the ergodic sector, see Fig. 9(b). Remarkably, in the m → 0 limit, the ergodic sector 47) for N F = 500 and 100 distinct realizations. The colored lines correspond to the analytical results of Eqs. ( 56) and [START_REF] Marčenko | Distribution of eigenvalues for some sets of random matrices[END_REF]. We find perfect agreement with the theoretical predictions in each case.

supports GOE statistics instead, see Fig. 9(a), a point we elaborate on further at the end of this section. Moreover, the relative weight of the ergodic and Poisson sectors changes as a function of m, although their dimensions always add up to 2N F . As a consequence, the spectra of Ω N N and Ω N N coexist in the interval m ∈ (0; 1/2) and it is the variation of their relative contributions that causes the crossover from GUE to Poisson statistics as m decreases. In the limit m → 0, the dimension of N converges to 2N F and the spectral statistics of Ω becomes Poissonian.

These two limits (weak and strong dissipation) have implications in the spectral statistics of the steady state at finite g. As depicted in Fig. 7, as g increases from -∞ the spectral correlations of the steady state exhibit a perturbative crossover from Poisson to GUE statistics for all m. A plateau is reached when g is of the order of the inverse mean level spacing of the Hamiltonian. The plateau extends to infinity in the case m > 1/2, in agreement to the previous discussion of the limit g → ∞. On the other hand, for m < 1/2, the plateau is of finite length and, for sufficiently large g, the effects of the nullspace of Γ become relevant. Therefore, the curve 〈r〉 decreases again to meet its expected value in the limit g → ∞, that is, the result of the overlap of the two independent spectra previously mentioned.

Having established the main features of the dependence of the spectral correlations of the steady state on g and m, we now discuss a curious aspect of the strong dissipation limit, namely, that the ergodic sector displays GOE statistics in the limit m → 0, contrarily to a the expectation of GUE statistics discussed above. As shown in Fig. 9(a) and (b), for finite N F , there is a GUEto-GOE crossover as m decreases from m = 1/2 to m → 0. In the thermodynamic limit, GUE statistics are attained for any finite m, while GOE statistics hold only in the strict limit m → 0 5 . In Appendix F, we show that the spectrum of χ N N can be expanded in powers of m and, to lowest nontrivial order (∼ m), it coincides with the spectrum of a matrix from class CI [START_REF] Altland | Nonstandard symmetry classes in mesoscopic normal-superconducting hybrid structures[END_REF], which has the same level correlation as the GOE. We thus conclude that inside the phase m < 1/2 the spectral correlations of the steady state gradually change with m due to two competing phenomena: a RMT-to-Poisson crossover due to the relative weight of the ergodic and nonergodic sectors and a GUE-to-GOE crossover in the ergodic sector.

Conclusions and outlook

In this paper, we studied the single-body spectral and steady-state properties of a fermionic random quadratic Liouvillian. We studied the spectral support and boundary of the single-body Liouvillian spectrum, the spectral gap ruling the approach to the steady-state, and the single-body distribution, purity, and single-body level statistics of the steady state. Our analysis focused on the phase transition observed in the single-body spectrum and its repercussions in the steady-state properties. More precisely, in phase I, the spectral and steady-state properties of quadratic and fully-random Liouvillians are qualitatively similar: the spectrum is formed by a single connected component; the gap grows linearly with dissipation strength; the steady-state purity is monotonic with dissipation strength; and there is a nonergodic-to-ergodic crossover in the steady-state level statistics, the full steady state being ergodic for sufficiently strong dissipation. In phase II, there are qualitative differences: the single-body spectrum splits into two disconnected components at a finite system-environment coupling strength; the spectral gap closes for strong dissipation; the purity is non-monotonic with dissipation strength; and the steady-state decouples into an ergodic and a nonergodic sectors.

In summary, our work identifies a regime of universal random Markovian dissipation but also illustrates the possibility of nonergodic behaviour in quadratic open quantum systems (see also Ref. [START_REF] Sá | Integrable nonunitary open quantum circuits[END_REF]) and the potential to suppress dissipation even in the presence of strong systemenvironment coupling. A natural extension of this work is to ask whether these nonergodic features survive interactions. Their robustness could be addressed, for instance, with the SYK Lindbladian [START_REF] Sá | Lindbladian dissipation of strongly-correlated quantum matter[END_REF][START_REF] Kulkarni | Lindbladian dynamics of the Sachdev-Ye-Kitaev model[END_REF][START_REF] García-García | Keldysh Wormholes and Anomalous Relaxation in the Dissipative Sachdev-Ye-Kitaev Model[END_REF]). Moreover it is also not clear whether the nonergodic features of the steady state survive in the thermodynamic limit N F → ∞ at large but finite g (i.e., whether the red line in Fig. 1(a) is smoothly connected to the rest of the diagram). This interesting question would require a more detailed finite-size scaling analysis and is deferred to future work. Other interesting possibilities are to consider bosonic Liouvillians and non-Markovian generators. Finally, further work is needed to determine whether the properties of stochastic Markovian dissipative models that we described are also present in more realistic models of open quantum systems.

Funding information

A Vectorization and third quantization

In this appendix we explicitly show the equivalence between our vectorization scheme and the more standard third quantization, first introduced in Ref. [START_REF] Prosen | Third quantization: a general method to solve master equations for quadratic open Fermi systems[END_REF].

In Ref. [START_REF] Prosen | Third quantization: a general method to solve master equations for quadratic open Fermi systems[END_REF], it was noted that an orthonormal basis under the Hilbert-Schmidt inner product

〈O 1 |O 2 〉 = 2 -N f Tr O †
1 O 2 can be defined on the (vectorized) space of operators by considering the set of all possible vectors

P α = f α 1 1 f α 2 2 ... f α 2N F 2N F , (59)
where α i ∈ {0, 1} and { f i } i represents a set of 2N F Majorana fermions satisfying the Clifford algebra { f j , f k } = 2δ j,k .They can be obtained from the creation and annihilation operators, c † j and c j , through

f i = 2 j U i j C j , (60)
where U is defined in Eq. [START_REF] Kitaev | Hidden correlations in the Hawking radiation and thermal noise[END_REF]. At this point, creation and annihilation operators can be naturally defined on the vectorized space through

d j P α = δ α j ,1 f j P α = 1 2 f j P α -ΠP α Π f j d † j P α = δ α j ,0 f j P α = 1 2 f j P α + ΠP α Π f j , (61)
where Π = e iπN is the fermionic parity. The Liouvillian becomes quadratic when written as a function of these creation and annihilation operators and its diagonalization proceeds in a very similar manner to our case (see Appendix B).

Our goal is to find a linear transformation from the elements of ã to d j and d † j . In order to do this, we first note that

(Π ⊗ ) d j P α =    -1 2 f j P α Π + P α Π f j = 4N F k=1 Q (1)
jk ãk P α , Π ⊗ Π T P α = P α 1 2 f j P α Π -P α Π f j = - 4N F k=1 Q (2)
jk ãk P α , Π ⊗ Π T P α = -P α , (Π ⊗ ) d † j P α =    -1 2 f j P α Π -P α Π f j = 4N F k=1 Q (2)
jk ãk P α , Π ⊗ Π T P α = P α 1 2 f j P α Π + P α Π f j = - 4N F k=1 Q (1)
jk ãk P α , Π ⊗ Π T P α = -P α , (62)
where

Q (1) = - 1 2 U USJ and Q (2) = - 1 2 U -USJ . ( 63
)
It is now easy to see that, after building the matrices

Q 1 = Q (1) Q (2) = - 1 2 U 0 0 U SJ -SJ , Q -1 = - Q (2) Q (1) = 1 2 U 0 0 U -SJ SJ , (64)
we can relate ã with

D = {d 1 , . . . , d 2N F , d † 1 , . . . , d † 2N F } T through    Q † 1 (Π ⊗ ) D = ã, if Π ⊗ Π T P α = P α , Q † -1 (Π ⊗ ) D = ã, if Π ⊗ Π T P α = -P α . (65)
The two subspaces that appear in the previous equation correspond to the two parity sectors of Ref. [START_REF] Prosen | Third quantization: a general method to solve master equations for quadratic open Fermi systems[END_REF]. In fact, if Π ⊗ Π T P α = P α (respectively, Π ⊗ Π T P α = -P α ), P α contains an even (respectively, odd) number of Majorana fermions.

We can now apply this transformation to Eq. ( 15) to reproduce the results in Ref. [START_REF] Prosen | Third quantization: a general method to solve master equations for quadratic open Fermi systems[END_REF]:

L = - i 2 ã † L ã - i 2 Tr(K) =    -i 2 D † L (1) D D -i 2 Tr(K), if Π ⊗ Π T P α = P α -i 2 D † L (-1) D D -i 2 Tr(K), if Π ⊗ Π T P α = -P α , (66)
where

L (γ) D = Q γ LQ † γ =   H (U) -i 2 N (U) + N (U) T -iδ γ,1 N (U) -N (U) T -iδ γ,-1 N (U) -N (U) T H (U) + i 2 N (U) + N (U) T   , (67)
and γ ∈ {-1, 1} is an eigenvalue of Π ⊗ Π T . The matrices H (U) (Hermitian and anti-symmetric) and N (U) (Hermitian) of the previous expression are defined in Sec. 2.

It is now clear that Eqs. (66) and (67) in the even-parity sector (γ = 1) yield the same result as Eq. ( 22) of Ref. [START_REF] Prosen | Third quantization: a general method to solve master equations for quadratic open Fermi systems[END_REF] 6 . We have thus successfully established the equivalence between both paths to vectorization.

B Spectrum and steady state of fermionic quadratic Liouvillians

In this appendix, we prove the statements made in Sec. 2 about the spectrum and steady-state of quadratic fermionic Liouvillians. In Appendix B.1 we show how to write the single-body matrix of the Liouvillian using our vectorization scheme, proving Eq. [START_REF] Kawasaki | Topological phases protected by shifted sublattice symmetry in dissipative quantum systems[END_REF]. Then, in Appendix B.2 we show how to diagonalize this single-body matrix and, thus, prove Eq. [START_REF] Jaynes | Information Theory and Statistical Mechanics[END_REF]. Finally, in Appendix B.3, we connect the steady-state and the single-body spectrum, proving Eq. ( 26).

i j C † i H i j C j = 1 4
i j H (U) i j w i w j .

B.1 Single-body Liouvillian

In this section, we start by showing how to arrive at Eq. ( 15) from the vectorization of the Lindblad equation, given in Eqs ( 6)- [START_REF] Prosen | Spectral theorem for the Lindblad equation for quadratic open fermionic systems[END_REF]. First, the von Neumann generator,

-i Ĥ, ρ = - i 2 i j H i, j C † i C j ρ -ρC † i C j , (68)
becomes, after vectorization,

- i 2 i j H i, j C † i C j ⊗ -⊗ C T j C † T i |ρ〉 = - i 2 i j H i, j C † i C j ⊗ + ⊗ C † T i C T j |ρ〉 = - i 2 ã † H 0 0 -J H T J ã |ρ〉 . ( 69
)
Particle-hole symmetry, SH T S = -H, was implicitly used in the first equality. One can proceed in a similar way to determine the vectorized version of the jump term in the Lindblad equation:

g µ Lµ ρ L † µ = g µi j l * µ,i l µ, j C i ρC † j → 1 2 ã † 0 -g N S SJ -g JSN 0 ã |ρ〉 . (70)
Finally, the dissipative contribution,

-g µ 1 2 ρ, L † µ Lµ = -g µi j 1 2 l µ,i l * µ, j C † i C j ρ + ρC † i C j (71) becomes - 1 2 ã † gΓ B 0 0 g JΓ T B J ã - i 2 Tr (K) |ρ〉 , (72)
where we used Tr(H) = Tr(SHS) = -Tr(H) = 0. Note that similarly to the case of H, we can always chose the adjoint fermionic creation and annihilation operators to satisfy particle-hole symmetry, which, in the vectorized space, reads as SA T S = -A, where S is defined in Eq. [START_REF] Prosen | Quantization over boson operator spaces[END_REF].

Summing Eqs. ( 69), (70), and (72), we can finally write the Liouvillian single-particle matrix:

L = - i 2 ã † H -i gΓ B -i g N S SJ -i g JSN -J (H + i gΓ B ) T J ã - i 2 Tr (K) , (73)
concluding the proof of Eq. [START_REF] Kawasaki | Topological phases protected by shifted sublattice symmetry in dissipative quantum systems[END_REF].

B.2 Single-body spectrum

To diagonalize the Liouvillian, we must look at transformations of the form ã -→ b = U ãU -1 , so that the canonical anticommutation relations are preserved. Let us consider the matrix exp{-i with ∆ = -S∆ T S (note that ∆ can always be chosen to satisfy particle-hole symmetry). It follows that

i j 1 2 ã † i ∆ i j ãj , ã † k = i ã † i ∆ ik , (74)
which implies that

e i 2 i j ã † i ∆ i j ãj ã † e -i 2 i j ã † i ∆ i j ãj = ã † e i∆ . ( 75
)
Defining b = R -1 ã and b T = ã † R, where R = e i∆ , we find:

L = - i 2 b T R -1 LR b - i 2 Tr (K) . ( 76
)
We have reduced the diagonalization of the Liouvillian to the determination of the matrix e -i∆ that diagonalizes L. Note that in general b T = b † , all we know is that { bi , b j } = δ i j . Also, the particle-hole symmetry of ∆ leads to the following restriction for the choice of e -i∆ :

SR T S = Se i∆ T S = e -i∆ = R -1 . ( 77
)
Any R that diagonalizes L can be made to satisfy Eq. ( 77) by reordering columns, if necessary. To see this, suppose first that v L is a left eigenvector of L, i.e., L T ν L = τ ν ν L . Then,

L Sν L = -S L T ν L = -τ ν Sν L . (78)
If {τ (i) } i is a set of 2N F eigenvalues of L with, for example, non-positive imaginary part and {ν

with

R 1 = ν (1) R , . . . , ν (2N ) R , Sν (1)
L , . . . , Sν

(2N ) L ( 80
)
and

P 23 =    0 0 0 0 0 0 0 0 0 0 0 0    (81)
does indeed satisfy Eq. (77). At this point, all that is left to complete the diagonalization of the Liouvillian is to determine the matrix R 1 explicitly. To this end, we first note that L can be made block upper triangular after application of the transformation U = 1 2 -JS JS :

U -1 LU = K 2i gΓ B 0 K † . ( 82
)
Note that the previous equation implies the eigenvalues of L and K coincide. It is straightforward to see that the matrix L can be finally diagonalized by application of another change of basis, implemented by

V = R X R † -1 0 R † -1 , ( 83
)
where X is the solution of

X K † -K X = 2i gΓ B . (84)
Indeed,

(U V ) -1 L(U V ) = R -1 KR 0 0 R -1 KR † = D K 0 0 D * K , (85)
Recalling that K = Hi gΓ , with Γ a positive semidefinite matrix, all the eigenvalues of K are non-positive, which implies that

v (1) R , . . . , v (2N F ) R = 1 2 R JSR . ( 86
)
To calculate the corresponding left eigenvectors we note that, after choosing the correct normalization, X L X R = , where X L (X R ) is a matrix whose rows (columns) are the left (right) eigenvectors of L. Thus, inverting U V and keeping the first 2N F rows gives us the desired result:

v (1) L , . . . , v (2N F ) L = 1 2 1 + X T R T -1 JS -X T R T -1 (87)
and, therefore,

R 1 = 1 2 R S 1 + X T R T -1 JSR -J -X T R T -1 . (88)
Finally, the Liouvillian assumes the form

L = - i 2 b T D L b - i 2 Tr (K) , (89)
where D L = R -1 LR obeying the particle-hole symmetry:

S D L S = SR T L T R -1 T S = -R -1 LR = -D L . (90)
Therefore, we can write Diagonal (D L ) = {β (1) , -β (1) , β (2) , -β (2) }, where β (1) and β (2) are both vectors of N F distinct eigenvalues of K. Actually, Eq. ( 89) can be further simplified by noting that it follows from Eq. (77

) that Sb = SR -1 ã = ã † R T = b , (91)
which means that b = {b (1) , b (1) , b (2) , b (2) } and b = {b (1) , b (1) , b (2) , b (2) }, where b (1) and b (2) are both vectors of N F annihilation operators and b (1) and b (2) of the corresponding creation operators. Equation (89) can be finally reduced to

L = - i 2 2 k=1 N F j=1 β (k) j b (k) j b (k) j -b (k) j b (k) j - i 2 Tr (K) = - i 2 2 k=1 N F j=1 β (k) j b (k) j b (k) j + i 2 2 k=1 N j=1 β (k) j -Tr (K) = - i 2 2N F j=1 β j b j b j , (92)
where b = {b (1) , b (2) } and β = {β (1) , β (2) }. The levels -iβ j /2 are the single-particle eigenvalues of the Liouvillian. The many-body eigenvalues are immediately obtained as

Λ n = - i 2 j n j β j , (93)
where n = {n 1 , . . . , n 2N F } and n j = 0, 1.

B.3 Steady state

From the discussion above, all elements of β have nonpositive imaginary part, which implies that if Im(β i ) = 0 holds for all β i , then the steady state, which satisfies L(ρ NESS ) = 0, is unique and it is annihilated by all b i . To see this, we note that L † ( ) = 0 or, in vectorized notation, 〈 | L = 0. Thus, 〈 | and |ρ NESS 〉 form a biorthogonal left-right eigenvector pair. Using the anticommutation relations of b i and Eq. ( 92), this implies that

〈 | b i L |O〉 = -〈 | Lb i |O〉 + i 2 β i 〈 | b i |O〉 , (94)
or equivalently,

〈 | b i L - i 2 β i |O〉 = 0, (95)
for all vectorized operators O. Since the single-body spectrum of the Liouvillian is contained in the lower-half plane and, by assumption Im(β i ) = 0, we have Li 2 β i = 0 and, hence, Eq. (95) implies that, for all i, 〈 | b i = 0 and, consequently, b i |ρ NESS 〉 = 0.

For quadratic systems, the steady state is Gaussian and completely determined by its 2N F ×2N F correlation matrix,

χ i j = Tr ρ NESS C i C † j = 2N k=1 S jk 〈 | ãi ãk |ρ NESS 〉 , (96)
where i ∈ {1, ..., 2N F }. Remarkably, χ i j is also easily determined in terms of the single-body matrices K and Γ B . Indeed, we can write ã as a combination of b and b

χ i j = 2N F k=1 4N F α,β=1 S jk R iα R kβ 〈 | bα bβ |ρ NESS 〉 , (97)
and use the relations

〈 | b i b j |ρ NESS 〉 = 〈 | b i b j |ρ NESS 〉 = 〈 | b i b j |ρ NESS 〉 = 0 and 〈 | b i b j |ρ NESS 〉 = δ i j (98)
along with Eq. (88) to simplify Eq. (97):

χ i j = 2N k,n=1 S jk R (1) in R (2) kn = 1 2 (1 + X ) i j , (99)
with R (1) , . . . , R (4) particle-hole blocks of the matrix R 1 , i.e., R 1 = R (1) R (2) R (3) R (4) , and X defined through Eq. ( 84). The single-particle matrix X thus completely determines the steady state.

C Resolvent of antisymmetric Hermitian random matrices

In this appendix, we show that the resolvent (also known as the Green's function) of anti-symmetric random Gaussian matrices (class D) coincides with that of GUE and GOE matrices (classes A and AI, respectively), in the large-N F limit. We will use the method of moments, in which we find the exact leading-order behaviour of the pth moment of a matrix Q and then infer the resolvent through the relation

G(z) = 1 N F Tr 1 z -Q = 1 N F ∞ p=0 1 z p+1 〈Tr Q p 〉 . ( 100
)
Let us first compute the resolvent of a GUE matrix Q, with probability distribution P(Q) ∼ exp{-N F /2 Tr M 2 }. All odd moments vanish. All even moments can be related to the second moment (the propagator),

〈Q ab Q cd 〉 = 1 N F δ ad δ bc (101)
through Wick's theorem, i.e., by summing over all possible pair contractions of the indices a, b, c, etc. For instance, the second moment is trivially

Tr Q 2 = a,b 〈Q ab Q ba 〉 = N F , (102)
while the first nontrivial moment, the fourth, is given by:

Tr Q 4 = a,b,c,d 〈Q a b Q bc Q cd Q d a 〉 = a,b,c,d (〈Q a b Q bc 〉 〈Q cd Q da 〉 + 〈Q ab Q da 〉 〈Q bc Q cd 〉 + 〈Q ab Q cd 〉 〈Q bc Q da 〉) = N F 2 + 1 N 2 F . ( 103
)
We can associate to each Wick contraction a perfect matching of the p matrices Q in the trace. Schematically, for p = 4:

Tr Q 4 = 〈Tr QQQQ〉 = 〈Tr QQQQ〉 + 〈Tr QQQQ〉 + 〈Tr QQQQ〉.
Then, non-crossing (or planar) perfect matchings (e.g., the ones corresponding to the first two contractions in the fourth moment) contribute with a factor N F to the trace, while each crossing in the perfect matching suppresses the contribution of that contraction by 1/N 2 F (e.g., the third contraction in the fourth moments has a single crossing). The number of non-crossing perfect matching of p elements is well-known to be the pth Catalan number, C p , and, hence, we conclude that 〈Tr Q p 〉 = N F C p in the large-N F limit. Equation (100) can then be resummed.

Let us now turn to symmetric and antisymmetric random matrices, Q ± = (Q ± Q T )/ 2. Q + is a GOE matrix, while Q -belongs to class D, as considered in the main text. As before, we can express each moment of Q ± in terms of the propagator of Q using Wick contraction. The second moment is given by

Tr Q 2 ± = Tr Q 2 ± Tr QQ T = Tr Q 2 ± ab Q ab Q ab = N F ± 1. ( 104
)
We see that the contribution due to (anti)symmetrization is subleading in 1/N F . This carries over to higher moments: any contraction of Q and Q T is suppressed by 1/N F . We note that, incidently, the corrections due to (anti)symmetrization are less suppressed than those arising due to nonplanarity. We can proceed similarly for the fourth moment, which is given by:

Tr Q 4 ± = 1 2 Tr Q 4 ± 4 Tr Q 3 Q T + 2 Tr Q 2 (Q T ) 2 + Tr QQ T QQ T . ( 105
)
The only nonsupressed contractions (either by symmetrization or nonplanarity) are

1 2 〈Tr QQQQ〉 + 〈Tr QQQ T Q T 〉 = Tr Q 4 . ( 106
)
For the sixth moment, we find

1 4 〈Tr QQQQQQ〉 + 3〈Tr QQQQQ T Q T 〉 = Tr Q 6 , (107)
and similarly for higher moments. We conclude that Tr Q p ± = 〈Tr Q p 〉 and, hence, their resolvents coincide.

D Spectral gap from the asymptotic analysis of Eq. (36)

In this appendix, we extract the leading-order behavior of the gap from an asymptotic analysis of Eq. ( 36), which we reproduce here for convenience,

y 2 = - m 4g x - g/2 1 -4g x + m 4x - 1 4g 2 ,
in the limits g → ∞ and g → 0. This procedure can be employed as an alternative to the perturbation theory of Sec. 3.2 and yields not only the leading scaling with g but also the exact prefactor.

Let us fix the notation first. We say that A(g) B(g) if the leading order terms of A and B are equal (more rigorously, lim g→∞ A(g)/B(g) = 1). We are therefore interested in computing C and α in x g (g) -C g α , where x g is (minus) the gap and C ≥ 0. We can now proceed as described before, setting y = 0 and replacing x by -C g α in order to estimate the asymptotic behavior in Eq. [START_REF] Denisov | Universal Spectra of Random Lindblad Operators[END_REF]. After some manipulations becomes:

± m 4C g -(α+1)/2 + 4mC g (α+1)/2 1 2 -m g -C g α - m 4C g -α - 1 4 g -1 . ( 108
)
If none of the exponents of g in the previous equation match, then every term must vanish identically and the equation becomes trivial. There are five different values of α for which some of the exponents match: α = 1, α = 0, α = -1/3, α = -1, and α = -3.

• At α = -3 the leading order term in the right hand side is -m/(4C)g 3 and it is unmatched in the left hand side, which implies that m = 0. Since the gap at m = 0 trivially vanishes, we must rule out α = -3.

• Similarly, for α = 0, we obtain that m = 1/2 and then C = 0, which also means that α = 0 is not allowed.

• For α = -1, Eq. ( 108) becomes

1 2 -m - m 4C = 0 ⇔ C = m/4 1/2 -m . ( 109
)
Note that this equation is only valid for m < 1/2 as C must be positive.

• If we instead assume α = -1/3, then m = 1/2 and

± 4mC = -m/(4C) ⇔ C 3/2 = ± m/8. Since C is positive, C = m 1/3 /4 = 2 -1/3 /4. • Last but not least, for α = 1, ± 4mC = 1/2 -m -C ⇔ C = 1 ± 2m 2 /2. This solution
exists for all values of m.

By inspection of the graphics present in Fig. 2 of the main text, we can assign one of the asymptotic behaviours above to each of the intersection points of the boundary with y = 0.

• For m < 1/2, the system is in phase II and thus there are two disconnected regions in the spectrum and four intersection points. The ones that delimit the region of eigenvalues with smaller real part correspond to the two solutions for the case α = 1. The other two must correspond to the other solution available, at α = -1, which means that they have the same scaling behaviour with g. Since the gap is given by (minus) the intersection point with largest real part, we conclude that for m < 1/2, Gap ∼ m/4

1/2-m g -1 .

• For m ≥ 1/2 there is just a single connected region (phase I) and thus there are just two intersection points. At m = 1/2, Gap ∼ 2 -1/3 g -1/3 /4.

• For m > 1/2, Gap ∼ 1 -2m 2 g/2, in agreement with the Marchecko-Pastur formula, Eq. ( 40).

The behaviour of x g (g) in the limit g → 0 can also be extracted from Eq. (108) (note that now A(g) B(g) means that lim g→0 A(g)/B(g) = 1). The non-positive solutions of α cannot represent the asymptotic behaviour of the gap, since it is clear from the expression for K(g) that the gap converges to 0 in the limit g → 0. Therefore, we are left with α = 1, for which the condition ± m/(4C) = -1/4 -m/(4C) must be verified. Since C > 0, this equation becomes m/(4C) = 1/2 ⇔ C = m and thus Gap ∼ mg, in agreement with perturbation theory.

E Perturbative steady state spectrum

In this appendix, we compute perturbatively the steady-state correlation matrix and its spectrum, in the limits of very weak and very strong dissipation. From Eqs. (84) and (99) in Appendix B.3, we can write the following equation for the correlation matrix of the steady state:

χ K † -Kχ = i g N . ( 110
)
We are interested in determining the solution to Eq. (110) in the limits g -→ 0 and g -→ ∞.

E.1 Weak dissipation

In the weak dissipation limit, we expand χ in powers of g, χ = χ (0) + χ (1) g + χ (2) g 2 + . . . , and plug it in Eq. ( 110). Comparing terms order by order, we obtain that:

χ (0) , H = 0,
i χ (1) , H + Γ χ (0) + χ (0) Γ = N . (111)

Let {v i } i be an orthonormal eigenbasis of H (the element v j is to be understood as a column vector). Since χ (0) commutes with H, {v i } i can always be chosen to also form an eigenbasis of χ (0) . Thus, denoting the eigenvalues of H and χ (0) as ε i and λ i , respectively, the second equality in Eq. ( 111) yields

-i χ (1) i,i ε i -χ (1) i,i ε i + 2λ i v † i Γ v i = v † i N v i ⇐⇒ λ i = 1 2 v † i N v i v † i Γ v i , ( 112
)
where χ

(1)
i,i is short-hand notation for v † i χ (1) v i . Note that we used that both H and χ (0) are Hermitian and, consequently, their spectrum is real. We thus find that lim g→0

χ(g) = χ (0) = i λ i v i v † i , (113)
with λ i given by Eq. ( 112). Note that, in general, we could have to be more careful here. If some of the v † i Γ v i were zero, this result would not be valid and we would have to look at terms of higher order in g in the power expansion of Eq. (110). This could only happen if v i belongs to the nullspace of Γ , since Γ is a positive semidefinite matrix. However, since in the present paper we work with random and independent H and Γ , the probability of v † i Γ v i = 0 is zero as the nullspace of Γ is a set of measure zero.

E.2 Strong dissipation

In the strong dissipation limit, χ is expanded in powers of 1/g, χ = χ (0) + χ (1) /g + χ (2) /g 2 + . . . and the comparison of terms of the same order in Eq. (110) yields Γ χ (0) + χ (0) Γ = N , Γ χ (1) + χ (1) Γ = i χ (0) , H , Γ χ (2) + χ (2) Γ = i χ (1) , H . (114) Let {γ i } i and {w i } i be the set of eigenvalues and orthonormal eigenvectors of Γ , respectively. Then, it is easy to see that from the first equality in Eq. (114), we can obtain

γ i + γ j χ (0) i, j = w † i N w j ⇔ χ (0) i, j = w † i N w j γ i + γ j , if γ i + γ j = 0, (115)
where χ (0) i, j = w † i χ (0) w j . Since Γ is positive semidefinite, γ i ≥ 0. If all γ i are positive, then Eq. (115) completely determines χ (0) = lim g→∞ χ(g).

However, if γ i +γ j = 0, Eq. (115) trivially holds for all χ (0) i, j and we have to look at higher order terms to compute it. Note that the nullspace of the matrix Γ is contained in the nullspace of N . In which implies:

2 β H N N jβ γ -1 β H N N β j λ j = αβ γ -1 α + γ -1 β H N N jβ χ (0) N N βα H N N α j . (123)
In the preceding equations, the Greek indices label the eigenvectors that do not belong to the nullspace of Γ and H N N jβ = (w N ) † j H w N β . We can now use Eq. (117) to simplify the above result,

2 β H N N jβ γ -1 β H N N β j λ j = αβ H N N jβ N N N βα γ α γ β H N N α j , (124)
and, solving for λ j :

λ j = 1 2 (w N ) † j H N N Γ -1 N N N N N Γ -1 N N H N N (w N ) j (w N ) † j H N N Γ -1 N N H N N (w N ) j . ( 125
)
In conclusion, we found that

lim g→∞ χ(g) = αβ w N α w N † α N N N w N β γ α + γ β w N † β + i λ i (w N ) i (w N ) † i . (126)
Considering the case m < 1/2, we now show that Eq. ( 117) can be considerably simplified. Noting that N = Γ + Γ B , it is easy to see that

χ (0) N N i j = 1 2 δ i, j + w N † i Γ B w N j γ i + γ j . ( 127
)
If we now decompose l µ, j = l R µ, j + il I µ, j into its real and imaginary parts, Γ and Γ B can be rewritten as

Γ j,k = M µ l R µ, j l R µ,k + l I µ, j l I µ,k = 2M µ y µ, j y µ,k (Γ B ) j,k = i M µ l R µ, j l I µ,k -l I µ, j l R µ,k = 2M µν y µ, j J µ,ν y ν,k , (128)
where J = 0 i -i 0 and y = l R l I . Since we assume that 2M < 2N F , and due to Eq. (127), only the subspace spanned by the eigenvectors of Γ associated to non-zero eigenvalues is relevant. We therefore restrict all the following analysis to this subspace. We start by changing coordinates to the basis defined by y α , writing the eigenvectors of Γ as w N i = α w N i,α y α .

The action of Γ on a vector x j = α xα y α, j that belongs to this subspace can thus be written as

implying that Γµ,α = k y µ,k y α,k is the matrix representation of Γ in the new basis 7 . Similarly, it is easy to see that Γ B in the new basis reads ΓB = J Γ . Since the matrix that implements this change of basis is not orthonormal, the metric after the transformation is no longer the identity. In fact, it becomes Γ , which means that it must be included in all inner products computed in this basis. This allows us to rewrite Eq. (127) as

χ (0) N N i j = 1 2 δ i, j + w N T i Γ J Γ w N j γ i + γ j = 1 2 δ i, j + γ i γ j w N T i J w N j γ i + γ j . ( 130
)
Note that, since the eigenvectors w N i are normalized, we have

w N T i Γ w N i = 1 ⇐⇒ w N T i w N i = 1 γ i , ( 131
)
which means that the norm of w N i under the usual inner product is 1/ γ i . For convenience, we perform the transformation w N i → w N i / γ i , leading to our final expression for χ (0) N N i j :

χ (0) N N i j = 1 2 δ i, j + γ i γ j w N T i J w N j γ i + γ j . (132)
F GOE statistics for the steady state in the limit g → ∞, m → 0

In this subsection, we are interested in studying the limit m → 0 of Eq. (132) for the case of a random Liouvillian (sampled as described in Sec. 2). We show that the steady state supports GOE statistics in this limit. Assuming M and N F are sufficiently large, but such that M 2N F , we can write

Γ = + M 2N F Γ (1) , (133)
where, from the central limit theorem, the spectrum of Γ (1) , γ

i (with

γ i = 1+ mγ (1)
i ), is convergent in the limit M → ∞ and N F → ∞ and of order one. Plugging this in Eq. ( 132) and Taylor expanding it up to second order in m = M /(2N F ), we obtain:

χ (0) N N i j = 1 2 δ i, j + 1 2 wT i J w j - m 16 γ (1)
i -γ where we defined χ (w)

N N = 1 2 + 1 2 J - m 16 Γ (1) 2 J + J Γ (1) 2 -2 Γ (1) J Γ (1) . ( 135
)
We conclude that the diagonalization of χ (w)

N N provides the eigenvalues of χ (0)

N N up to first order in m.

Figure 1 :

 1 Figure 1: Schematic representation of the main results of this paper, the single-particle spectral and steady-state properties of random quadratic Liouvillians. (a) Phase diagram in the 1/g versus m plane. Phases I and II are separated by the critical line g c (m) (dashed line). In the limit g → ∞, the phase transition occurs at m = 1/2, with the spectral gap vanishing for m < 1/2 (red line). (b, c) Single-body spectrum in the complex plane, with a single connected component in phase I (b) and two disconnected components in phase II (c). The spectral boundary (red line) can be obtained analytically. (d, e) Steadystate properties for g → ∞ and m > 1/2 (phase I). The steady-state has a single sector with the single-particle effective Hamiltonian eigenvalues, ω, distributed according to RMT (d) and GUE statistics for the spacing ratio r (e). (f, g) Steady-state properties for g → ∞ and m < 1/2 (phase II, red line). The steady-state spectrum splits into an ergodic and a nonergodic sectors, in which the effective-Hamiltonian eigenvalues follow, respectively, RMT and a normal distribution (f). The r statistics are Poisson in the nonergodic sector and interpolate from GUE to GOE in the ergodic sector as m → 0 (g).

Figure 2 :

 2 Figure 2: Single-body spectrum of the Liouvillian for m = 0.2 and g = 1 (a), g ≈ 5.2 (b), and g = 20 (c). The the 2N F blue points are obtained numerically by exact diagonalization (for a single realization with N F = 4000) and are bounded by the red curve, which is parametrized by Eq. (36).

Figure 3 :

 3 Figure 3: Average spectral gap as a function of the dissipation strength g and number of jump operators m. Left: Gap as a function of g for three different values of m: m = 0.25, 0.50, and 0.75. The full colored curves give the numerical results, the dashed black curve gives the exact result in the limit N F → ∞ obtained by solving Eq. (36). Right: Gap as function of m in limit g → ∞. The red dots correspond to the numerical results for g = e 10 and black curve to the exact result of Eq. (40). In both panels, the numerical results were obtained by exact diagonalization, for N F = 500 and 400 realizations, and match perfectly the analytical predictions.

Figure 4 :

 4 Figure 4: Spectrum of the effective Hamiltonian Ω for N F = 500, m = 0.25 (a) and m = 0.75 (b), and g = 10 -10. The histograms were obtained by exact diagonalization for 100 disorder realizations, with the correlation matrix obtained by the method described in Appendix E. The blue line represents a normal distribution with variance 2/M , as predicted by Eq. (51).

Figure 5 :

 5 Figure 5: Spectrum of the effective Hamiltonian Ω, for N F = 500, m = 0.25 (a) and m = 0.75 (b), g = 10 6. The histograms were obtained by exact diagonalization for 100 disorder realizations, with the correlation matrix obtained by the method described in Appendix E. In (a), the spectrum splits into two independent sectors. The central region corresponds to the nonergodic sector (subspace N ), and the distribution is Gaussian, as shown in the left inset, where we compare the eigenvalues ω centered at their mean and normalized by their standard deviation [ ω = (ω -〈ω〉)/σ ω ] against a normal distribution with unit variance (red line). The remaining eigenvalues belong to the ergodic sector and are perfectly described by Eq. (46), see right inset (red line). In (b), there is no nonergodic sector, and all eigenvalues are distributed according to Eq. (46).

  (a). In Fig.5(b) (phase I, m > 1/2), all eigenvalues belong to a single ergodic sector.

Figure 6 :

 6 Figure 6: Reduced purity as a function of the dissipation strength g for seven different values of m = 0.10-0.95, as given by Eq. (53). These results were obtained by exact diagonalization, with the correlation matrix obtained by the method described in Appendix E, for N F = 60 and by averaging over 400 different realizations.

Figure 7 :

 7 Figure 7: Mean level-spacing ratio 〈r〉 as a function of as a function of the dissipation strength g for six different values of m = 0.1-0.75. These results were obtained by exact diagonalization, with the correlation matrix obtained by the method described in Appendix E, for N F = 60 and by averaging over 400 different realizations.

Figure 8 :

 8 Figure 8: Level spacing ratio distribution r for N F = 500, m = 0.25 (a) and m = 0.75 (b), in the limit g → 0. The histograms were obtained by exact diagonalization, with the correlation matrix obtained by Eq. (41), for 100 different realizations. The colored lines correspond to the analytical results of Eqs. (56) and[START_REF] Marčenko | Distribution of eigenvalues for some sets of random matrices[END_REF]. We find perfect agreement with the Poisson prediction, confirming nonergodic behavior.

Figure 9 :

 9 Figure 9: Level spacing ratio distribution r in the RMT sector for m = 0.025, the RMT (b) and Poisson (c) sectors for m = 0.25, and for m = 0.75 (d), in the limit g → ∞. The histograms in panels (a), (b), and (d) were obtained by exact diagonalization of Eq. (46) for N F = 3000 and 1000 different realizations. The histogram of panel (b) was obtained by solving Eq. (47) for N F = 500 and 100 distinct realizations. The colored lines correspond to the analytical results of Eqs. (56) and[START_REF] Marčenko | Distribution of eigenvalues for some sets of random matrices[END_REF]. We find perfect agreement with the theoretical predictions in each case.

  This work was supported by Fundação para a Ciência e a Tecnologia (FCT-Portugal) through grants No. SFRH/BD/147477/2019 (LS) and UID/CTM/04540/2019 (PR). TP acknowledges ERC Advanced grant 694544-OMNES and ARRS research program P1-0402. This project was funded within the QuantERA II Programme that has received funding from the European Union's Horizon 2020 research and innovation programme under Grant Agreement No 101017733. ADL acknowledges support by the ANR JCJC grant ANR-21-CE47-0003 (TamEnt).

R

  } i are the corresponding left and right eigenvectors, respectively, then one can easily check that the matrix R = R 1 P 23 ,

j

  Γ i, j x j = 2M µ y µ, j α k y µ,k y α,k xα = 2M µ y µ, j α Γµ,α xα ,

  j +O m 3/2 = wT χ (w) N N w+O m 3/2 , (134)

With the important caveat that here we studied the single-particle properties, whereas for non-quadratic models many-body properties have to be considered. Nonetheless, unconstrained fully-random Liouvillians also exhibit a connected spectrum, the gap (which is the same in the single-and many-body case) has the exactly same scaling as here, and the steady state exhibits an nonergodic-to-ergodic crossover and corresponding change in purity as function of g.

We assume a generic case where Liouvillian has no non-trivial Jordan blocks. For the most general treatment see, e.g., Ref.[START_REF] Prosen | Spectral theorem for the Lindblad equation for quadratic open fermionic systems[END_REF].

Also referred to as covariance matrix.

Note that Eq.[START_REF] Liao | Many-Body Level Statistics of Single-Particle Quantum Chaos[END_REF] differs from that in Refs.[START_REF] Janik | Non-hermitian random matrix models[END_REF] by some numerical factors, which have their origin in the different normalizations used. Using our conventions, H r is sampled from the probability distribution ∼ exp{-N F Tr H 2 r } and A from the distribution ∼ exp{-N F Tr AA T }.

Note that m → 0 corresponds to any value of M that is either fixed or grows slower than N F . The limit m → 0 therefore still covers a wide range of physical systems.

Note that these equations seem to differ by some factors, but they are just a consequence of different definitions of the Lindblad equation (a factor of 2 in the dissipation part) and the fact that Ĥ = 1

i j ã † i ∆ i j ãj }

Note that in the case where the entries y µ,k are independent real random variables, this transformation establishes an equivalence between the non-zero eigenvalues of a Wishart matrix with dimension 2N F > 2M and 2M degrees of freedom with the spectrum of another Wishart matrix with dimension 2M and 2N F degrees of freedom.

Acknowledgements

We thank Sergey Denisov and Kristian Wold for very fruitful discussions.

fact, since N and N S = SN T S are positive semidefinite matrices, w † i N w i ≥ 0 and w † i N S w i ≥ 0. If w i belongs to the nullspace of Γ , then w † i Γ w i = 0 ⇒ w † i N w i = 0 ⇒ N w i = 0. Suppose, for example, that, in Eq. ( 115), γ j = 0. Then, because of the property we have just showed, χ (0) i, j = 0. This means that, introducing the projector onto the nullspace of Γ , P N , and onto its orthogonal complement, we have P N , χ (0) N N = P N χ (0) P N = 0 and χ (0)

N N = P N χ (0) P N = 0. This is equivalent to the statement that χ (0) is block diagonal,

where χ (0)

N N = P N χ (0) P N and χ (0)

N N can be determined directly from Eq. (115),

with (w N ) i an eigenvector of Γ that does not belong to its nullspace, whereas, to compute χ (0)

N N , we must resort to the other equalities in Eq. (114). From the second one, we can write

where we used that

N N and H N N share the same eigenbasis but it does not allow us to compute its eigenvalues. To determine them, we start by observing that, from the third equality in Eq. (114),

Up to this point, we were completely free to choose a basis {(w N ) i } i for the nullspace of Γ as long as we kept it orthonormal. However, we will fix it now in such a way that H N N becomes a diagonal matrix:

i (note that it is always possible since H N N is Hermitian). Plugging this in Eq. (119), we arrive at:

Setting j = k, we can eliminate the new variable χ

(1)

N N and write a closed equation for the variables already present in Eq. ( 118):

We can now solve the second equality in Eq. (118) for χ

(1) N N (and its Hermitian conjugate), insert it in the last equality and use the fact that, since χ (0)

We now resort to first-order perturbation theory to simplify the previous equation while keeping it exact to first order in m. The eigenvalues of J are either -1 or 1, which means that the spectrum of χ (w) N N splits into two regions, one close to 0 and the other to 1. J is diagonalized by the matrix

Γ

Γ

Γ

the first order correction to the positive eigenvalues of J is given by the spectrum of

Γ

where

12 is a generic complex symmetric matrix. In the limit of m → 0 the spacings of the eigenvalues of χ (w) N N are thus determined by the spacing of the eigenvalues of V V † . Since the eigenvalues of V V † coincide with the positive eigenvalues of the chiral matrix

V V † belongs to class CI [START_REF] Altland | Nonstandard symmetry classes in mesoscopic normal-superconducting hybrid structures[END_REF][START_REF] Haake | Quantum signatures of chaos[END_REF], which has the same level statistics as the GOE. This finally explains the observation that the level correlations of Eq. (132) approach those of the GOE in the limit m → 0.