Random sparse generators of Markovian evolution and their spectral properties

TL;DR

This paper introduces a model of random sparse matrices as generators of Markovian evolution, analyzing how sparsity affects spectral properties and revealing universal behaviors in ultra-sparse regimes.

Contribution

It develops a new ensemble of sparse Markov generators, studies their spectral characteristics, and uncovers universal spectral statistics for ultra-sparse cases.

Findings

Sparsity closes the large spectral gap in non-sparse generators.

Eigenvalue distribution scales with sparsity parameter $\

,

Abstract

The evolution of a complex multi-state system is often interpreted as a continuous-time Markovian process. To model the relaxation dynamics of such systems, we introduce an ensemble of random sparse matrices which can be used as generators of Markovian evolution. The sparsity is controlled by a parameter $\varphi$, which is the number of non-zero elements per row and column in the generator matrix. Thus, a member of the ensemble is characterized by the Laplacian of a directed regular graph with $D$ vertices (number of system states) and $2 \varphi D$ edges with randomly distributed weights. We study the effects of sparsity on the spectrum of the generator. Sparsity is shown to close the large spectral gap that is characteristic of non-sparse random generators. We show that the first moment of the eigenvalue distribution scales as $\sim \varphi$, while its variance is $\sim \sqrt{\varphi}$. By using extreme value theory, we demonstrate how the shape of the spectral edges is determined by the tails of the corresponding weight distributions, and clarify the behavior of the spectral gap as a function of $D$. Finally, we analyze complex spacing ratio statistics of ultra-sparse generators, $\varphi = \mathrm{const}$, and find that starting already at $\varphi \geqslant 2$, spectra of the generators exhibit universal properties typical of Ginibre's Orthogonal Ensemble.

Tables

Table 1. Table 1: The distributions of the off-diagonal elements M i ​ j subscript 𝑀 𝑖 𝑗 M_{ij} of 𝒦 𝒦 \mathcal{K} (edge weights) and the corresponding distributions of the diagonal elements J l ​ l subscript 𝐽 𝑙 𝑙 J_{ll} of 𝒦 𝒦 \mathcal{K} and the corresponding constants C 𝐶 C and β 𝛽 \beta for the convergence of J l ​ l subscript 𝐽 𝑙 𝑙 J_{ll} to the Weibull distribution Ψ β subscript Ψ 𝛽 \Psi_{\beta} in Eq. ( 27 ). (*) constants obtained by a power-law approximation of the left tail of the gamma distribution.

off-diag. $𝒦=M_{ij}$	${\chi }_k^2$	uniform
diag. $𝒦=J_{ll}$	gamma$(\frac{k\phi }{2},2)$	Irwin-Hall
C	$\frac{2^{\phi }}{\phi !}$*	$\frac{1}{\phi !}$
$\beta$	$\frac{k}{2}\phi$*	$\phi$

Table 2. Table 2: (top) The normalizing parameters c 𝑐 c and d ​ ( D ) 𝑑 𝐷 d(D) for max 1 ≤ l ≤ D ⁡ J l ​ l subscript 1 𝑙 𝐷 subscript 𝐽 𝑙 𝑙 \max_{1\leq l\leq D}J_{ll} to converge to the Gumbel distribution, where J l ​ l subscript 𝐽 𝑙 𝑙 J_{ll} is gamma distributed with shape and rate parameter k 𝑘 k and θ 𝜃 \theta , respectively, see Eq. ( 40 ). (bottom) The relation between the χ 2 superscript 𝜒 2 \chi^{2} and the gamma distribution.

	gamma$(k,\theta )$
$c$	$\frac{1}{\theta }$
$d\left(D\right)$	$\theta \left(\log D+\left(k-1\right)\log \log D-\log \Gamma \left(k\right)\right)$
${\chi }_n^2$	gamma$\left(k=\frac{n}{2}\phi ,\theta =2\right)$

Table 3. Table 3: Mean and angle of spacing ratio distributions obtained with 10 2 superscript 10 2 10^{2} samples of random 10 4 × 10 4 superscript 10 4 superscript 10 4 10^{4}\times 10^{4} matrices rounded to the 4th digit. The matrix ensembles correspond to the ones shown in Figure  7 .

	GinOE	$\phi =1$	$\phi =2$	$\phi =3$
$-⟨\cos \theta ⟩$	0.7379	0.7871	0.7359	0.7372
$⟨r⟩$	0.2347	0.3516	0.2225	0.2284

Equations

$$\frac{d}{dt}P(t)=\mathcal{K}P(t),$$

$$K_{ii}=-\sum _{j\ne i}{K}_{ji}.$$

$$\mathcal{K}=\mathcal{M}-\mathcal{J},$$

$${\mathcal{K}}_{\mathrm{R}M}=-{\mu }_{0}D\cdot 1\mathrm{l}+{\sigma }_0\sqrt{D}(\mathcal{G}-\mathcal{D}),$$

$${\mathcal{K}}^{\prime }=\frac{1}{\sigma \sqrt{D}}(\mathcal{K}+{\mu }_{0}D\cdot 1\mathrm{l})$$

$$\mu (\lambda )=⟨\frac{1}{D}\sum _{i=1}^D{\lambda }_i⟩,$$

$$\mu (\lambda )=⟨\frac{1}{D}\mathrm{tr}(\mathcal{K})⟩=-\phi {\mu }_0,$$

$${\sigma }^2(\mathrm{Re}\lambda )$$

$${\sigma }^2(\mathrm{Im}\lambda )$$

$${\sigma }^2(\lambda )$$

$$={\sigma }^2(\mathrm{Re}\lambda )-{\sigma }^2(\mathrm{Im}\lambda ).$$

$${\sigma }^2(\lambda )$$

$$=\phi \left({\sigma }_0^2+\frac{\phi }{D}{\mu }_0^2-\frac{1}{D}{\sigma }_0^2\right).$$

$$\sigma (\mathrm{Re}\lambda )\gtrsim \sqrt{\phi }.$$

$$\alpha =\frac{|\mu (\mathrm{Re}\lambda )|}{\sigma (\mathrm{Re}\lambda )}.$$

$$\alpha \lesssim \sqrt{\phi }.$$

$$\alpha \approx c_1+c_2\sqrt{\phi },$$

$${\gamma }_{\ast }=\min \{|\mathrm{Re}{\lambda }_i|:\mathrm{Re}{\lambda }_i<0\}.$$

$${\gamma }_{\ast }=-{\lambda }_2=\underset{|v|=1,v\perp 1}{\min }{v}^t(-\mathcal{K})v,$$

