Spectral and steady-state properties of fermionic random quadratic Liouvillians

We study spectral and steady-state properties of generic Markovian dissipative systems described by quadratic fermionic Liouvillian operators of the Lindblad form. The Hamiltonian dynamics is modeled by a generic random quadratic operator, i.e., as a featureless superconductor of class D, whereas the Markovian dissipation is described by $M$ random linear jump operators. By varying the dissipation strength and the ratio of dissipative channels per fermion, $m=M/(2N_F)$, we find two distinct phases where the support of the single-particle spectrum has one or two connected components. In the strongly dissipative regime, this transition occurs for $m=1/2$ and is concomitant with a qualitative change in both the steady-state and the spectral gap that rules the large-time dynamics. Above this threshold, the spectral gap and the steady-state purity qualitatively agree with the fully generic (i.e., non-quadratic) case studied recently. Below $m=1/2$, the spectral gap closes in the thermodynamic limit and the steady-state decouples into an ergodic and a nonergodic sector yielding a non-monotonic steady-state purity as a function of the dissipation strength. Our results show that some of the universal features previously observed for fully random Liouvillians are generic for a sufficiently large number of jump operators. On the other hand, if the number of dissipation channels is decreased the system can exhibit nonergodic features, rendering it possible to suppress dissipation in protected subspaces even in the presence of strong system-environment coupling.


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 Gorini-Kossakowski-Sudarshan-Lindblad form [1][2][3], or just Lindblad form for short: 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 of 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 [4,5]. Nevertheless, it still finds many valuable applications in various subjects, namely in quantum optics and quantum computation [6].
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 [7][8][9][10][11][12][13][14][15]. 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Ĥ 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' matrixthe 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 steady-state correlators, can be straightforwardly computed from single-particle quantities. Moreover, the single-body spectrum can be identified with that of a non-Hermitian Hamiltonian [7], 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 [16,17], 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 [18]. The success of the approach induced a variety of different attempts to generalize it [19], culminating in the formulation of the celebrated Bohigas-Giannoni-Schmit conjecture [20] 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