$$v_i=\{\begin{array}{ll}\sqrt{1-\frac{1}{D}}& i=l\\ -\frac{1}{\sqrt{D(D-1)}}& i\ne l,\end{array}$$

$${\gamma }_{\ast }\le \underset{1\le l\le D}{\min }{J}_{ll}+O(D^{-1}).$$

$$-{\lambda }_D=\underset{|v|=1}{\max }{v}^t(-\mathcal{K})v,$$

$$-{\lambda }_D\ge \underset{1\le l\le D}{\max }{J}_{ll}.$$

$${\gamma }_{\ast }\approx \underset{1\le l\le D}{\min }{J}_{ll},$$

$$\delta {\gamma }_{\ast }=\frac{|{\gamma }_{\ast }-{\min }_{1\le l\le D}{J}_{ll}|}{{\gamma }_{\ast }}.$$

$$⟨\underset{1\le l\le D}{\min }{J}_{ll}⟩=D\int dxxf(x)(1-F(x){)}^{D-1}.$$

$$⟨{\gamma }_{\ast }⟩\approx D\int dxxf(x)(1-F(x){)}^{D-1}.$$

$$F(x)=C{x}^{\beta }\text{for }0\le x\le C^{1/\beta },$$

$$(DC{)}^{1/\beta }{m}_D\to {\Psi }_{\beta }\text{in law}.$$

$$⟨{\gamma }_{\ast }⟩\approx ⟨m_D⟩\approx \Gamma \left(1+\frac{1}{\phi }\right)(\phi !{)}^{1/\phi }{D}^{-1/\phi }$$

$$⟨{\gamma }_{\ast }⟩\sim D^{-1/\phi }.$$

$$D\approx {\left[\frac{\Gamma \left(1+\frac{1}{\phi }\right)}{c}\right]}^{\phi }\phi !.$$

Full text

I Introduction

Continuous-time Markov chains (CTMCs) [1] provide a popular framework to model dynamics of multi-state systems in diverse fields ranging from physics, chemistry, and biology [2, 3, 4] to economics [5, 6] and game theory [7, 8]. CTMCs are used to model chemical reactions [9, 10, 11, 12, 13, 14, 15], gene regulation processes [16, 17, 18, 19, 20], quantum dynamics (approximated by rate equations) [21, 22, 23, 24, 25], evolutionary game dynamics [8, 26, 27], and many other processes. CTMCs are also the key element of such celebrated models of statistical physics as contact processes [28, 29, 30], zero-range processes [31, 32] and exclusion processes like ASEP [31, 2, 33, 34, 35, 36, 37, 38, 39]. In some fields, CTMCs are known under the names ‘classical Markovian master equations’ or ‘rate equations’.

A continuous-time Markovian evolution in finite discrete space consisting of $D$ states can be specified with a transition rate matrix $\cal{K}$ [1], which is a generator of Markovian evolution. (It is called ‘Kolmogorov operator’ in Ref. [40]). The equation governing the evolution of a probability vector $P(t)$, defined on the state space,

[TABLE]

has the formal solution, $P(t)=\exp(t\mathcal{K})P_{0}$, where $P_{0}=P(0)$ is the initial probability vector. The evolution of $P(t)$ is thus fully determined by the generator $\cal{K}$, especially by its spectral properties. The fact that the operator $\mathcal{M}_{t}=\exp(t\cal{K})$ should map a non-negative vector onto another non-negative vector while preserving $\ell_{1}$-norm, means that $\cal{K}$ satisfies a set of constraints and these constraints have an effect on its spectral properties [1].

In order to model the evolution of a complex system with CTMCs, we would have to first design a specific Kolmogorov operator. Taking into account the large variety of existing models, it would be beneficial to figure out universal properties of $\cal{K}$-generators, i.e., properties that are typical rather than specific to a particular model. The first step in this direction is to define random ensembles of generators. A similar situation arises in the case of unitary time-continuous evolution, where the corresponding generators (quantum Hamiltonians) were explored and classified by using the powerful toolbox of random matrix theory (RMT) [41, 42, 43, 44, 45]. Implementation of this idea resulted in the creation of Quantum Chaos theory [46, 47, 48] which made - and is still making – a strong impact on many-body quantum physics, both theoretical [49] and experimental (see, e.g., Ref.  [50]).

Recently, RMT-based approaches were developed to analyze spectral properties of random generators of open quantum (Lindblad operators) [51, 52, 53, 54] and classical (Kolmogorov operators) [55, 56, 40] Markovian evolution. The considered generators, both quantum and classical, were on purpose constructed in a completely random way - up to the constraints that make them legitimate generators. In the case of Kolmogorov operators, this means that they are represented by dense matrices [55, 40]. It was shown that the spectral density of such operators represents a free sum of a uniform disc and a Gaussian distribution which results in a distinctive spindle-like shape [40], as shown in Figure 2 (a). This density is universal, in the sense that the particular way the operators are sampled does not affect the shape of the spindle (but may affect its position on the real axis and its overall scaling) [40].

In contrast to the random Kolmogorov operators, for most applications and known models, the corresponding $\cal{K}$-generators are represented by sparse matrices. This is a consequence of locality and other topological constraints imposed on the allowed transitions in the state spaces of the models. For many-component or many-particle systems, elements of the generator matrix typically represent changing multiple (or all) components of the system simultaneously, e.g., for an exclusion process, a generic matrix element could represent correlated hopping of many particles. Since such processes are usually absent in physically (biologically, economically,…)-motivated models, most elements of the corresponding $\cal{K}$-matrices are zeros.

Sparsity affects the spectra, $\{\lambda_{i}\}$, $i=1,2,...,D$, of the corresponding generators. Most noticeably, the spectral gap, i.e., the distance between $\lambda_{1}=0$ and the eigenvalue closest to it, $\gamma_{*}=\min\{|\operatorname{Re}\lambda_{i}|\}$, does not grow with the increase of number of states $D$. This is in sharp contrast to the case of dense random generators  [51, 53, 54] . In Figure 1 we contrast the case of dense random stochastic matrices (a) with various model systems described by sparse Markov generators (b-f).

The large gap of dense random generators implies that even the slowest decaying mode of a generic initial probability vector converges rapidly to the equilibrium state, the relaxation time (inverse of the spectral gap) decreasing inversely in the state space size $D$. In contrast, physical generators of CTMCs in general have spectral gaps and relaxation times functionally depending on $D$ very differently than (anti-)linearly; see, e.g. [57, 58, 59, 60] for the example of the exclusion process. This difference of behavior suggests that it is more suitable to model physical CTMC generators by sparse rather than dense random matrices.

Our motivation is to refine the RMT approach to random Kolmogorov operators by including sparsity which is characteristic to physically relevant $\cal{K}$-generators. We specify an ensemble of random matrices of fixed sparsity $\varphi$ as an ensemble of negative combinatorial Laplacians of random regular directed graphs. The sparsity is controlled by the vertex degree $\varphi$ which is equal to the number of non-zero elements per row and column of the generator matrix. In graph terms, this means that each vertex has exactly $\varphi$ incoming and $\varphi$ outgoing edges. The non-zero elements (edge weights) are taken to be random, positive, independent, and identically distributed (iid).

A similar setup was studied in Ref. [56], where an ensemble of oriented Erdős-Rényi graphs [62], parameterized with edge probability distribution $p(D)$, was used. The vertex distribution in this case is binomial-distributed [62], and not constant as in our case. However, one might expect similar behavior in the $D\to\infty$ limit with the correspondence $p(D)=\varphi/D$. The authors of Ref. [56] considered the regime $Dp(D)\gg(\log D)^{6}$, which they found to have the same universal properties as in the non-sparse case. In this work, we consider sparsity beyond this limit, including specifically $\varphi\sim D^{0}$ (vertex degree not growing with $D$) and $\varphi\sim\log D$.

In this paper, we investigate the dependence of spectral properties of the sparse Kolmogorov operators on sparsity parameter $\varphi$, number of states (dimension of the state space) $D$, and on the edge weight distribution, i.e., on the distribution of the nonzero elements of matrix $\mathcal{K}$. Explicit results are mostly derived for $\chi_{2}^{2}$ and uniform weight distributions; however, these results can be adapted to other weight distributions.

We consider both the bulk of the spectral distribution and its edges.

As for the bulk, we focus on its position $\mu$ (the mean of the corresponding eigenvalue distribution) and its variances along the real and imaginary axes (standard deviations of the distribution of the real and imaginary eigenvalue parts, respectively). The first variance estimates the spread of the relaxation rates, while the second one gives the timescales of the oscillations during the relaxation.

As for the edges, we address the spectral gap and the extent of the spectrum along the real axis (the real part of the eigenvalue with largest absolute real part). These determine respectively the slowest and fastest time scales of relaxation to the steady state. The spectral gap is of physical interest for multi-state Markov processes, see, e.g., [57, 58, 59, 60, 63, 64, 65, 66, 67] for the ASEP and [68, 69] for the contact process. The horizontal extent, in addition to its interpretation as the fastest timescale for CTMCs, is also relevant in the graph theory interpretation, e.g., to quantify the computational complexity of the community detection problem [70, 71] and the max-cut [72, 73] problems.

We demonstrate that the position and variance of the spectral bulk of sparse Kolmogorov operators scale as $\sim\varphi$ and $\sim\sqrt{\varphi}$, respectively. These characteristics do not depend on $D$ but on the first and second moments of the weight distribution. The dependence of the spectral edges on the weight distribution is less straightforward and highly non-universal. In particular, we show that, in the regime of high sparsity, $\varphi\ll D$, the spectral gap (horizontal extent) depends only on the left (right) tail of the weight distribution. We evaluate the dependence of the spectral gap on $\varphi$ and $D$ for weight distributions with exponential and power-law tails.

We consider the cases of $\chi_{2}^{2}$ and uniform weight distributions in detail. For these distributions, we find that the spectral gap and the horizontal extent of the spectrum can be approximated by the largest and smallest diagonal entry of the generator matrix, respectively. Using the conjecture that this correspondence holds in general, we use extreme value theory (EVT) [74, 75] to analytically derive dependencies of the spectral edges on $\varphi$ and $D$. In particular, we infer that the distributions of spectral edges only depend on the tails of the weight distributions.

Finally, we analyze correlations of the eigenvalues of sparse Kolmogorov operators. We show that, for $\varphi\geq 2$, the complex spacing ratio distributions [61] of the spectral bulks follow the distribution typical to Ginibre’s Orthogonal Ensemble.

The paper is organized as follows. In Section II we introduce an ensemble of sparse random Kolmogorov operators. We analyze the bulk of the spectral distributions of the ensembles in Section III. In Sections IV and V we address the spectral gap and the horizontal extent of the spectrum, respectively. A discussion on correlations between eigenvalues in terms of the complex spacing ratio follows in Section VI. We conclude with a summary of our results in Section VII. Appendices contain information on the models whose spectra are presented in Figure 1, details of the sampling of sparse random Kolmogorov operators, and the details of analytical derivations.

II Random sparse Kolmogorov operators

In this section, we first recall some basic properties of Kolmogorov operators and review the case of full random $\mathcal{K}$-matrices. We then define an ensemble of random sparse operators. In what follows, matrices will be referred to by calligraphic letters (e.g., $\mathcal{K}$) while their elements will be denoted by non-calligraphic letters (e.g., $K_{ij}$).

II.1 Basic information

In order to be qualified as a Kolmogorov operator, a matrix $\cal{K}$ has to fulfill two conditions, (i) all its off-diagonal elements have to be real and non-negative, $K_{ij}\geq 0$, $i\neq j$, and (ii) the sum over every column should be zero. The latter is fulfilled by setting all the diagonal elements as

[TABLE]

The first condition guarantees the preservation of the non-negativity of a vector during the evolution induced by Eq. (1), while the second one guarantees the preservation of the $\ell_{1}$-norm of the vector.

The spectrum of $\cal{K}$ is in general complex. Since $\cal{K}$ maps real vectors onto real vectors, the spectrum is invariant under complex conjugation, so all complex eigenvalues come in conjugated pairs. The spectrum contains at least one eigenvalue $\lambda_{1}=0$ with right eigenvector corresponding to the steady state. By virtue of the Perron-Frobenius theorem [76, 77, 78], the components of the steady state vector can be chosen to be non-negative, which makes it, after normalization, a probability vector. Any Kolmogorov operator can be represented in terms of a real non-negative matrix, $\cal{M}$, $M_{ij}\geq 0$,

[TABLE]

where elements of the diagonal matrix $\cal{J}$ are $J_{jj}=\sum_{i}M_{ij}$.

We now briefly review the case of dense (non-sparse) random Kolmogorov operators [55, 56, 40]. Elements $M_{ij}>0$ are i.i.d. sampled from a distribution with density $p(x)$ and first two moments $\mu_{0}=\int xp(x)dx$ and $\sigma^{2}_{0}=\int(x-\mu_{0})^{2}p(x)dx$. The particular choice of distribution does not play an essential role (provided that it is not very pathological). For example, we could sample a matrix $\cal{Z}$ from Ginibre’s Unitary Ensemble (GinUE) and then square its elements, $M_{ij}=|Z_{ij}|^{2}$ [40]). The matrix $\cal{M}$ is then full in the sense that, with probability $1$, all its elements are different from zero.

The elements of the matrix $\mathcal{M}$ are i.i.d., thus, in the asymptotic limit, its spectral density is a uniform disk of radius $\sqrt{D}\sigma_{0}$, with the center at [math]. In the dense limit, the elements of $\cal{J}$ are sums of $D$ independent random variables, so its elements can be approximated with Gaussian-distributed random variables with mean $D\mu_{0}$ and variance $D\sigma_{0}^{2}$.

Following the RMT approach [40], the Kolmogorov operator in Eq. (2) can be modelled as

[TABLE]

where $\cal{G}$ is a member of Ginibre’s Orthogonal Ensemble (GinOE) and $\mathcal{D}$ is a diagonal matrix. Elements of $\mathcal{G}$ and $\mathcal{D}$ are sampled from the normal distribution of zero mean and unit variance. Here $\sigma_{0}\sqrt{D}\cdot{\cal G}$ models $\cal{M}$ while $\cal{J}$ is approximated as $\mu_{0}D\cdot{\mathchoice{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm 1\mskip-5.0mul}}+\sigma_{0}\sqrt{D}\cdot{\cal D}$.

The spectral density of the non-trivial part, ${\cal K}^{\prime}={\cal G}+{\cal D}$, is a free convolution of a disk and a Gaussian distribution along the real axis, which results in a spindle-like shape. Figure 2 (a) presents both the spectrum of a single random dense Kolmogorov operator and histogram obtained with $100$ samples.

Alternatively, we can state that the spectral density of the rescaled generator

[TABLE]

is expected, in the asymptotic limit, to be the $D$-independent spindle (”an additive Gaussian deformation of the circular law”, according to Ref. [56]).

The spectrum of the random non-sparse generator has a large gap which scales as $D$, as seen in Figures 1 and 2. We will see that this feature is strongly affected when we introduce sparsity.

II.2 Ensemble of sparse random Kolmogorov operators as a set of oriented graphs

The operator ${\cal K}$ described in the introduction can be considered as the negative Laplacian of a random directed graph with positive, iid edge weights, without self-loops, and with fixed vertex degree equal to $\varphi$.

For example, the graph corresponding to the ${\cal K}$ generator of a process of a particle hopping on a $d$-dimensional hypercubic lattice with periodic boundary conditions and random hopping rates is a particular (to the nearest-neighbor connections) realization of the ensemble with $\varphi=2d$. Figure 1 (d) shows an example spectrum for $d=2$.

The regularity of the graphs ensures that, with probability $1-O(D^{-\varphi-1})$, they are strongly connected as long as $\varphi\geq 2$ [79]. Strong connectivity is a good feature because it means that the matrix ${\cal K}$ is not of block-diagonal structure and the state space is not partitioned into disconnected subsets. As there is only one strongly connected component, there is only one absorbing component. This implies that the multiplicity of the zero eigenvalue is one and the steady state is unique. Finally, every state in the state space is reachable from every other state. The steady state, therefore, has all states populated.

Some physical models motivating this study presented in Figure 1 are - except for the contact process - all strongly connected. The contact process is only effectively strongly connected. It has two strongly connected components, where one is a single vertex and the other includes the remaining $D-1$ vertices. The giant component is the only absorbing component and consequently, the steady state is unique.

We will focus on sparse generators with $\varphi\geq 2$ and will discuss $\varphi=1$ in Section VI.

The physical models presented in Figure 1 motivate us to focus on two types of dependencies of $\varphi$ on the matrix size $D$, namely $\varphi=\mathrm{const}$ and $\varphi\sim\log D$. For generators of single particle hopping models - an example is shown in Figure 1 (d) - the average number of non-zero elements per column and row is constant and independent of $D$. It increases logarithmically with $D$ in many-body hopping models such as the ASEP or the contact process, Figure 1 (b), (c), and (e). There is no simple dependence of $\varphi$ on $D$ in the gene transcription model, Figure 1 (f), as the matrix size $D$ is controlled by multiple parameters, see Appendix A.

What can we say about spectral densities of the ultra-sparse $\mathcal{K}$-generators, with $\varphi=\mathrm{const}$? A ’naive’ adjustment of the RMT approach, which consists in describing the elements of a sparse matrix $\cal{M}$ with probability density function $\widetilde{p}(x)=(1-\frac{\varphi}{D})\delta(x)+\frac{\varphi}{D}p(x)$, re-scaling the mean and variance accordingly, and then using the RMT model, Eq. (4), would not work here for two reasons. First, the spectral densities of such sparse matrices cannot be approximated with members of ’dense’ RMT ensembles. Second, the Central Limit Theorem no longer applies and the entries of matrix $\mathcal{J}$ cannot be approximated with normal random variables (the entries become distribution-specific).

III Position and width of the bulk of the spectrum

In this section, we analyze the dependence of the position and horizontal width of the bulk of the spectrum on the sparsity parameter $\varphi$ and the matrix dimension $D$. We first provide (Subsections III.1 and III.2) expressions and bounds for the position and the width, characterized respectively by the mean $\mu(\lambda)$ of all eigenvalues and the standard deviation $\sigma(\operatorname{Re}\lambda)$ of the real parts of the eigenvalues. These results are expressed in terms of the mean and standard deviation of the weight distribution (distribution of non-zero elements of the Kolmogorov operator ${\cal K}$), denoted by $\mu_{0}$ and $\sigma_{0}$ respectively.

Since the most prominent effect of sparsity is to reduce the parametrically large gap seen in the full random case, it is instructive to analyze the ratio $\alpha=|\mu(\lambda)|/\sigma(\operatorname{Re}\lambda)$. This quantity provides insight into the distance of the bulk of the spectrum from the origin, relative to the size of the bulk. Subsection III.3 is devoted to an analysis of the ratio $\alpha$.

Numerical results presented in this section are obtained by sampling edge weights from the $\chi_{2}^{2}$ and the standard uniform distribution.

The spectrum of dense generators ($\varphi=D-1$) consists of two distinct parts - an eigenvalue $\lambda_{1}=0$ and the rest of the eigenvalues forming the spectral bulk away from the imaginary axis, as shown in Figure 1 (a) and Figure 2 (a). In contrast, the bulk of the spectrum is much closer to the imaginary axis for $\varphi\ll D$, as seen in Figure 2 for (b) $\varphi=\sqrt{D}$, (c) for $\varphi=\log D$ and (d) for $\varphi=3$. For $\varphi=\sqrt{D}$, the bulk of the spectrum is visibly separated from the zero, as in the dense case. In fact, the spectral boundary is given by the same spindle (properly rescaled). Whether the spectral distribution is separated from zero for $\varphi=\log D$ and $\phi=3$ is difficult to say with certainty from the available numerical data ($D\approx 8000$).

III.1 Position

The position of the spectral bulk of ${\cal K}$ can be identified with the mean $\mu(\lambda)$ of eigenvalues $\lambda_{i}$,

[TABLE]

where the average $\langle\dots\rangle$ is taken over the ensemble of random Kolmogorov operators described in Section II. Because the eigenvalues are either real or come in complex conjugate pairs, the mean of the spectral bulk is real, $\mu(\lambda)=\mu(\operatorname{Re}\lambda)$.

A simple calculation, presented in Appendix C, shows that $\mu(\lambda)$ can be expressed as

[TABLE]

The averaging $\langle\dots\rangle$ over the matrix ensemble in Eq. (6) and Eq. (7) is, in principle, not needed since typicality is expected, i.e., for large enough $D$, a single sample will display all the spectral features of the ensemble. This is because the quantity $\frac{1}{D}\operatorname{tr}(\mathcal{K})$ is concentrated around its average $\left\langle\frac{1}{D}\operatorname{tr}(\mathcal{K})\right\rangle$ for increasing $D$, as shown in Appendix C.

For the four different dependencies of $\varphi$ on $D$ shown in Figure 2, Eq. (7) implies the following: For $\varphi=\mathrm{const}$, the mean is independent of the matrix size $D$. For $\varphi=\log D$ ($\varphi=\sqrt{D}$) the mean decreases logarithmically with $D$ (as $\sim\sqrt{D}$) and for $\varphi=D$ the mean decreases linearly with $D$ as is expected for the dense generators [55].

In Figure 2, the location $\mu(\lambda)$ of generator matrices $\cal{K}$ is indicated with a red dot in each panel. The real part of the dot resides in the bulk of the spectrum for every dependence of $\varphi$ on $D$ shown in Figure  2.

III.2 Horizontal width

In Section III.1 we investigated where the bulk of the spectrum is located in the complex plane. We now analyze the width of the distribution. We are especially interested in the horizontal width.

We characterize the width of the bulk spectrum, both in the real and imaginary directions, $\operatorname{Re}\lambda$ and $\operatorname{Im}\lambda$, using the estimated variances

[TABLE]

where we used the fact that $\sum_{j=1}^{D}\lambda_{j}$ is real.

Because the eigenvalues appear in complex conjugate pairs, $\sigma^{2}(\operatorname{Re}\lambda)$ and $\sigma^{2}(\operatorname{Im}\lambda)$ are related to the estimated complex pseudo-variance via

[TABLE]

The estimated pseudo variance lower bounds the estimated variance of the real parts of the eigenvalues, $\sigma^{2}(\lambda)\leq\sigma^{2}(\operatorname{Re}\lambda)$.

The complex pseudo variance can be analytically calculated for the ensemble of random generator matrices as

[TABLE]

Details of the calculation are provided in Appendix C. The bound of the estimated real variance by the pseudo variance together with Eq. (III.2) leads to the asymptotic lower bound of $\sigma(\operatorname{Re}\lambda)$ in terms of the sparsity parameter $\varphi$. As $1\leq\varphi\leq D$, the estimated horizontal width of the bulk spectrum cannot grow asymptotically slower than $\sqrt{\varphi}$,

[TABLE]

Numerically, we find that the bound in Eq. (12) is asymptotically sharp for $\varphi\ll D$, as shown in Figure 3 through the ratio $\alpha$ of mean $\mu(\operatorname{Re}\lambda)$ and width $\sigma(\operatorname{Re}\lambda)$. The collapse of the data points in Figure 3 (c) implies that $\sigma(\operatorname{Re}\lambda)\sim\sqrt{\varphi}$.