Main results
In the thermodynamic limit N F → ∞, the single-particle properties of random quadratic Liouvillians, specified by Eqs. (1) and (2), 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 [52][53][54][55].
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, Figure 1: Schematic representation of the main results of this paper, the singleparticle 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) Steady-state 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 ratior (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's eigenvalues follow, respectively, RMT and a normal distribution (f). Ther statistics are Poisson in the nonergodic sector and interpolate from GUE to GOE in the ergodic sector as m → 0 (g).
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, signaling 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 behaviors 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, fully mixed, 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 [40,41] 1 . The nonergodic behavior below the transition is, however, not accessible to these fully-random 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 as 2 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 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, 3 defined as

Vectorization and adjoint fermions
To study the eigenvalue problem of the quadratic Liouvillian superoperator, it is convenient to recast 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 N F -dimensional vectors and mixed states and operators by 2 N F × 2 N F matrices. Alternatively, we can see mixed states and 1 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 exactly the same scaling as here, and the steady state exhibits a nonergodic-to-ergodic crossover and a corresponding change in purity as a function of g. 2 We assume a generic case where Liouvillian has no nontrivial Jordan blocks. For the most general treatment see, e.g., Ref. [9]. 3 Though this is certainly true for finite systems, for systems of infinite dimension, the eigenmodes with eigenfrequencies close to the gap can add up and lead to an algebraic relaxation to the steady state. operators as 2 2N F -dimensional vectors over a tensor product of two copies of the fermionic Fock space, while superoperators are represented by 2 2N F × 2 2N F matrices. More explicitly, if |α〉 and |β〉 are basis elements in the fermionic Fock space andÔ 1,2 Fock-space operators, we map Following this procedure, the Liouvillian [Eq.
(1)] is mapped to: 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 [7]. To enforce the canonical anticommutation relations in the vectorized representation also, we define the vector 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., S andS implement the particle-hole in the original and vectorized spaces, respectively. This vectorization scheme is equivalent to third quantization [7]. 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 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 the particle-hole-symmetric and antisymmetric combinations and the non-Hermitian single-particle effective Hamiltonian In Appendix B.1 we show that where the single-particle Liouvillian is given by with Note that the matrix L also satisfies particle-hole symmetry,S L TS = −L.
As shown in Appendix B.2, we can (almost) always find a change of basis that renders the Liouvillian diagonal, 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:

Steady state
From the diagonal form (17), 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 4 , which satisfies the particle-hole symmetry Alternatively, we can define an effective thermal-like HamiltonianΩ as This parametrization is convenient as it automatically takes care of the normalization and positive-definiteness of the steady-state density matrix. Moreover, because the steady state is Gaussian,Ω can be fully characterized by the single-particle matrix Ω, with normalization Z = Tr e which is related to the correlation matrix through the matrix relation: 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 steady-state Lyapunov equation for the correlation matrix: Alternatively, χ can be constructed more efficiently using the right and left eigenvectors of the single-particle matrix L, see Eq. (100) below and Ref. [7]. 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 From it, we define the Majorana operators where C is the Nambu vector defined after Eq. (2), satisfying f † i = f i and the anticommutation relation In this basis, K is transformed into a more suitable matrix, K (U) , where H (U) is an anti-symmetric matrix, and Γ (U) is a symmetric matrix, with N (U) = l (U) T l (U) * and l (U) = l U T .

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 [56]. The simplest way to achieve this is to draw a matrix H int from the Ginibre orthogonal ensemble (GinOE) [57], i.e., sample a real matrix from the probability distribution ∼ exp{−N F Tr H † int H int } and then set 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) [57] 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 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 [44,45], with separate sets of modes decaying at different rates, interspersed by periods of metastability [46].
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 in detail, both numerically and analytically, as a function of g and m.

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 = H − i gΓ . In Refs. [52][53][54][55], 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 [52], supersymmetry [53], diagrammatics [54], or free probability [55], 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 where m = M /(2N F ) and x + i y represents a point in the complex plane 5 . 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 where 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): 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 [52][53][54][55] 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 behaviors 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 single-body 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 [52][53][54][55], see Fig. 2(b), and the single-body 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: 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 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 [52,53,55]: 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 from splitting 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 with an increasingly small absolute real part and the other aligned with the real axis with an increasingly large absolute real part. Physically, the existence of these two distinct regions of eigenvalues means that the group of modes associated with the region with a 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 parts are populated. Eventually, those also fade away and the system reaches 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. (37). 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. (37), 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 behavior 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 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 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. In both panels, the numerical results were obtained by exact diagonalization, for N F = 500 and 400 realizations, and match perfectly the analytical predictions.

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 that follows from the Marchenko-Pastur law [58]. Indeed, we have lim drawn from a real Wishart ensemble with rank 2M , and hence the gap coincides the hard edge of the Marchenko-Pastur distribution. The critical behavior 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 . 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. (36) to determine the gap's scaling behavior 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 behavior 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 [38][39][40]. Because the single-body gap coincides with the manybody 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. 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 on 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. [40]) 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 long-time limit. Because of its Gaussian nature (see Sec. 2), we will base our discussion on 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 single-body operators, which describe the occupation probabilities of different singleparticle 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 singleparticle 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: Because of the particle-hole symmetry of χ, Eq.  (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. (113) of Appendix E, which we rewrite in the Majorana basis by performing the unitary transformation U, Eq. (27): 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 The eigenvalues of Ω, in turn, are: Each l (H) 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 , 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 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 behavior 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. (116) of Appendix E: 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. (116) since γ i +γ j = 0 for some pairs (i, j). More precisely, χ becomes a matrix that acts separately in two subspaces (see Appendix E): subspaceN spanned by the eigenvectors of Γ (U) with a corresponding non-zero eigenvalue, wN 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. Anticipating the results found below, we callN the RMT or ergodic sector and N the Poisson or nonergodic sector.
The component of χ that acts inN (i.e., in the sector with no vanishing γ i ), χNN , is directly obtained from Eq. (116) by replacing w i with wN 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. (126), which can be rewritten as where we used that Γ NN . 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 first note that H (U) N N is anti-symmetric and thus we can diagonalize it with the matrix OU, for an orthogonal matrix O. As a consequence, Eq. (47) can be rewritten as where H is a random variable with average σ = 1 and variance 2σ 2 = 2, by applying the central limit theorem we conclude that Note that, despite hidden in the notation, the variables X and Y are not independent. However, it does not pose any problem as, in the thermodynamic limit, 2M µ=1 γ −1 µ = 2M 〈γ −1 〉 and therefore For ω i , this yields In summary, for m > 1/2, the eigenvalues of Ω are correlated according to RMT, as dictated by Eq. (46). 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. (52) that lim M →∞ χ N N = /2, signaling that the Poisson sector is fully mixed.
In Fig. 5, we plot the spectrum of the single-particle effective Hamiltonian obtained numerically from ED. In Fig. 5  where averaging over the appropriate random ensemble is understood. While 〈γ −1 〉 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(a). In Fig. 5(b) (phase I, m > 1/2), all eigenvalues belong to a single ergodic sector.

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 applies to the purity: 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 which we dub reduced purity, and where P FM = 2 −N F is the purity of the fully-mixed state. The reduced purity, which 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 a 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 behaviors 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 that 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 between 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.
To conclude this subsection, we note that even though the reduced purity completely determines the steady state when it is equal to 0 (χ = /2), it provides only partial information when assuming other values. In particular, since it is not a proper measure of entanglement, it would be interesting to check whether the entanglement of the steady state also obeys a similar nonmonotonic behavior for m < 0.5. As the steady state is mixed, the entanglement entropy is also not a good measure of entanglement and one needs to resort to other more suitable quantities such as the mutual information or the negativity [59,60]. However, this analysis falls out of the scope of this paper and we leave it for future work.

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,r n [61,62]. Let s n be the sequence of differences between consecutive eigenvalues of Ω. Then,r n is defined as [61] r n = min s n+1 s n , s n s n+1 .
This quantity has the advantage of being independent of the local level density and bounded (0 <r n < 1). If the steady state is ergodic, the spectral correlations of Ω coincide with those of 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 [62]: ≈ 0.603, and 〈r〉 Poi = 2 log 2−1 ≈ 0.386. (59) (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 ofr 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(c). However, the results for m < 1/2 show a crossover from GUE statistics at m = 1/2 to Poisson statistics as m → 0. Again, this behavior 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. (52)]. This prediction is confirmed in Fig. 9(b). In contrast, in the  57) and (58). We find perfect agreement with the theoretical predictions in each case. ergodic sector (N ), the eigenvalues of ΩNN 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 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 ΩNN 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 ofN 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 with 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 the expectation of GUE statistics discussed above. As shown in Fig. 9(a) and (b), for finite N F , there is a GUE-to-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 6 . In Appendix F, we show that the spectrum of χNN can be expanded in powers of m and, to lowest nontrivial order (∼ m), it coincides with the spectrum of a matrix from class CI [56], 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: an 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 singlebody 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 behavior in quadratic open quantum systems (see also Ref. [63]) 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 [47][48][49]. 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.

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. [7].
In Ref. [7], it was noted that an orthonormal basis under the Hilbert-Schmidt inner product can be defined on the (vectorized) space of operators by considering the set of all possible vectors 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 where U is defined in Eq. (27). At this point, creation and annihilation operators can be naturally defined on the vectorized space through 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 where It is now easy to see that, after building the matrices we can relateã with D = {d 1 , . . . , The two subspaces that appear in the previous equation correspond to the two parity sectors of Ref. [7]. 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. [7]: where L (γ)

B Spectrum and steady state of fermionic quadratic Liouvillians
In this appendix, we prove the statements made in Sec. 2 about the spectrum and steadystate of quadratic fermionic Liouvillians. In Appendix B.1 we show how to write the singlebody matrix of the Liouvillian using our vectorization scheme, proving Eq. (15). Then, in Appendix B.2 we show how to diagonalize this single-body matrix and, thus, prove Eq. (17). Finally, in Appendix B.3, we connect the steady-state and the single-body spectrum, proving Eq. (26).

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)- (9). First, the von Neumann generator, becomes, after vectorization, 7 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 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: Finally, the dissipative contribution, where we used Tr(H) = Tr(SHS) = − Tr(H) = 0. Note that similarly to the case of H, we can always choose the adjoint fermionic creation and annihilation operators to satisfy particle-hole symmetry, which, in the vectorized space, reads asSA TS = −A, whereS is defined in Eq. (10).

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 2 i jã † i ∆ i jã j } with ∆ = −S∆ TS (note that ∆ can always be chosen to satisfy particlehole symmetry). It follows that which implies that Definingb = R −1ã andb ′T =ã † R, where R = e i∆ , we find: We have reduced the diagonalization of the Liouvillian to the determination of the matrix e −i∆ that diagonalizes L. Note that in generalb ′T ̸ =b † , all we know is that {b i ,b ′ j } = δ i j . Also, the particle-hole symmetry of ∆ leads to the following restriction for the choice of e −i∆ : Any R that diagonalizes L can be made to satisfy Eq. (78) 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, with