III.3 Ratio of mean and horizontal width

In this section, we combine the information of the location of the spectrum given by Eq. (6) and the horizontal width of the bulk given by Eq. (8) into the ratio

[TABLE]

This quantifies how close the bulk spectrum is, relative to its size, to the stationary value $\lambda_{1}=0$. i.e., to the imaginary axis. For $\alpha=O(1)$ the estimated width of the bulk is of the same order as the estimated mean, thus the spectrum is located close to 0. For $\alpha\gg 1$ the estimated mean is much bigger than the horizontal width of the bulk and the bulk of the spectrum is far away from [math].

The analytical result for the estimated mean of the spectrum, Eq. (7), together with the asymptotic bound on the standard deviation of the real parts of the spectrum, Eq. (III.2), imply the following asymptotic bound on $\alpha$

[TABLE]

Numerically, we observe that the bound in Eq. (14) is asymptotically tight for $\varphi\ll D$, i.e.

[TABLE]

for constants $c_{1}$ and $c_{2}$. Since $\mu(\lambda)$ scales linearly with $\varphi$, this behavior is consistent with $\sigma(\operatorname{Re}\lambda)\sim\sqrt{\varphi}$, stated previously. The constants are found to be $c_{1}\approx 0.15$ ($\approx 0.1$) and $c_{2}\approx 0.84$ ($\approx 1.3$) for the $\chi_{2}^{2}$ (uniform) distribution.

Numerical results for $\alpha$ are summarized in Figure 3. For each combination of $\varphi$ and $D$, $\alpha$ is averaged over $n$ samples of random generators such that $nD=50^{\prime}000$. The weight distribution is the $\chi_{2}^{2}$ distribution in (a) and in the lower part of (c), and is the uniform distribution in $[0,1]$ in (b) and in the upper part of Figure 3 (c). We have found that these results are qualitatively the same for exponentially distributed edge weights.

In Figure 3 (a,b) we show the value of $\alpha$ as a function of $D$ and $\varphi$. On the $x$-axis $D$ varies in steps of $10^{3}$ between $10^{3}$ and $10^{4}$. We observe that $\alpha$ increases with $\varphi$ and is independent of $D$, as predicted by Eq. (15). In Figure 3 (c) we show $\alpha$ as a function of $\varphi$ for different dependencies of $\varphi$ on $D$. In all the cases, values of $\alpha$ collapse onto the black solid line given by Eq. (15).

For $\varphi\sim D$, the ratio $\alpha$ scales as $\sim\sqrt{D}$, thus recovering the parametrically large gap in the non-sparse case. For constant $\varphi$, the location of the bulk relative to its size is constant and independent of $D$, i.e, if measured relative to the size of the bulk, the bulk does not move away from the imaginary axis with increasing $D$. We have thus quantified how sparsity cures one of the less physical aspects of the non-sparse random model of Markov generators.

IV Spectral gap

In this and the following section, we will consider the spectral edges, namely, the locations of the eigenvalues nearest and farthest from the imaginary axis. In this section we will investigate the spectral gap $\gamma_{*}$ of $\mathcal{K}$,

[TABLE]

The spectral gap $\gamma_{*}$ is asymptotically, approximately bounded by the right extent of the bulk $|\mu(\lambda)|-\sigma(\lambda)$, which depends on $\varphi$ as $\sim\varphi-\sqrt{\varphi}\sim\varphi$. So for constant $\varphi$, the spectral gap is bounded from above, while for $\varphi$ increasing with $D$ the spectral gap can increase with $D$.

Here the edge weights are distributed according to the $\chi_{2}^{2}$ and the standard uniform distributions. We first demonstrate that, for $\varphi=\mathrm{const}$, the average spectral gap $\langle\gamma_{*}\rangle$ decreases as $D^{-1/\varphi}$, while $\langle\gamma_{*}\rangle$ is constant if $\varphi$ increases logarithmically with $D$. We then show that the spectral gap is well approximated by the smallest (in magnitude) diagonal term of $\mathcal{J}$ ($\mathcal{K}$) and use the theory of extreme values to underpin the numerical observations. The results are then generalized to weight distributions with power-law left tails in that for constant $\varphi$ the average spectral gap decreases as a power-law in $D$ and the crossover from decreasing to increasing $\langle\gamma_{*}\rangle$ happens when $\varphi\sim\log D$.

IV.1 Numerical results

In Figure 4 we show the average spectral gap $\langle\gamma_{*}\rangle$ for edge weights distributed as $\chi_{2}^{2}$ (a-c) and according to the standard uniform distribution (d-f). For every combination of $\varphi$ and $D$, the average of the spectral gap is estimated with $100$ samples. In Figure 4 (a) and (d) we show $\langle\gamma_{*}\rangle$ as a function of $D$ for different dependencies of $\varphi$ on $D$. The average spectral gaps for constant $\varphi=3,5,8,13$ (presented with colored circles) clearly follow a power-law scaling with $D$.

In Figure 4 (b) and (e) we show the average spectral gap $\langle\gamma_{*}\rangle$ as a function of $\varphi$ and $D$. The black dashed lines are contour lines of constant $\langle\gamma_{*}\rangle$. They are near straight lines, showing that for a logarithmic increase of $\varphi$ in $D$ the spectral gap is constant.

We show the average spectral gap $\langle\gamma_{*}\rangle$ as a function of $D$ for $\varphi=\frac{4}{5}\log D+8$ in Figure 4 (a) and $\varphi=\frac{7}{10}\log D+8$ in (d) as black diamonds. These dependencies of $\varphi$ on $D$ agree well with the top dashed contour lines in (b) and (e), respectively. The average spectral gap of $\varphi$ depending logarithmically on $D$ is constant in Figure 4 (a) and (d).

IV.2 Gap $\approx$ minimum of $\mathcal{J}$

Let us assume for a moment that the generator matrix $\mathcal{K}$ is hermitian with eigenvalues $\lambda_{D}\leq\dots\leq\lambda_{2}<\lambda_{1}=0$. Then $\mathbb{1}=(1,\dots,1)^{t}$ is the eigenvector with eigenvalue 0 and all other eigenvectors are orthogonal to it. By the Courant-Fischer theorem [80]

[TABLE]

where the minimum runs over all vectors $v\in\mathbb{R}^{D}$, which have Euclidean norm $|v|=1$ and are perpendicular to $\mathbb{1}$. Choosing $1\leq l\leq D$ arbitrary and $v$ as (see Appendix D for more details)

[TABLE]

a simple calculation shows that (at least for $\varphi\ll D$)

[TABLE]

Similarly, by using the Courant-Fisher theorem, for the eigenvalue with largest magnitude $\lambda_{D}$ we find

[TABLE]

and with $v$ as the $l$-th vector of the standard basis of $\mathbb{R}^{D}$

[TABLE]

Under some mild conditions on random weights $K_{ij}$, a result from Ref. [71] shows that the inequality Eq. (21) becomes an equality in the large $D$ limit with probability approaching 1. Motivated by this observation and the bound from Eq. (19), we expect a similar asymptotic tightness for Eq. (19). However, it is an open question whether the result from Ref. [71] applies to the bound of the spectral gap, Eq. (19). Further, the proof presented in Ref. [71] makes use of the Central Limit Theorem for the diagonal elements $J_{ll}$ of $\mathcal{J}$, and so the corresponding result does not apply to the case of constant or logarithmically increasing (with $D$) sparsity parameter $\varphi$.

Nevertheless, the above arguments allow us to conjecture that in the limit of large $D$ and $\varphi\ll D$ the following

[TABLE]

holds for general, non-hermitian random generator matrices $\mathcal{K}$, with iid and non-exotic weight distributions. We support our conjecture with numerical data presented in Figures 5 (a) and (b). We quantify the approximation in Eq. (22) by the relative error between the spectral gap $\gamma_{*}$ and the minimum $\min_{1\leq l\leq D}J_{ll}$ of the diagonal of $\mathcal{J}$,

[TABLE]

Figure 5 shows $\langle\delta\gamma_{*}\rangle$ as a function of $\varphi$ and $D$ for the $\chi_{2}^{2}$ distribution and the standard uniform distribution. The average relative error is at least two orders of magnitude smaller than the average spectral gap shown in Figure 4 (b) and (e). For increasing $D$, the approximation in Eq. (22) improves. Thus, the approximation in Eq. (22), works well in the case $\varphi\ll D$.

IV.3 Extreme value theory

The distribution of the right-hand side of Eq. (22) can be tackled with the theory of extreme values. As all non-zero entries of $\mathcal{M}$ (edge weights) are identically and independently distributed random variables, so are the diagonal entries of $\mathcal{J}$. Let the cumulative distribution function (cdf) of the diagonal entries $J_{ll}$ of $\mathcal{J}$ be denoted by $F$ and its probability density function by $f(x)=\frac{d}{dx}F(x)$. If the edge weights are distributed according to a $\chi^{2}$ distribution (or any gamma distribution) the cdf $F$ of $J_{ll}$ is a gamma distribution function. If the edge weights are uniformly distributed, $F$ is an Irwin-Hall distribution function, see Table 1. The expected value of $\min_{1\leq l\leq D}J_{ll}$ is given in terms of $F$ (and $f$) by

[TABLE]

Eq. (22) and Eq. (24) imply that

[TABLE]

We demonstrate the validity of Eq. (25) with Figure 4 (a), where the solid lines, given by Eq. (25), perfectly match numerically sampled average spectral gap $\langle\gamma_{*}\rangle$. In the next section, we will use the theory of extreme values to handle the integral in Eq. (25).

IV.3.1 Power-law tail distributions

Let us consider first the case $\varphi=\mathrm{const}$ and increasing $D$. By the Fisher-Tippet-Gnedenko (or ’extreme value’) theorem [75], $\min_{1\leq l\leq D}J_{ll}$ converges in law, under some mild assumptions on the distribution of $J_{ll}$ and properly renormalization, to the Weibull distribution. The Weibull cumulative distribution function is given by $\Psi_{\beta}(x)=e^{-x^{\beta}}$, where $\beta>0$ and the support is on the positive real line.