JS JS
: Note that the previous equation implies that 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 where X is the solution of Indeed, Recalling that K = H − i gΓ , with Γ a positive semidefinite matrix, all the eigenvalues of K have a non-positive imaginary part, which implies that v (1) 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: and, therefore, Finally, the Liouvillian assumes the form where D L = R −1 LR obeying the particle-hole symmetry: 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. (90) can be further simplified by noting that it follows from Eq. (78) that (1) and b (2) are both vectors of N F annihilation operators and b ′(1) and b ′(2) of the corresponding creation operators. Equation (90) can be finally reduced to 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 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. (93), this implies that or equivalently, 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 L− i 2 β i ̸ = 0 and, hence, Eq. (96) 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, 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 ′ and use the relations along with Eq. (89) to simplify Eq. (98): with R (1) , . . . , R (4) particle-hole blocks of the matrix R 1 , i.e., (4) , and X defined through Eq. (85). 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 antisymmetric 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 behavior of the pth moment of a matrix Q and then infer the resolvent through the relation 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), 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 while the first nontrivial moment, the fourth, is given by: We can associate each Wick contraction with 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 matchings of 2p elements is well-known to be the pth Catalan number, C p , and, hence, we conclude that Tr Q 2p = N F C p in the large-N F limit. Equation (101) 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 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, incidentally, 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: The only unsuppressed contractions (either by symmetrization or nonplanarity) are For the sixth moment, we find 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, 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. (36). After some manipulations, it becomes: 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 on the right-hand side is −m/(4C)g 3 and it is unmatched on the left-hand side, which implies that m = 0. Since the gap at m = 0 trivially vanishes, we must rule out α = −3.
• For α = −1, Eq. (109) becomes Note that this equation is only valid for m < 1/2 as C must be positive.
By inspection of the graphics present in Fig. 2 of the main text, we can assign one of the asymptotic behaviors 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 behavior with g. Since the gap is given by (minus) the intersection point with the largest real part, we conclude that for m < 1/2, Gap ∼ m/4 • 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.
The behaviour of x g (g) in the limit g → 0 can also be extracted from Eq. (109) (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 behavior 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. (85) and (100) in Appendix B, we can write the following equation for the correlation matrix of the steady state: We are interested in determining the solution to Eq. (111) 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. (111). Comparing terms order by order, we obtain that: 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. (112) yields 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 with λ i given by Eq. (113). 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. (111). 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.
However, if γ i + γ j = 0, Eq. (116) 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 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. (116), γ j = 0. Then, because of the property we have just shown, χ (0) i, j = 0. This means that, introducing the projector onto the nullspace of Γ , P N , and onto its orthogonal complement, PN , we have χ This is equivalent to the statement that χ (0) is block diagonal, where χ with (wN ) i an eigenvector of Γ that does not belong to its nullspace, whereas, to compute χ 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: H N N = i ε i (w N ) i (w N ) † i (note that it is always possible since H N N is Hermitian). Plugging this in Eq. (120), 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. (119): We can now solve the second equality in Eq. (119) for χ (1) NN (and its Hermitian conjugate), insert it in the last equality, and use the fact that, since χ 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. (128) as Note that, since the eigenvectors wN i are normalized, we have which means that the norm of wN i under the usual inner product is 1/ γ i . For convenience, we perform the transformation wN i → wN i / γ i , leading to our final expression for χ 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. (133) 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 writẽ where, from the central limit theorem, the spectrum ofΓ (1) , γ (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. (133) and Taylor expanding it up to second order in m = M /(2N F ), we obtain: where we defined We conclude that the diagonalization of χ 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) NN splits into two regions, one close to 0 and the other to 1. J is diagonalized by the matrixŨ where V =Γ V V † belongs to class CI [56,64], which has the same level statistics as the GOE. This finally explains the observation that the level correlations of Eq. (133) approach those of the GOE in the limit m → 0.