For distributions of $J_{ll}$ with power-law left tail, the renormalization of $\min_{1\leq l\leq D}J_{ll}$ for convergence to the Weibull distribution is well known, see e.g. Theorem 3.3.2, page 137 in Ref. [75]. We use a version modified to our case. Let a positive random variable $X$ have cdf $F$ with $\beta$-power left tail, i.e.

[TABLE]

where $C>0$ is a constant. Further, let $m_{D}=\min_{1\leq l\leq D}X_{l}$, where the $X_{l}$ are iid copies of $X$. Then

[TABLE]

The Irwin-Hall distribution has a left power-law tail given by $F(x)=\frac{x^{\varphi}}{\varphi!}$ for $0\leq x\leq 1$. The constants for the Irwin-Hall distribution are listed in table 1.

We assume that the convergence in Eq. (27) is not only in distribution but that the renormalized moments of $m_{D}$ converge as well. If the convergence of the moments is sufficiently fast, then Eq. (27) together with Eq. (22) imply

[TABLE]

when the weight distribution (distribution of non-zero off-diagonal elements of $\mathcal{K}$) is such that the diagonal of $\mathcal{J}$ has a power-law left tail and the coefficients $C$ and $\beta$ are given by $C=1/\varphi!$ and $\beta=\varphi$.

Finally, we consider the case that the weight distribution is uniform. We observe that the approximation in Eq. (28) works very well in this case. The solid lines in Figure 4 (d) are given by the right-hand side of Eq. (28) and they match the numerically calculated average spectral gap.

Eq. (28) implies that, for constant $\varphi=\mathrm{const}$ and increasing $D$, the average spectral gap decreases as

[TABLE]

In Figure 4 (f) we show that the numerically retrieved power-law exponents of the average spectral gap, Figure 4 (d), match the scaling in Eq. (29).

We find that the large deviation result is not only valid for constant $\varphi$ and increasing $D$ but also for $\varphi$ increasing logarithmically with $D$; see Figure 4 (d). This allows us to estimate the crossover from decreasing to increasing spectral gap. Let $c$ denote a constant and let $\langle\gamma_{*}\rangle=c$. Then by Eq. (28)

[TABLE]

In Figure 4 (e) the contour lines of constant average spectral gap $c$ perfectly line up with the functional dependence of $D$ on $\varphi$ through Eq. (30) shown as white dots.

To find $\varphi$ as a function of $D$ such that the average spectral gap is constant, we assume that $\varphi$ is reasonably large and approximate $\Gamma\left(1+\frac{1}{\varphi}\right)\approx 1$ and by Stirling’s formula $(\varphi!)^{1/\varphi}\approx\frac{\varphi}{e}$. Denoting $y=\log\frac{\varphi}{ce}$ and rearranging Eq. (30) gives us

[TABLE]

which can be inverted by the Lambert $W$ function. Resubstituting $\varphi=cee^{y}$ we arrive at

[TABLE]

which for $\log D\geq ce^{2}$ behaves as [81]

[TABLE]

where $\eta(D)\to 0$ slowly, as $\eta(D)\sim(\log\log D)^{-1}$. So in the limit $1\ll\varphi\ll D$ the crossover from decreasing to increasing spectral gap happens at $\varphi\sim\log D$ with corrections of the order $\log\log D$. This confirms our numerical observations that the average spectral gap $\langle\gamma_{*}\rangle$ appears to be constant for $\varphi\sim\log D$ in the range of matrix sizes $D$ we considered.

IV.3.2 Approximate power-law distributions

If the weight distribution is a $\chi^{2}$ or exponential distribution, the diagonal elements of $\mathcal{J}$ are distributed according to Gamma distribution, see table 1. The left tail of the Gamma distribution only follows approximately a power-law. Approximating the left tail by a Taylor expansion, we obtain constants $C$ and $\beta$ presented in Table 1. Especially, for the $\chi_{2}^{2}$ distribution, we presented so far in the main text the power-law approximation of the gamma distribution and the large deviation result in the previous subsection suggest that the average spectral gap $\langle\gamma_{*}\rangle$ decreases for constant $\varphi$ and increasing $D$ as a power in $D$ with exponent given $-1/\varphi$, see Eq. (29).

In Figure 4 (c) we present the numerically calculated exponents of the power-law decrease of $\langle\gamma_{*}\rangle$, for $\chi_{2}^{2}$ weight distribution, with $D$ and compare it to the prediction $-1/\varphi$. We find excellent agreement for small $\varphi\leq 5$. For larger $\varphi$ the deviation between the numerical exponent and $-1/\varphi$ is visible but the agreement is still good.

A quantitative comparison between the numerically calculated spectral gap $\langle\gamma_{*}\rangle$ and the EVT prediction by a power-law approximation of the left tail of the gamma distribution resulted in poor agreement. As the expected minimum value of the diagonal of $\mathcal{K}$ perfectly agrees with $\langle\gamma_{*}\rangle$, we attribute the disagreement to the power-law approximation of the left tail and slow convergence of Eq. (27) for diagonal elements of $\mathcal{J}$ distributed according to the gamma distribution.

IV.4 Summary

We presented numerical and analytical arguments that, for the weight distributions considered, the average spectral gap decreases as a power-law for constant $\varphi$ and increasing $D$ with exponent given (approximately) by $-1/\varphi$. The crossover between decreasing and increasing spectral gap happens at $\varphi\sim\log D$, with $\log\log D$ corrections, for uniform weight distribution. For $\chi_{2}^{2}$ distributed edge weights the crossover was observed at $\varphi\sim\log D$. If $\varphi$ increases with $D$ faster than $\log D$ then the average spectral gap increases.

The presented results generalize. Let us assume that the spectral gap is well approximated by the smallest (in magnitude) diagonal of $\mathcal{J}$, at least in the regime of large $D$ and $\varphi\ll D$. Then, after appropriate renormalization, the distribution of the spectral gap is given by the limiting extreme value distribution of the diagonal elements of $\mathcal{J}$. Thus the classification of functional dependencies of the spectral gap on $\varphi$ and $D$ with respect to weight distributions reduces to the classification of limiting extreme value distributions and renormalizations. Extensive research has been conducted on the latter and the renormalizations of a lot of common distributions are well known [74, 75]. Thus the presented approach allows the calculation of the distribution of the spectral gap for broad classes of weight distributions.

V Horizontal extent (largest absolute real part)

In this section we investigate the horizontal extent $\tilde{\gamma}$ of the spectrum given by the eigenvalue with largest absolute real part

[TABLE]

We focus on the averaged horizontal extent $\langle\tilde{\gamma}\rangle$. We show that for $\varphi\sim\log D$ the average horizontal extent increases logarithmically with $D$ for $\chi_{2}^{2}$ or uniformly distributed edge weights. For constant $\varphi$ and increasing $D$ the dependence of $\langle\tilde{\gamma}\rangle$ is qualitatively very different for the two distributions. For the $\chi_{2}^{2}$ distribution $\langle\tilde{\gamma}\rangle$ increases logarithmically, while for the uniform distribution, the average horizontal extent converges to $\varphi$ as a power-law in $D$. Ultimately, this is because the support of the uniform distribution is bounded, while the right tail of the $\chi_{2}^{2}$ distribution extends to infinity.

The structure of this section follows closely the one from Section IV. We first present numerical results demonstrating the above statements. We then argue that the horizontal extent is given by the largest, in magnitude, diagonal element of $\mathcal{K}$ and invoke again EVT to analytically underpin the functional dependencies of $\langle\tilde{\gamma}\rangle$ on $\varphi$ and $D$.

V.1 Numerical results

In Figure 6 we show the average horizontal extent as a function of $D$ for constant $\varphi$ and $\varphi\sim\log D$ for edge weights distributed according to a $\chi_{2}^{2}$ (a,c) and the standard uniform distribution (b,d). In (a) the dependence of $\langle\tilde{\gamma}\rangle$ on $D$ for constant $\varphi$ shows a clear logarithmic increase with $D$ for $\chi_{2}^{2}$ distributed edge weights. In contrast, for the uniform distribution, the average horizontal extent increases with $\varphi$ as a power-law, see (b). The power-law behavior sets in for small $\varphi$ only for larger $D$. For $\varphi=4$ and $\varphi=5$ deviations from the straight lines in Figure 6 (b) are visible for $D<10^{5}$ and $D<10^{4}$, respectively. The average horizontal extent for constant $\varphi=2$ and $\varphi=3$ is not shown. We found that it does not converge to $\varphi$ in the range of matrix sizes $D$ we investigated.

For $\varphi\sim\log D$ the dependence of $\langle\tilde{\gamma}\rangle$ on $D$ is logarithmic for both the $\chi_{2}^{2}$ and the uniform distribution, as shown in Figure 6 (c) and (d).

In the remainder of this section, we will present analytic arguments similar to Section IV. We will explain the difference of the dependence of $\langle\tilde{\gamma}\rangle$ on $D$ for constant $\varphi$ between $\chi_{2}^{2}$ and uniform-like distributions. We show that $\left\langle\tilde{\gamma}\right\rangle\sim\log D$ for both distributions and $\varphi\sim\log D$.

V.2 Extent $\approx$ maximum of $\mathcal{J}$

By the Perron-Frobenius theorem the spectrum of $\mathcal{K}$ is confined to the ball centered around $\min_{i}K_{ii}<0$ with radius $r=|\min_{i}K_{ii}|$. Thus $2\max_{1\leq l\leq D}J_{ll}\geq|\operatorname{Re}\lambda|$ for all eigenvalues $\lambda$, so

[TABLE]

For symmetric generator matrices $\mathcal{K}$ we showed in Section IV.2 that

[TABLE]

and stated a result from [71] that for symmetric random generator matrices under mild conditions on the weights $K_{ij}$, $\max_{1\leq l\leq D}J_{ll}$ concentrates around the largest eigenvalue in magnitude, $\tilde{\gamma}$. This together with the upper bound by the Perron-Frobenius theorem Eq. (35) leads to our conjecture that the concentration of $\max_{1\leq l\leq D}J_{ll}$ around $\tilde{\gamma}$ in the symmetric case extends to the non-hermitian case as well

[TABLE]

A concentration result similar to the one in [71] for non-hermitian random generator matrices $M$ has to the best of our knowledge not appeared in the literature.

To quantify the deviation in Eq. (37) we introduce the relative error of $\tilde{\gamma}$ and $\max_{1\leq l\leq D}J_{ll}$

[TABLE]

In Figure 5 (c) and (d) we show the average relative error $\langle\delta\tilde{\gamma}\rangle$ as a function of $D$ and $\varphi$. If the edge weights are $\chi_{2}^{2}$ distributed then for $2\leq\varphi\leq 20$ and $10^{3}\leq D\leq 10^{5}$ the average relative error is smaller than $\approx 10^{-3}$ and decreases with increasing $D$. Thus Eq. (37) is a good approximation for large $D$ and $\varphi\ll D$ and the error appears negligible in the limit of large $D$. For uniformly distributed edge weights the average relative error $\langle\delta\tilde{\gamma}\rangle$ is smaller than $10^{-1}$ for $2\leq\varphi\leq 20$ and $10^{3}\leq D\leq 10^{5}$ and for $\varphi\geq 4$ decreases with $D$. For $2\leq\varphi\leq 3$, the error does not seem to decrease for increasing $D$. We conclude that Eq. (37) is an excellent approximation for large $D$ and $4\leq\varphi\ll D$.

V.3 Extreme value theory

Recall that the diagonal elements of $\mathcal{J}$ are iid random variables. Similar to the minimum extreme value statistics, the expected value of $\max_{1\leq l\leq D}J_{ll}$ is

[TABLE]

where we denoted the cdf of the diagonal elements $J_{ll}$ of $\mathcal{J}$ by $F$ and the pdf by $f=\frac{d}{dx}F$. A numerical calculation of the integral in Eq. (39) for $\chi_{2}^{2}$ distributed edge weights is shown in Figure 6 (a) and (c) and compared to the average horizontal extent $\langle\tilde{\gamma}\rangle$. The quantities agree excellently.

The remainder of this section is devoted to employing the Fisher-Tippet Gnedenko or extreme value theorem to $\max_{1\leq l\leq D}J_{ll}$ and thus analytically calculate the integral in Eq. (39).

V.3.1 Gamma distribution

Recall that if the edge weights are distributed according to a $\chi^{2}$ distribution then the diagonal elements of $\mathcal{J}$ are gamma distributed. The maximum of $D$ gamma distributed iid random variables $X_{l}$ converges in law to a standard Gumbel distribution $\operatorname{Gum}$ [71],

[TABLE]

with parameters $c$ and $d(D)$ given in table 2 for the gamma and $\chi^{2}$ distributions. The cdf of the Gumbel distribution is $x\to e^{-e^{-x}}$ with mean $\gamma$, where $\gamma$ denotes the Euler-Mascheroni constant, not to be confused with the horizontal extent $\tilde{\gamma}$. The assumption that the first moment converges and the convergence is fast enough together with Eq. (37) yields

[TABLE]

For constant $\varphi$ and increasing $D$ the dominant contribution of $d(D)$ is $2\log D$ for the $\chi^{2}$ distribution. Thus the increase is expected to be logarithmic. This is qualitatively consistent with numerical calculations of the average horizontal extent of random generator matrices $\mathcal{K}$ with $\chi_{2}^{2}$ distributed edge weights and constant $\varphi$ shown in Figure 6 (a). There $\tilde{\gamma}$ increases logarithmically with $D$. Quantitatively, the deviation between the average horizontal extent and the right-hand side of Eq. (41) is not small. The deviation decreases for increasing $D$ (not shown). We attribute the slow convergence to a sub-optimal choice of parameters $c$ and $d(D)$, as the right-hand side of Eq. (39) agrees perfectly with the numerically calculated $\langle\tilde{\gamma}\rangle$.

Let us assume that Eq. (41) is valid for $\varphi$ increasing logarithmically. Note that for the $\chi^{2}$ distribution, the rate parameter of the corresponding gamma distribution is linear in $\varphi$. Thus for large enough $\varphi$ by Stirling’s formula, the dominant term in Eq. (41) is logarithmic in $D$. Hence the average horizontal extent $\langle\tilde{\gamma}\rangle$ should increase logarithmically for $\varphi\sim\log D$. This is again qualitatively confirmed by numerical results shown in Figure 6 (c), where $\langle\tilde{\gamma}\rangle$ as a function of $D$ for $\varphi\sim\log D$ increases logarithmically with $D$.

V.3.2 Power-law tail distributions with bounded support

If the distribution of the diagonal of $\mathcal{J}$ has bounded support and the right tail decreases as a power-law, then we can reuse the extreme value result from Section IV.3. For a random variable $X$ with right support endpoint $x_{0}$ and cdf $F$ with power-law right tail, i.e.

[TABLE]

then $m_{D}=\max_{1\leq i\leq D}X_{l}$, where $X_{l}$ are iid copies of $X$, converges, properly renormalized, in law to a Weibull distribution

[TABLE]

Again, assuming that the first moment converges as well and the convergence is fast enough we get for edge weights distributed according to the standard uniform distribution,

[TABLE]

We find excellent numerical agreement of the right-hand side of Eq. (44) with the average horizontal extent $\langle\tilde{\gamma}\rangle$ for $\varphi\geq 4$. In Figure 6 (b) we show $\langle\tilde{\gamma}\rangle$ as a function of $D$ for fixed $\varphi$. The solid lines denote the right-hand side of Eq. (44). They agree perfectly for large enough $D$ and $\varphi\geq 4$. For $4\leq\varphi\lessapprox 6$ and small $D$ the agreement is still reasonable but deviations are clearly visible. Thus for fixed $\varphi\geq 4$ and increasing $D$, $\langle\tilde{\gamma}\rangle$ converges to $\varphi$ as a power in $D$ with exponent $-1/\varphi$,

[TABLE]

Numerically we find that Eq. (44) is valid for $\varphi$ increasing with $D$ logarithmically, see Figure 6 (d). There we show the average horizontal extent as a function of $D$ for $\varphi=\log D$. It increases logarithmically with $D$. The logarithmic increase can be justified analytically by extending Eq. (44) beyond constant $\varphi$. In the limit of large enough $\varphi$ we approximate $\Gamma(1+1/\varphi)\approx 1$ and by Stirling’s formula $(\varphi!)^{1/\varphi}\approx\varphi/e$ and get

[TABLE]

Thus in the limit of large $\varphi$ the average horizontal extent increases as $\sim\varphi\sim\log D$.

V.4 Summary

We showed numerically and analytically that the horizontal extent increases logarithmically for $\chi_{2}^{2}$ and uniformly distributed edge weights if $\varphi\sim\log D$. For constant $\varphi\geq 4$ and uniformly distributed edge weights the horizontal extent increases to $\varphi$ as $\sim\varphi-D^{-1/\varphi}$, while $\langle\tilde{\gamma}\rangle$ increases logarithmically for constant $\varphi$ and $\chi_{2}^{2}$ distributed edge weights.

The difference of the dependence of the average horizontal extent on $\varphi$ between the $\chi_{2}^{2}$ and uniform distribution goes back to the difference of the right tails. When edge weights are uniformly distributed the diagonal has bounded support and a power-law right tail, while it has unbounded support and an exponentially decaying right tail for $\chi_{2}^{2}$ distributed edge weights.

Similar to the spectral gap the limiting distribution of the horizontal extent is given by the limiting extreme value distribution of the diagonal elements of $\mathcal{K}$, under the assumption that the largest (in magnitude) diagonal of $\mathcal{K}$ is approximating $\tilde{\gamma}$ well enough. Thus the classification of the horizontal extent with respect to weight distributions reduces to the classifications of convergence in extreme value theory.

VI Complex Spacing ratios

So far we considered the marginal distribution of eigenvalues of sparse random generator matrices. But correlations between the eigenvalues are also of interest. Correlations between eigenvalues of real spectra are often quantified with the distribution of consecutive level spacings or their ratios. The latter avoids the need to unfold the corresponding spectrum [82, 83] and has been generalized to complex eigenvalues in the recent work [61]. The complex spacing ratio (CSR) of eigenvalue $\lambda$ of matrix $\mathcal{K}$ is defined as

[TABLE]

where $\lambda^{NN}$ and $\lambda^{NNN}$ denote the closest, by the Euclidean distance, and second closest eigenvalue (of $\mathcal{K}$) to $\lambda$, respectively. By definition, the density of CSRs is supported on the unit disk on the complex plane.

If eigenvalues $\lambda$ are uncorrelated, the CSR density is uniform. Eigenvalues of generic random matrix ensembles are typically correlated and feature mutual repulsion. This leads to vanishing CSR density at $z=0$ and $z=1$. According to [84], complex level spacings categorize random matrix ensembles in three universality classes. Generic random matrices fall into one of these classes according to their symmetries. The random generators considered in this paper have real entries so they obey the same symmetry as real Ginibre matrices (GinOE).

In Figure 7 we show the CSR densities of (a) GinOE members with Gaussian entries (b-d) and sparse random generators with $\chi_{2}^{2}$ distributed edge weights and $\varphi=1,2,3$. The densities are estimated from $100$ samples for $D=10^{4}$. We also checked that the obtained densities are independent of the weight distribution. As suggested in Ref. [61], we avoid eigenvalues close to the real line (by excluding all eigenvalues from the strip $\operatorname{Im}\lambda<10^{-14}$) when sampling CSR densities.

The CSR density of GinOE matrices shown on Figure 7 (a) exhibit typical depletion at $z=0$ and $z=1$. In Ref. [40], it was shown that the CSR density obtained for dense random Kolmogorov operators agrees well with the distribution shown in Figure 7 (a). The CSR density of sparse generators with sparsity $\varphi\geq 2$ (c,d) agrees remarkably well with the GinOE case.

The CSR density for $\varphi=1$ is anomalous, see Figure 7 (b). It has an extremely high density around $z=-1$ while being nearly flat on the rest of the unit disk. In this ultimate case, the operator can be presented as

[TABLE]

where ${\cal V}$ is a diagonal matrix (with elements distributed according to, e.g., $\chi_{2}^{2}$) and ${\cal P}$ is a circulant permutation matrix corresponding to a cyclic unit shift. The spectrum of ${\cal P}-{\mathchoice{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm 1\mskip-5.0mul}}$ lies on a circle of unit radius centered at $\lambda=-1$ and constitutes a set of equidistant roots of unity. This spectrum is slightly deformed and split into several loops by multiplication of ${\cal P}-{\mathchoice{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm 1\mskip-5.0mul}}$ with ${\cal V}$. Away from $\lambda=-1$, ${\cal V}$ dominates which results in the appearance of a real ’tail’; see Fig. 8.

In graph theory terms, such a sparse random graph fragments into a set of disjoint elementary cycle graphs. The independence of the spectra of these cycles leads to the flatness of the density away from $z=-1$, while the elementary cycle structure of the connected components is responsible for the CSR peak at $z=-1$.

To quantify the distance between CSR densities, we use the average length $\langle r\rangle$ and the average cosine of the angle $-\langle\cos\theta\rangle$ of spacing ratios, where $\langle\dots\rangle$ again denotes the average over the random matrix ensemble [61]. We numerically estimate $\langle r\rangle_{\operatorname{GinOE}}\approx 0.7379$ and $-\langle\cos\theta\rangle_{\operatorname{GinOE}}\approx 0.2347$ for 100 $10^{4}\times 10^{4}$-matrices. These agree well with $\langle r\rangle$ and $-\langle\cos\theta\rangle$ for $\varphi=2$ and $\varphi=3$, as shown in Table 3. We found similar results for $\varphi>3$ (not shown). In contrast, the corresponding quantities for $\varphi=1$ deviate substantially from $\langle r\rangle_{\operatorname{GinOE}}$ and $-\langle\cos\theta\rangle_{\operatorname{GinOE}}$, as also shown in Table 3. We conclude that, for $\varphi\geq 2$, correlations between eigenvalues of sparse random Kolmogorov operators agree with correlations of eigenvalues of GinOE matrices.

VII discussion

VII.1 Summary of results

Motivated by the incapability of dense random Kolmogorov operators to capture spectral features of model Markov processes, we introduced and analyzed an ensemble of sparse random Kolmogorov operators. We showed that, if the number of non-zero elements per column (and row) $\varphi$ increases with the matrix size $D$, the bulk of the spectrum is shifted away from the stationary eigenvalue [math] in the limit of large matrix size $D$. This is independent of the weight distribution, i.e. of the distribution of the nonzero matrix elements.

In contrast, the spectral edges depend on the tails of the weight distribution. The tails of the weight distribution determine, together with $\varphi$, the tails of diagonal elements of generator matrices. We numerically showed that the spectral edges are well approximated by the extremes of the diagonal elements. From extreme value theory it follows that for diagonal distributions with power-law left tails (this includes among others edge weights being uniform, exponential, $\chi^{2}$, gamma or beta distributed), the average spectral gap decreases as a power-law in $D$ for fixed $\varphi$, is constant for $\varphi\sim\log D$ and increases, whenever $\varphi$ increases with $D$ substantially faster than $\log D$.

A similar approach was used to calculate the horizontal extent, given by the eigenvalue with the smallest real part. We linked the horizontal extent to the largest diagonal element (in magnitude) of the generator matrix and used extreme value theory to calculate the latter.

Finally, we showed that complex spacing ratio distributions of generator matrices with $\varphi\geq 2$ follow the distribution typical of Ginibre’s Orthogonal Ensemble, while there is a strong anomaly for $\varphi=1$.

VII.2 Open questions

(1) We have introduced sparsity to model $\mathcal{K}$-generators of physical Markov processes, and have used the sparsity to tune spectral features of the generators. There are other ways of providing random matrices with a structure that models physical constraints (e.g., locality). E.g., one could consider banded matrices [85, 86, 87, 88, 89, 90, 91, 92, 93, 94] or matrices with decaying off-diagonal terms [95, 96, 93] or temperature based models [97]. These are alternate routes to tuning spectral features. To the best of our knowledge, generators of CTMCs with such structures have not yet been considered.

(2) The application of extreme value theory to find the limiting distribution of the spectral edges relied on the observation that the spectral edges are well approximated by the minimum and the maximum of the diagonal of the generator matrix. By the Courant-Fisher theorem, the extremes of the diagonal are upper and lower bounds, respectively, for symmetric generators. In this case, a concentration of the largest eigenvalue in magnitude around the maximum of the diagonal was shown in [71]. An analytical treatment of general non-symmetric generators and the spectral gap is to the best of our knowledge not known. We hope that our results motivate a rigorous investigation of the connection between the spectral edges and the diagonal of the generator matrix.

(3) Generators of CTMCs have real entries and thus their eigenvalues are real or come in complex conjugate pairs. In the investigation of correlations between eigenvalues, we left out real eigenvalues. The appearance of a large number of real eigenvalues in the spectrum of non-Hermitian matrices is a phenomenon of wide interest [98, 99, 100, 101, 102, 103, 55, 104, 105, 106, 107]. For real Ginibre matrices, the average number of real eigenvalues is $\sim D^{-1/2}$ [99, 100, 101] while for dense generators of CTMCs, it is substantially larger [55]. We observed that the fraction of real eigenvalues is larger for small $\varphi$ and smaller for larger $\varphi$ (not presented). Understanding of the functional dependence of the number of real eigenvalues for sparse CTMC generators is an interesting problem.

(4) We focused on the location and extent of the bulk spectrum as well as the spectral edges. One could inquire about the evolution of other features of the spectral distribution as a function of sparsity, e.g., about the envelope of the spectral distribution. In [40] the spectral density of dense random CTMCs was described by the convolution of two asymptotically free matrices, leading to the prominent spindle shape. Free probability arguments break down for sparse random CTMCs. Analytical tools which have been employed to calculate the spectral density of sparse, random matrices include replica tricks [108, 109, 110, 111, 112], single defect and effective medium approximations [113, 114, 115], supersymmetry-based techniques [116, 117] and the cavity approach [118, 119, 120, 112]. Spectral properties of symmetric, sparse, random CTMCs have been investigated with the cavity method [121, 122, 123] and with supersymmetric approaches [117]. Investigations of the spectral density of non-symmetric sparse, random Kolmogorov operators with the above methods might be an interesting objective.

(5) We have considered sparse generators of CTMCs based on strongly connected, sparse random graphs. It is an open question whether our results can be generalized to other sparse graph ensembles. One potential avenue to explore are directed Erdös-Renyi (dER) graphs.

In dER graphs, the probability of an edge connecting any two vertices is $0<p\leq 1$. For a dER graph to be strongly connected with a high probability, the value of $p$ must exceed $\sim\log D/D$ [124, 125]. As a result, the average degree of the vertices must increase logarithmically with $D$ to ensure strong connectivity. Consequently, the range of constant average vertex degree and increasing vertex number $D$ is excluded.

Nonetheless, modifying the dER graph by enforcing a minimum (in- and out-) degree $\geq 2$ guarantees strong connectivity with high probability [79]. Exploring the spectral properties of CTMC generators based on dER graphs may represent a promising next step towards generalizing our results.

(6) Finally, there is an interesting question: What could ’sparsity’ mean in the quantum limit? Namely, what is ’sparsity’ for Lindblad operators?

Here we start from the genetic link which allows to obtain a generator of a classical (quantum) Markovian evolution as a properly normalized non-probability (trace) preserving stochastic map (channel). Eq. (3), where $\cal{M}$ is a no-probability-preserving map and $\cal{J}$ takes care of normalization, illustrates this link in the case of Kolmogorov operators. It seems to be intuitive that sparsity $\varphi$ of a stochastic matrix can be associated with rank $r$ of a quantum channel [126].

Thus, in the ultimate limit $\varphi=r=1$, quantum versions of stochastic maps – that are permutations now – are rank-one channels – that are unitaries. It is tempting to extend this analogy beyond the limit $\varphi=r=1$ and state that mixed-unitary channels (convex combinations of unitaries) are quantum versions of bistochastic matrices (that are, according to Birkhoff, convex combinations of permutations [127]).

However, there is also notion of double stochastic (or ”unital”) channels [126]. The two classes, mixed unitaries and bistochastic channels, are not identical: there are double stochastic channels that are outside of the convex hull of unitaries [128]. What class to associate with classical bistochastic maps is then ’a matter of taste’  [129]. To resolve the dichotomy, one could rely on the concept of super-decoherence [40] and state that all channels which have classical bistochastic matrices as their fully decohered versions, are quantum analogues of the matrices. In this case the broader class of double stochastic channels is chosen [129].

The superdecoherence-based reasoning can also be applied to generators. In this case the unitary (Hamiltonian) part of a Lindblad operator does not play any role since it vanishes in the limit of complete decoherence [40] and the quantum ’sparsity’ is defined by the rank of the dissipative part (the minimal number of jump operators).

Interestingly, Lindblad operators of ultra-low rank $r=1$ were considered in Refs. [52] and [130]. Features similar to ones we detected for ultra-sparse Kolmogorov operators were reported (e.g., the spectral gap is defined by a real-valued outlier with position independent of $D$).

Acknowledgements.

The authors thank Gernot Akemann, Alexander van Werden, Ioana Dumitriu, Tomaz̆ Prosen, Karol Życzkowski, and Wojciech Tarnowski for helpful discussions and comments. This research is supported by the Deutsche Forschungsgemeinschaft through SFB 1143 (project-id 247310070) [GN and MH] and Research Council of Norway, project IKTPLUSS-IKT og digital innovasjon - 333979 (as a part of the ERA-NET project ”DQUANT: A Dissipative Quantum Chaos perspective on Near-Term Quantum Computing”) [SD].

Appendix A Models in Figure 1

In the following, we denote spin creation and annihilation operators as $\sigma^{+}$ and $\sigma^{-}$, respectively. We split the generator matrix $\mathcal{K}$ into an off-diagonal matrix $\mathcal{M}$ and a diagonal matrix $\mathcal{J}$ such that $\mathcal{K}=\mathcal{M}-\mathcal{J}$ and the diagonal entries of $\mathcal{J}$ are the sums of the columns of $\mathcal{M}$.

In Figure 1 (b) we show the TASEP on a ring with $L=12$ sites and staggered hopping amplitudes. The $\mathcal{M}$ matrix is given by [61]

[TABLE]

where $p_{j}=1$ if $j$ is even and $p_{j}=0.2$ if $j$ is odd.

In Figure 1 (c) we show the ASEP on a chain of length $L=12$ with open boundary conditions and next nearest neighbor hopping. The $\mathcal{M}$ matrix is given by

[TABLE]

In Figure 1 (d) we show the spectrum of a single particle hopping on a $65\times 65$ grid with periodic boundary conditions and random hopping amplitudes. The $\mathcal{M}$ matrix is given by

[TABLE]

where $\langle\dots\rangle$ denotes summation over nearest neighbors and $p_{(i,j)\to(i^{\prime},j^{\prime})}$ are randomly uniformly chosen between 0 and 1 under the constraint that $p_{(i,j)\to(i^{\prime},j^{\prime})}=1-p_{(i^{\prime},j^{\prime})\to(i,j)}$. This diffusion model can of course be extended to many particles, but we choose to show the single-particle sector here.

In Figure 1 (e) we show the spectrum of a contact process [28] on a chain with $L=12$ sites and open boundary conditions. The master equation is generated by $-H$, where $H$ is given by

[TABLE]

and

[TABLE]

Finally, in Figure 1 (f) we show the spectrum of the generator matrix $\mathcal{K}$ of a gene transcription model taken from [20]. The following master equations model the accumulation and release of mechanical strain of DNA during transcription. The parameters chosen for the spectral data in Figure 1 (f) are the mRNA transcription rate $r=2$ and decay rate $\lambda=0.05$, the maximum number of transcripts until no further strain can be put on the DNA $m_{c}=10$, the relaxation rate of the DNA string $g=0.05$ and a maximum number of transcription events $m_{\max}=400$ to make the generator matrix $M$ finite. By $m$ we denote the number of current transcripts and by $\alpha$ the number of transcripts made since the last relaxation event. Then for $0\leq m\leq m_{\max}$ and $1\leq\alpha\leq m_{c}-1$ the master equation reads

[TABLE]

while for $\alpha=0$ we have

[TABLE]

and for $\alpha=m_{c}$

[TABLE]

Appendix B Sampling of sparse random generators

In order to obtain a sparse random generator matrix, our approach involves first sampling a random directed graph with $D$ vertices and both in- and out-vertex degrees of $\varphi$. Subsequently, the non-zero elements of the corresponding adjacency matrix are sampled from a common positive distribution. This procedure results in the off-diagonal matrix $\mathcal{M}$. The random Markov generator matrix is then constructed as $\mathcal{K}=\mathcal{M}-\mathcal{J}$, where $\mathcal{J}$ is a diagonal matrix with diagonal elements equal to the sums of the columns of $\mathcal{M}$.

The random directed graph is generated by iteratively connecting each vertex to $\varphi$ other vertices, while rejecting edges if the corresponding vertex already has $\varphi$ incoming edges. For the final vertices, it may not be feasible to connect to other vertices without violating the constraint of $\varphi$ incoming edges for each vertex. In such cases, the entire process is restarted. To mitigate the risk of restarting the procedure, we reduce the probability of connecting to a vertex that already has a high degree. Following this approach, we find that we rarely need to restart the algorithm for the matrix sizes and vertex degrees $\varphi$ examined in this study.

To compute the eigenvalues of the Markov matrices, we utilize an exact diagonalization method, while the Arnoldi method is employed to calculate the spectral gap. We deem an eigenvalue to have converged once the norm of the residuals of the Schur vectors is less than $10^{-12}$.

Appendix C Analytical results for the bulk spectrum

In this section, we will derive the analytical results of the estimated mean $\mu(\lambda)$ in Eq. (7) and the estimated pseudo-variance in Eq. (III.2) in the main text and show that $\frac{1}{D}\sum_{j=1}\lambda_{j}$ concentrates around its average $\langle\dots\rangle$.

Denote by $\iota$ the function $\iota:\{1,\dots,\varphi\}\times\{1,\dots,D\}\to\{1,\dots,D\}^{2}$ with $\iota(l,j)=(i,j)$ where $i$ is the $l$th non-zero index in column $j$ in $\mathcal{M}$. Note that $\iota(l,j)=(i,j)$ implies $i\neq j$ and $l\to\iota(l,j)$ is injective for fixed $j$. Further, let in this appendix the location of the bulk be denoted as

[TABLE]

and the pseudo-variance as

[TABLE]

Here we explicitly do not include the averaging over the matrix ensemble $\langle\dots\rangle$ in contrast to the main text.

C.1 Location

The average value with respect to $\langle\dots\rangle$ of the location $\mu(\lambda)$ can then be computed as

[TABLE]

where we used that $\left\langle K_{\iota(l,j)}\right\rangle=-\mu_{0}$. This is Eq. (7) in the main text. Similar,

[TABLE]

Although the off-diagonal elements of $\mathcal{K}$ are weakly dependent because of the constraint that the number of non-zero elements per row and column has to equal $\varphi$, the non-zero elements $K_{\iota(l,j)}$ are independent. Hence, $\left\langle K_{\iota(l_{1},j)}K_{\iota(l_{2},j)}\right\rangle=\left\langle K_{\iota(l_{1},j)}\right\rangle\left\langle K_{\iota(l_{2},j)}\right\rangle$ and $\left\langle K_{\iota(l_{1},j_{1})}\right\rangle\left\langle K_{\iota(l_{2},j_{2})}\right\rangle$, so

[TABLE]

where we used that the second moment $\left\langle K_{\iota(l,j)}^{2}\right\rangle$ equals $\sigma_{0}^{2}+\mu_{0}^{2}$. This implies that

[TABLE]

The right-hand side vanishes for increasing $D$ and $\varphi$ growing slower with $D$ than linear. Relatively to $\left\langle\mu(\lambda)\right\rangle$ the typical deviation of $\mu(\lambda)$ from its average value always vanishes for either increasing $D$ or $\varphi$, as

[TABLE]

C.2 Complex pseudo-variance

The first term in the averaged pseudo-variance given by Eq. (C) can be calculated as

[TABLE]

We proceed with $\sum_{i=1}^{D}\langle K_{ii}^{2}\rangle$ in Eq. (C.2) and get

[TABLE]

The former sum in Eq. (C.2) is given by

[TABLE]

where again we used that $\langle K_{\iota(l,i)}^{2}\rangle=\sigma_{0}^{2}+\mu_{0}^{2}$, while the latter sum in Eq. (C.2) is

[TABLE]

Combining Eq. (60) and Eq. (C.2) we get

[TABLE]

Now, we are left with calculating $\sum_{i\neq j}\langle K_{ij}K_{ji}\rangle$, the second term in Eq. (C.2),

[TABLE]

where the $\overline{\iota}$ denotes swapping the first and second component. Note that $K_{\overline{\iota(l,i)}}$ is not necessarily a non-zero entry of $\mathcal{K}$, hence $K_{\overline{\iota(l,i)}}$ and $K_{\iota(l,i)}$ depend weakly on each other. In the large $D$ limit we can assume that the dependence is sufficiently weak and we treat $K_{\overline{\iota(l,i)}}$ and $K_{\iota(l,i)}$ as independent, thus $\left\langle K_{\overline{\iota(l,i)}}K_{\iota(l,i)}\right\rangle=\mu_{0}\left\langle K_{\overline{\iota(l,i)}}\right\rangle$. By the assumed independence the mean of every entry in the $i$th row, except the diagonal, is $\left\langle K_{\overline{\iota(l,i)}}\right\rangle=\frac{\varphi}{D}\mu_{0}$. Hence,

[TABLE]

Collecting the above results we arrive at

[TABLE]

The second term of the averaged pseudo-variance in Eq. (C) has been calculated in the previous subsection,

[TABLE]

Finally, we can evaluate

[TABLE]

which is Eq. (III.2) in the main text.

Appendix D Bound of spectral gap for symmetric $M$

In this section, we give the proof of Eq. (19). Let $\mathcal{K}=\mathcal{M}-\mathcal{J}$ be a symmetric generator matrix. By Eq. (17) we have to show that $v^{t}\mathcal{K}v\leq\min_{1\leq l\leq D}J_{ll}+O\left(D^{-1}\right)$ for the vector $v$ given

[TABLE]

where $1\leq l\leq D$ is arbitrary. It is easy to see that $|v|=1$ and $v\perp v_{1}$. So we proceed with

[TABLE]

Note that any summand in Eq. (D) where either $i=j=l$ or $i\neq l$ and $j\neq l$ is zero. Inserting the definition of $v$ we get

[TABLE]

After collecting all the prefactors in Eq. (D) the spectral gap is upper-bounded by

[TABLE]

where we denote $\tilde{J}_{ll}=\sum_{i\neq l}M_{il}$. As the number of non-zero elements of $\mathcal{M}$ in every row and column is the same, the distribution of $J_{ll}$ and $\tilde{J}_{ll}$ coincide. In the limit of large $D$, $J_{ll}$ and $\tilde{J}_{ll}$ are independent. Thus we can approximate $\gamma_{*}\leq J_{ll}+O\left(D^{-1}\right)$ at least for $\varphi\ll D$. As the index $l$ was chosen arbitrarily we get

[TABLE]

which is Eq. (19) in the main text.